1 Introduction and Problem Definition
Computer simulations have been widely and successfully used to understand and predict physical behaviour in biological, medical, technical, and many other applications. A central entity in simulation technology is the computational domain on which the solution is sought numerically. Commonly, the computational domain coincides with the spatial extent of the physical object under consideration. In many applications, this leads to a time-dependent spatial domain and—thinking for example of valves or bearings—there are also various applications where changes its topology over time.
Several approaches have been developed to handle the time variant topology of . Our focus is on mesh-based approaches, which can be broadly organised in three groups, namely, boundary conforming, quasi boundary conforming with residual gap, and non-boundary conforming approaches. Non-boundary conforming (non-interface fitted) approaches find a way to shift the complication of the topology change from the discretisation to the integration. In the field of fluid-structure contact interaction (FSCI) a number of methods has recently been developed. Ager et al.  employ a fixed Eulerian mesh as fluid discretization and handle boundary motion and topology changes of the physical fluid domain with an approach based on CutFEM . This includes a suitable method for the numerical integration on boundary intersected elements , a ghost penalty stabilization to overcome issues of very small cuts, and a Nitsche-based weak imposition of fluid boundary (coupling) conditions. Kapada et al. propose a computational framework based on Cartesian hierarchical B-spline grids for the fluid discretization . Therein, integration over intersected elements is handled via sub-triangulation or uniform subdivision of the rectangular grid; a ghost penalty stabilization and Nitsche’s method for fluid boundary conditions is used. Both aforementioned formulations bear similarities with the immersogeometric variational framework proposed previously by Kamensky et al. .
The second class of approaches use a boundary conforming discretization, but avoid the actual topology change of the fluid domain with a small residual gap. Ensuring the residual gap with a penalty force or a displacement restriction in the mesh motion, space-time  or ALE methods  are employed to handle the deforming fluid domain.
Finally, boundary conforming discretizations for spatial computational domains with topology changes can be obtained by reformulating the problem as space-time problem with contiguous computational domain. For two-dimensional spatial computational domains, this approach has been implemented with unstructured finite volume  and finite element meshes .
In this work, we will follow the boundary conforming approach and choose the space-time domain as computational domain. The concept is visualised in Figure 1 with the spatial domain collapsed onto the horizontal axis. Between and , the spatial domain splits into two unconnected parts, however, the computational space-time domain is contiguous. A boundary conforming discretization with pentatope elements leads to a simplex space-time (SST) mesh, which is fully unstructured in space and time.
To generate such four-dimensional simplex meshes, several similar approaches have been developed recently. The robust meshing strategy proposed by Behr 
extrudes a tetrahedral mesh resulting in a four-dimensional mesh of tensor product elements. These hyperprisms are split into pentatopes with an element-wise Delauney triangulation of the perturbed nodal coordinates. In a further development, Karabelas and Neumüller replace the element-wise Delauney triangulation with a predefined decomposition of each hyperprism requiring a consistently numbered tetrahedral mesh. A third approach of Wang similarly employs a global node indexing scheme and extends it with a node insertion procedure to support local mesh density operations . The above mentioned strategies have in common that the four-dimensional mesh is based on an extruded tetrahedral mesh. An alternative approach to generate high-quality pentatope meshes could be based on Coxeter triangulations, with the currently open issue of generating boundary conforming meshes .
In most cases, SST meshes are employed to facilitate adaptive mesh refinement in space and time. Suitable pentatope mesh refinement procedures have been explored by Neumüller and Steinbach  and Grande . Anisotropic four-dimensional mesh adaptation is pioneered by Caplan et al [16, 17] and successfully employed in the solution of the advection diffusion equation . Further, recent application examples of four-dimensional SST meshes from the field of mathematics deal with parabolic evolution problems [19, 20] or a broader class of transient PDEs recast as constrained first-order system . In the field of computational engineering, adaptive temporal refinement of pentatope meshes is used for two-phase flow simulations —also combined with complex material laws such as the Carreau-Yasuda-WLF model  or the -rheology —as well as gas flow simulations in the piston ring pack of internal combustion engines .
In this work, we generate pentatope finite element meshes for spatial domains with time variant topology. Therefore, we combine the extrusion based approach by Behr  with a four-dimensional extension of the elastic mesh update method (EMUM). Originally proposed as elastic grid approach by Lynch , EMUM was refined and employed as automatic mesh moving scheme for the deforming spatial domain / stabilized space-time finite element formulation . Since then, it has been widely used and become a standard technique to handle moving domains in fluid-structure interaction (FSI) simulations [27, 28, 29]. Also, very recent FSI simulations rely on it [30, 31, 32]. Furthermore, EMUM has been used in the context of free-surface flows (for example in ) and to update finite element meshes according to prescribed boundary displacements (for example in ). The four-dimensional extension of EMUM (4DEMUM) presented in Section 2 allows us to obtain boundary conforming pentatope meshes of complex geometries. The geometry may have holes, and does not have to have a tensor product shape in any dimension. This means that the pentatope mesh can account for a time variant topology of the spatial domain. However, the mesh generation method is still limited to cases where the four-dimensional geometry can be obtained by extrusion of a (complex) three-dimensional geometry with subsequent elastic deformation (in the sense of 4DEMUM).
Now, we proceed as follows. In Section 2, we describe the strong and weak form for a finite element implementation of 4DEMUM. Next, Section 3 presents the embedding of 4DEMUM into the workflow of a four-dimensional space-time finite element simulation, followed by Section 4 which aims at experimentally validating the method, as well as demonstrating the particular potential of simulations on SST meshes for domains with time variant topology. Finally, Section 5 contains concluding remarks and a brief outlook.
2 Four-Dimensional Elastic Mesh Update Method
The basic concept of the elastic mesh update method is to update the node positions according to the deformations of an elastic solid. As outlined above, the method has been widely used to adjust two- or three-dimensional meshes to deforming boundaries. In the following, we extend the method to four-dimensional meshes. Therefore, we consider a virtual four-dimensional linear elastic solid occupying the region .
Figure (a)a displays two dimensions of an example of . Let’s assume the domain boundary consists of a Neumann part and a Dirichlet part
. For each degree of freedom the boundary is split such thatand . We can now formulate the four-dimensional elastostatic problem
with Dirichlet boundary conditions prescribed on . Where no Dirichlet values are prescribed, a homogeneous Neumann boundary condition is assumed.
The assumption of a homogeneous, isotropic, linear elastic material behaviour leads to a constitutive equation
where is an equivalent of the Cauchy stress tensor and and are the Lamé parameters. In this work, we choose . Further,
denotes the identity matrix andis the sum of the main diagonal entries of the strain tensor.
The solid’s deformation is measured with a linear strain model
where are displacements in the four spatial directions. The displacements are observed at the nodes of the extruded mesh and the nodal positions are described by their coordinates . Accordingly, the nabla operator collects the partial derivatives.
For an admissible test function space with functions that vanish on , and a suitable trial function space that respects the Dirichlet boundary conditions, the discretised weak form of Equation (1) can be stated as: Find , such that for all
For details on pentatope Lagrange finite elements the reader is referred to the previous publications [9, 10]. In short, the weak form above results in a linear equation system that is solved for the displacements . Finally, the node positions with coordinates are obtained by adding the displacements to the extruded mesh coordinates:
In this addition, we identify with . Two dimensions of a resulting mesh on the space-time domain are shown in Figure (b)b. The displacements in all four dimensions are coupled, so interior mesh nodes can move in any direction, even if non-zero boundary displacements are applied only to specific ones.
3 Simulation Workflow
An overview of the workflow of a four-dimensional SST finite element simulation is given in Figure 5. The arrows correspond to seven steps that are taken to complete a specific subtask in the workflow. When suitable, the arrows are labeled with the software employed to complete the task. The first four steps are considered pre-processing and produce the SST meshes covering the considered physical space-time domain . In detail, they are: First, identify a projection , such that complex features of are included in . In our work, this step is performed manually. For an example of such a projection result see Figure 39. As second step, generate a tetrahedral mesh covering , for example with GMSH . Then, according to Behr , extrude the mesh covering over an interval resulting in a mesh of tensor-product elements and apply an element-wise Delauney triangulation to generate a pentatope mesh covering (third step). The fourth steps is to use 4DEMUM as described in Section 2 to deform the pentatope mesh to cover . Note that the time dimension can be one of the initial three dimensions, and extrusion can be used to generate one of the spatial dimensions of the mesh.
With a boundary conforming discretization at hand, the fifth step is to perform the space-time simulation. Details of the simulation depend on the application example and are given in Section 4
. In any case, the solution of the finite problem is obtained on the unstructured space-time mesh. To visualise the results over time, a series of tetrahedral meshes covering the spatial domain at given time instances is generated. The node positions of these meshes are passed as query points to an efficient interpolation tool. The tool identifies the pentatope of the space-time mesh which contains the query point and performs a linear barycentric interpolation on the element (step 6). This interpolation is in accordance with the finite element approximation of the simulation. For domains with convex boundaries some of the query points can be located slightly outside of the space-time mesh. In such cases, the solution from the element with the closest center is extrapolated. As final step 7, the data on the tetrahedral meshes can be easily visualised with available tools such as ParaView . An alternative visualisation approach described by Karabelas and Neumüller 
is based on the element-wise intersection of the pentatope mesh with a hyperplane.
4 Application Examples
For the physical simulations in the following application examples, we distinguish between space and time coordinates. The spatial nabla operator collects the derivatives with respect to the spatial coordinates. In the following two examples, we consider the compressible Navier-Stokes equations
Here , , , , , and
are density, velocity vector, pressure, viscous stress tensor, total energy per unit mass, and heat flux vector, respectively. For convenience, the conservation variables are collected in the vector. We consider an ideal, calorically perfect gas, with the specific gas constant , the ratio of specific heats , and a Prandtl number of . With these assumptions, we can perform a change of variables from the conservation variables to the pressure-primitive variables . The latter are used as primary unknowns and the governing equations above are formulated as generalized advective-diffusive system
Therein, partial time derivatives are denoted with , partial derivatives in each of the spatial directions are denoted with , and Einstein summation convention applies to repeated indices. Details of the generalized advection matrices and diffusion matrices can be found in [9, 38].
To find a numerical solution to the compressible Navier-Stokes equations with space-time finite elements, a discretisation of the physical space-time domain is required. An example of such a discretisation can be seen in Figure 6. In practise, a procedure is required to numerically integrate over , over the spatial computational domain at the initial time , and over which is the temporal evolution of the spatial domain boundary . For each degree of freedom, the space-time boundary consist of non-overlapping Dirichlet and Neumann parts .
Based on the domains introduced above and for admissible test and trial function spaces and , the weak form of Equation (9) can be stated as: For given initial conditions , find , such that for all
The first integral collects terms of the residual that are multiplied with the test function vector and integrated over . In the second integral, the derivatives with respect to the spatial coordinates are shifted to the test function vector, leading to the third integral over the Neumann part of the space-time boundary . Therein, the boundary normal fluxes are evaluated with the outwards pointing surface normal as
The initial condition is weakly enforced with the integral over . The weak form is completed with a SUPG operator to overcome instabilities of the pure Galerkin formulation which occur in convection-dominated flow simulations. For more details on the space-time finite element scheme for the compressible Navier-Stokes equations, the reader is referred to .
4.1 Transient Gas Flow through Valve
In the following, we present a transient three-dimensional gas flow simulation through a valve. The fluid domain consist of a square channel with sides in the --plane and a length of in -direction. In the course of the simulation, the position of the rounded valve member determines the topology of the spatial computational domain (see Figure 19). During the valve cycle, the valve member is first lowered and later lifted again. This leads to a splitting of the fluid domain at and a reconnection at (see Figure (a)a). Additionally, the exit cross-section of the channel deforms from the initial square geometry into a rectangular cross-section with half the size. Detailed geometry specifications and boundary conditions are collected in Figure 9. Along the solid walls, no-slip boundary conditions are enforced, i.e., the velocity is set to zero and a wall temperature is prescribed. On the left open boundary, the pressure and temperature are given; on the right open boundary a pressure value is set. Prescribed temperature and pressure values are summarised in Table LABEL:tab:valveBc. The gas viscosity is modelled using Sutherland’s relation
with , K, K.
Following the workflow outlined in Section 3, we start with a projection of four-dimensional space-time geometry of the valve into the ---hyperplane. The resulting three-dimensional ---domain corresponds to the temporal evolution of the side view shown in Figure (a)a and includes the topology change of the flow domain. We discretise the ---domain with tetrahedral elements. This is sufficient to resolve a transient two-dimensional flow field as documented in . In the next step, the mesh is extruded in z-direction and triangulated with pentatopes, such that both sides of the square channel cross-section are split into 20 line elements. To account for the closing exit, we identify with and apply 4DEMUM with
prescribed as Dirichlet boundary condition on the domain boundary. Therein, denotes the Heaviside function. After the mesh update, the fourth dimension of the resulting mesh is interpreted as time and the SST flow simulation is performed. The simulation took 30 minutes of wall clock time on 240 cores using a distributed memory parallelisation based on MPI.
The results displayed in Figure 19 show the pressure, velocity, and temperature distribution in the closed, half open, and fully opened valve. In Figure (h)h, it can be clearly seen that the reduced cross-section of the closing exit accelerates the flow to the maximum velocity of . The influence of the flow on the temperature distribution is shown in Figure (i)i. In summary, we observe a laminar flow with a maximum Reynolds number and Mach number at . The density variations between and underline the necessity to consider compressibility effects.
4.2 Flow Inspired by Clamped Artery
As further application example, we simulate a transient flow through an artery which is temporarily sealed by a clamp and reopened. The considered artery section is long and has an approximately circular cross-section with a diameter of . The clamp centre is located at . Regarding the fluid domain, we assume that the artery volume is displaced by the clamp and returns to the initial approximately circular shape as the clamp is removed. See Figure 26 for several views on the problem geometry in initial and clamped configuration. As time frame for the closing and opening, we choose the duration of one cardiac cycle approximated by . For the first , the artery is in its initial shape. Over the next , the clamp is applied and seals the artery from until . From until , the artery is reopened and and for the last it is again in its initial shape.
The flow through the deforming domain is driven by a pressure gradient of 40 mmHg, which corresponds to the difference between the minimum and maximum aortic pressure during a cardiac cycle. On the left open boundary, a pressure is prescribed. On the right open boundary, a pressure is prescribed. We assume that the fluid enters the domain at normal body temperature . The colder clamp leads to a temperature variation of along the artery wall. The wall temperature is denoted by . A summary of the corresponding boundary conditions is given in Table LABEL:tab:arteryBc. Regarding the velocity degrees of freedom, we strongly enforce the no-slip boundary condition on the arterial wall and set the tangential velocity components on the open boundaries to zero. The constant fluid viscosity in this test case, poise = , is inspired by blood. However, the gas constant, ratio of specific heats and Prandtl number are chosen as in the previous example, leading to a fictitious fluid with a density of roughly . More accurate models consider blood as an isothermal generalised Newtonian fluid  or even include shear-thinning, viscoelasticity and thixotropy .
4.2.1 Steady Simulation
To validate our formulation for these settings, we perform – in a first step – a flow simulation through a straight circular duct. The results are collected in Figure 29. The velocity distribution is visualised with glyphs showing a paraboloid in qualitative agreement with the Poiseuille paraboloid obtained for an incompressible Hagen-Poiseuille flow. Note that neither on the inflow nor on the outflow boundary the velocity component in normal direction is prescribed. For the pressure driven incompressible flow through a straight circular pipe, there is the well-known analytical solution for the velocity component in axial direction , which reads
for the considered configuration. Figure (b)b shows that the compressible flow solutions are flatter than . The simulation is performed on a coarse, medium, and fine mesh, with , , and tetrahedral elements, respectively. The slight variation of the centerline axial velocity— on the coarse mesh (3.1% smaller than fine), on the medium mesh (0.8% smaller than fine), and on the fine mesh—indicates that the coarse mesh resolution is sufficient to obtain a solution within the range of engineering accuracy.
4.2.2 Transient Simulation
As second step, we perform a transient simulation on a straight duct with approximately circular cross-section. Starting point for the SST mesh generation is an unstructured tetrahedral ---Mesh. The mesh has a temporal resolution comparable to the one shown in Figure 39, but the influence of the clamp is excluded for now. Next, the tetrahedral mesh is extruded in -direction, such that the cross-section in the --plane forms a square. In -direction, 14 nodes are added during the extrusion, such that the spatial resolution of the resulting mesh is comparable to the coarse pipe mesh discussed above. To obtain the approximately circular cross-section, we perform the 4DEMUM with the space-time coordinates identified with . Further, the extruded mesh is shifted and scaled, such that . On the pentatope mesh boundary , we prescribe the displacements
which map the square cross-section in the --plane on an approximately circular shape (see Figure (b)b). In a final step before the transient simulation, the mesh is shifted and rescaled such that . The shifts are performed to allow for relatively simple expressions in the boundary conditions of 4DEMUM, and at the same time allow the fluid simulation to start at .
In this second simulation, a transient feature is introduced by ramping up the inflow pressure from the initial value to over the first . For the second , the pressure value is kept constant as shown in Figure (a)a. The flow velocity at the center of the inflow (Figure (b)b) closely follows the temporal evolution of the pressure value. Note that the computed flow velocity is slightly larger than zero at . This is in line with our numerical formulation, which weakly enforces the initial condition (compare Equation (10)). However, the transient nature of the problem is properly captured in the computed flow field.
With this test case, we also want to explain our choice of an approximately circular cross-section. A perfectly circular cross-section introduces an expected complication in 4DEMUM, i.e., the elements formerly in the corners of the square cross-section attain very large dihedral angles and eventually lead to a tangled mesh. We circumvent this problem by introducing the prefactor in Equation (15) to obtain a mesh with an approximately circular cross section in the --plane.
The “missing 10%“ towards the perfect circular cross-section have hardly any influence on the flow field as shown in Figure 35. A more quantitative comparison is presented in Figure 38. The parabolic velocity profile in radial direction (Figure (a)a) as well as the linear pressure decay along the pipe axis (Figure (b)b) are obtained independently of the approximation of the circular cross-section. Therefore, we consider an approximately circular cross-section of the artery in the following.
4.2.3 Transient Simulation with Topology Change
As third and final step, we consider the transient simulation with topology change of the spatial computational domain. We now include the topology change caused by the clamp in the ---mesh shown in Figure 39. The further pentatope mesh generation steps of extrusion, connectivity generation, and elastic deformation are performed as in the previous example; the boundary displacements are given in Equation (15).
The subsequent finite element flow simulation was performed in on 240 cores using a distributed memory parallelisation based on MPI. It took 16 minutes of wall clock time. Figure 46 presents the simulation results at , , and . For the fully clamped case, we obtain two separate domains and negligible flow velocities (Figure (a)a). In the absence of flow, the temperature distribution in the fluid is an interpolation of the temperature prescribed on the domain boundaries (Figure (b)b). Note that on the right most boundary () no temperature boundary condition is prescribed because this part is an outflow boundary during most of the simulation. However, from until , back-flow across this boundary introduces a small disturbance in the temperature field. When the artery is reopened to roughly half of the total diameter, a strong pressure gradient across the clamp region accelerates the flow in this area (Figure (c)c). Further at , the temperature distribution is strongly influenced by the flow field as well as the cooler clamp (Figure (d)d). In the open configuration (Figure (e)e and (f)f), we observe a linear pressure decrease from the inflow to the outflow and a parabolic velocity profile everywhere except for the clamp region. The comparison to the velocity profile of the approximately circular pipe shows that the flow speed computed for the clamped artery is slower at , , (Figure (a)a). This is to be expected as the same pressure gradient has to overcome the additional obstacle of the clamp region (Figure (b)b).
5 Conclusion and Outlook
In this paper, we presented the generation and application of four-dimensional SST meshes that allow for a boundary conforming discretisation of spatial domains with time variant topology. To produce pentatope meshes of complex geometries, the elastic mesh update method was extended to four dimensions. We described the integration of 4DEMUM into the simulation workflow. Namely, we described the steps from the space-time geometry to the SST mesh of the physical domain, as well as the post-processing steps to visualise the solution on the four-dimensional mesh as series of data sets on three-dimensional meshes. The workflow was successfully applied to two test cases featuring geometries of a valve and a clamped artery. The transient three-dimensional flow solutions validate the mesh generation in the sense that proper pentatope finite element meshes are obtained.
As an outlook, the enhanced meshing capabilities open up a path to parallel-in-time computations on complex domains. We applied domain-decomposition not only to the spatial domain, but to the complete space-time domain. Another promising application area are FSCI simulations with topology changes on boundary conforming meshes. For this application, the topology changes have to be included in the SST mesh, however, spatial and temporal position of the topology changes can be determined in the course of a coupled FSI simulation and adjusted using 4DEMUM.
-  Ager C, Schott B, Vuong AT, Popp A, Wall WA. A consistent approach for fluid-structure-contact interaction based on a porous flow model for rough surface contact. International Journal for Numerical Methods in Engineering 2019; 119(13): 1345-1378.
Burman E, Claus S, Hansbo P, Larson MG, Massing A. CutFEM: discretizing geometry and partial differential equations.International Journal for Numerical Methods in Engineering 2015; 104(7): 472–501.
-  Sudhakar Y, Moitinho de Almeida J, Wall WA. An accurate, robust, and easy-to-implement method for integration over arbitrary polyhedra: Application to embedded interface methods. Journal of Computational Physics 2014; 273: 393 - 415.
-  Kadapa C, Dettmer W, Perić D. A stabilised immersed framework on hierarchical b-spline grids for fluid-flexible structure interaction with solid–solid contact. Computer Methods in Applied Mechanics and Engineering 2018; 335: 472–489.
-  Kamensky D, Hsu MC, Schillinger D, et al. An immersogeometric variational framework for fluid–structure interaction: Application to bioprosthetic heart valves. Computer Methods in Applied Mechanics and Engineering 2015; 284: 1005 - 1053.
-  Sathe S, Tezduyar TE. Modeling of fluid–structure interactions with the space–time finite elements: contact problems. Computational Mechanics 2008; 43(1): 51.
-  Bogaers AE, Kok S, Reddy BD, Franz T. An evaluation of quasi-Newton methods for application to FSI problems involving free surface flow and solid body contact. Computers & Structures 2016; 173: 71–83.
-  Rendall TC, Allen CB, Power ED. Conservative unsteady aerodynamic simulation of arbitrary boundary motion using structured and unstructured meshes in time. International Journal for Numerical Methods in Fluids 2012; 70(12): 1518–1542.
-  von Danwitz M, Karyofylli V, Hosters N, Behr M. Simplex space-time meshes in compressible flow simulations. International Journal for Numerical Methods in Fluids 2019; 91(0): 29-48.
-  Behr M. Simplex space-time meshes in finite element simulations. International Journal for Numerical Methods in Fluids 2008; 57(9): 1421–1434.
-  Karabelas E, Neumüller M. Generating admissible space-time meshes for moving domains in (d + 1) dimensions. In: Langer U, Steinbach O. , eds. Space-Time Methods Berlin, Boston: De Gruyter. 2019 (pp. 185 - 206).
-  Wang L. Discontinuous galerkin methods on moving domains with large deformations. PhD thesis. University of California, Berkeley; 2015.
-  Choudhary A, Kachanovich S, Wintraecken M. Coxeter triangulations have good quality. Mathematics in Computer Science 2020; 14(1): 141–176.
-  Neumüller M, Steinbach O. Refinement of flexible space–time finite element meshes and discontinuous Galerkin methods. Computing and Visualization in Science 2011; 14(5): 189–205.
-  Grande J. Red–green refinement of simplicial meshes in dimensions. Mathematics of Computation 2019; 88(316): 751–782.
-  Caplan PC, Haimes R, Darmofal DL, Galbraith MC. Anisotropic geometry-conforming d-simplicial meshing via isometric embeddings. Procedia Engineering 2017; 203: 141–153.
-  Caplan PC, Haimes R, Darmofal DL, Galbraith MC. Four-Dimensional Anisotropic Mesh Adaptation. Computer-Aided Design 2020; 129: 102915.
-  Caplan PCD. Four-dimensional anisotropic mesh adaptation for spacetime numerical simulations. PhD thesis. Massachusetts Institute of Technology, Cambridge; 2019.
-  Langer U, Neumüller M, Schafelner A. Space-Time Finite Element Methods for Parabolic Evolution Problems with Variable Coefficients. In: Apel T, Langer U, Meyer A, Steinbach O. , eds. Advanced Finite Element Methods with Applications: Selected Papers from the 30th Chemnitz Finite Element Symposium 2017Cham: Springer International Publishing. 2019 (pp. 247–275).
Steinbach O, Yang H. Space-time finite element methods for parabolic evolution equations: discretization, a posteriori error estimation, adaptivity and solution. In: Langer U, Steinbach O. , eds.Space-Time Methods Berlin, Boston: De Gruyter. 2019 (pp. 207 - 248).
-  Voronin K, Lee CS, Neumüller M, Sepulveda P, Vassilevski PS. Space-time discretizations using constrained first-order system least squares (CFOSLS). Journal of Computational Physics 2018; 373: 863–876.
-  Karyofylli V, Frings M, Elgeti S, Behr M. Simplex space-time meshes in two-phase flow simulations. International Journal for Numerical Methods in Fluids 2018; 86: 218-230.
-  Karyofylli V, Wendling L, Make M, Hosters N, Behr M. Simplex space-time meshes in thermally coupled two-phase flow simulations of mold filling. Computers & Fluids 2019; 192: 104261.
-  Gesenhues LG. Advanced methods for finite element simulation of rheology models for geophysical flows. Dissertation. RWTH Aachen University, Aachen; 2020.
-  Lynch DR. Unified approach to simulation on deforming elements with application to phase change problems. Journal of Computational Physics 1982; 47(3): 387–411.
-  Johnson AA, Tezduyar TE. Mesh update strategies in parallel finite element computations of flow problems with moving boundaries and interfaces. Computer Methods in Applied Mechanics and Engineering 1994; 119(1): 73–94.
-  Hübner B, Walhorn E, Dinkler D. A monolithic approach to fluid–structure interaction using space–time finite elements. Computer Methods in Applied Mechanics and Engineering 2004; 193(23-26): 2087–2104.
-  Bazilevs Y, Calo VM, Hughes TJ, Zhang Y. Isogeometric fluid-structure interaction: theory, algorithms, and computations. Computational Mechanics 2008; 43(1): 3–37.
-  Bazilevs Y, Takizawa K, Tezduyar TE. Computational fluid-structure interaction: methods and applications. John Wiley & Sons . 2013.
-  Spenke T, Hosters N, Behr M. A multi-vector interface quasi-Newton method with linear complexity for partitioned fluid–structure interaction. Computer Methods in Applied Mechanics and Engineering 2020; 361: 112810.
-  La Spina A, Förster C, Kronbichler M, Wall WA. On the role of (weak) compressibility for fluid-structure interaction solvers. International Journal for Numerical Methods in Fluids 2020; 92(2): 129–147.
-  Liu J, Yang W, Lan IS, Marsden AL. Fluid-structure interaction modeling of blood flow in the pulmonary arteries using the unified continuum and variational multiscale formulation. Mechanics Research Communications 2020; 107: 103556.
-  Zwicke F, Eusterholz S, Elgeti S. Boundary-conforming free-surface flow computations: Interface tracking for linear, higher-order and isogeometric finite elements. Computer Methods in Applied Mechanics and Engineering 2017; 326: 175–192.
-  Wendling L, Behr M, Hopf A, Kraemer F, Weber C, Turner P. CFD Simulation of Oil Jets for Piston Cooling Applications Comparing the Level Set and the Volume of Fluid Method. SAE International Journal of Advances and Current Practices in Mobility 2019; 1(2019-01-0155): 550–561.
-  Geuzaine C, Remacle JF. Gmsh: A 3-D finite element mesh generator with built-in pre-and post-processing facilities. International Journal for Numerical Methods in Engineering 2009; 79(11): 1309–1331.
-  Fernández-Fernández JA. Mesh Projector. https://github.com/JoseAntFer/MeshProjector; 2019.
-  Ayachit U. The ParaView Guide: A Parallel Visualization Application. Clifton Park, New York: Kitware Inc. . 2015.
-  Xu F, Moutsanidis G, Kamensky D, et al. Compressible flows on moving domains: Stabilized methods, weakly enforced essential boundary conditions, sliding interfaces, and application to gas-turbine modeling. Computers & Fluids 2017; 158: 201-220.
-  Pacheco DRQ, Steinbach O. A continuous finite element framework for the pressure Poisson equation allowing non-Newtonian and compressible flow behavior. International Journal for Numerical Methods in Fluids 2020. doi: 10.1002/fld.4936
-  Owens RG. A new microstructure-based constitutive model for human blood. Journal of Non-Newtonian Fluid Mechanics 2006; 140(1): 57 - 70.
-  White FM. Viscous fluid flow. McGraw-Hill New York . 2006.