We evaluate a number of different finite-element approaches for fluid–structure (contact) interaction problems against data from physical experiments. This consists of trajectories of single particles falling through a highly viscous fluid and rebounding off the bottom fluid tank wall. The resulting flow is in the transitional regime between creeping and turbulent flows. This type of configuration is particularly challenging for numerical methods due to the large change in the fluid domain and the contact between the wall and the particle. In the finite-element simulations, we consider both rigid body and linear elasticity models for the falling particles. In the first case, we compare the results obtained with the well-established Arbitrary Lagrangian–Eulerian (ALE) approach and an unfitted moving domain method together with a simple and common approach for contact avoidance. For the full fluid–structure interaction (FSI) problem with contact, we use a fully Eulerian approach in combination with a unified FSI-contact treatment using Nitsche's method. For higher computational efficiency, we use the geometrical symmetry of the experimental setup to reformulate the FSI system into two spatial dimensions. Finally, we show full three-dimensional ALE computations to study the effects of small perturbations in the initial state of the particle to investigate deviations from a perfectly vertical fall observed in the experiment. The methods are implemented in open-source finite element libraries, and the results are made freely available to aid reproducibility.

Flows containing particles, i.e., particulate flows or particles settling in a fluid, have many industrial and biological applications. Examples range from the transport of blood cells in blood flows1 to the simulation of fluidized bed reactors.2 

We shall consider single elastic spherical particles falling freely in a viscous fluid and rebounding off the bottom wall of the fluid domain at Reynolds numbers in the transitional regime between creeping and turbulent flows. The multiphase and fluid–structure interaction (FSI) problem with solid contact posed by the settling in the fluid and rebounding off a wall is challenging from both an analytical and a numerical perspective.

From the theoretical point of view, the correct model for the transition to contact with the bottom wall is not yet fully understood. In the case where a rigid solid is assumed, most flow models lead to results contradicting real world observations. For example, if a creeping flow is assumed such that the linear Stokes equations are applicable, then contact can only occur under singular forces, cf. Ref. 3. When the non-linear incompressible Navier–Stokes together with no-slip boundary conditions are taken for the fluid model, then contact cannot occur and it is impossible to release contact.4 This can however be overcome; if the boundary condition is modified to a free-slip condition,5 the rough nature of the surface is taken into account6 or the fluid is taken to be compressible.7 If the solid model is changed to take the elasticity of the body into account, then it is currently assumed that even with perfectly smooth boundaries and incompressibility, rebounding without contact can occur due to the storage of energy in the elastic solid.8–10 This has been refined recently (in the Stokes setting),11 where it has been shown that the internal storage of energy is not sufficient but that additionally a change in the “flatness” is necessary to achieve physically meaningful rebound without topological contact.

For numerical methods, the challenge lies in the discretization of the resultant FSI system.12 It consists of a free boundary value problem with a moving interface. The most well-established method for this is the Arbitrary Lagrangian–Eulerian (ALE) approach.13 This approach leads to very efficient and accurate computations in situations where the method is usable. However, its usage is limited as it breaks down when deformations with respect to the reference configuration become too large and when contact occurs.14 The lattice-Boltzmann method is an alternative approach that is also able to describe the behavior of settling particles.15 To deal with large deformations, overlapping mesh techniques have been developed.16 Here, the background fluid domain and the region around the structure are meshed separately so that the fluid–solid interface is resolved. The two meshes are then coupled using unfitted approaches. This then allows a hybrid approach, where the solid and the near fluid are treated using the ALE framework, while the remaining fluid is treated in Eulerian coordinates.17 To overcome both large deformations and contact, fully Eulerian approaches have lately become the focus of research. In the case of rigid bodies, a number of different approaches have been considered, for example, based on fictitious domain methods using Lagrange multipliers,18 XFEM type approaches,19,20 and most recently, CutFEM approaches21 using Nitsche's method.22,23 Here, a major issue remains of achieving a realistic rebound effect since an artificial contact/lubrication force is added to the equation governing the motion of the solid to prevent the overlap of the solid regions.18 Nevertheless, topological changes appear to be unproblematic for the CutFEM type approaches.24 

Considering full fluid–structure interactions, immersed approaches have become popular in recent years.25–28 Here, the fluid and the solid are treated in their natural Eulerian and Lagrangian coordinate systems, respectively, and the sub-domains are meshed separately. The two meshes are then coupled by means of Nitsche's method25,27 or using Lagrange multipliers.26,28,29 Another possibility to handle large deformations and contact is by using fully Eulerian approaches, where both the solid and fluid equations are formulated in the Eulerian coordinate framework, which simplifies the coupling within monolithic algorithms.30–34 All these approaches are however relatively new and require further development with respect to accuracy and robustness.

The aforementioned methods have been applied to different test cases for numerical validation, and a priori error estimates are also available in most cases. The established benchmarks for fluid-structure interaction problems such as those in Ref. 35 completely avoid contact since the methods that handle contact remain relatively new. For rigid-body motion, most numerical studies are interpreted qualitatively or compared to artificial, analytically derived solutions. Especially in the cases where artificial forces are introduced in order to avoid contact of rigid solids, real validation is near impossible as this introduces model parameters for which there is no a priori knowledge on a good choice. However, a number of FSI methods have recently become available that are able to resolve contact.36–40 This then raises the question of how well the different modeling and discretization approaches depict the behavior of contact and rebound observed in physical experiments.

In this work, we take recently published data from experiments where different solid spherical particles were allowed to settle in a viscous fluid.41,42 We then use a rigid-body ALE, a rigid-body Eulerian CutFEM, and a fully Eulerian FSI approach to simulate the scenarios presented by the physical experiments. This aims to show the validity and the applicability of these different approaches to the different aspects/problems posed by this process. Furthermore, we will illustrate how spatially reduced models are able to capture the behavior in comparison to full three-dimensional computations. To the best of our knowledge, there is currently no comparable benchmark that considers such a multiphase flow/fluid-structure interaction problem with contact, which is validated against experimental data.

The remainder of this paper is structured as follows: In Sec. II, we describe the considered problem, that is, a description of the physical experiment, the mathematical models used to describe the experiment, the specific setups we will simulate, and the quantities used to compare the numerical simulations with the experimental data. Section III then briefly covers the reduced formulation we apply to increase the computational efficiency in our numerical methods. The numerical computations are then presented in Sec. IV; we present the details of the different numerical approaches in Subsection IV A, and the results are then presented in Subsection IV B. We discuss the conclusions from these results in Sec. V and consider the aspects that remain open. Furthermore, we define and compute two simplified setups in  Appendix B, designed to help others reproduce the presented computational results.

We describe the experimental setup used to gather the data and the mathematical model that we will use to reproduce the behavior observed in the experiments, and we define relevant quantities used to compare the results quantitatively.

The experiments in Ref. 42 capture the settling and impact process of spherical particles with different sizes and densities in a cylindrical tank. The latter contains a liquid mixture consisting of glycerine and water at equal volume fractions. At the bottom of the cylindrical tank, which is made of acrylic glass for optical access, a massive steel anvil serves as the impact object. Moreover, the cylinder is surrounded by a rectangular container, filled with refractive index matching liquid, to compensate for optical distortion coming from the curved cylinder walls. The filling level allows the observation of the particle settling along a vertical distance of 140 mm–160 mm, depending on the particle size. Initially, each particle is held in place and submerged in the liquid by using a vacuum tweezer. The particle is then released by switching off the vacuum pump. The particle is tracked during the settling process and the impact on the steel anvil, including the rebound, using a high-speed CMOS-camera. This acquires shadow images at a frame rate of 1000 fps and a scale factor of 8.89 pixel/mm. An image processing algorithm coded in MATLAB yields the in-plane particle coordinates as a function of time and allows us to extract the instantaneous particle settling velocity.

The resulting data are available via Mendeley Data.41 These data are the basis for our comparison and validation of the numerical code.

We consider a bounded domain Ω R d , with d ∈ {2, 3}, over a finite, non-empty, time interval [0, Tend]. This is divided into a d-dimensional fluid region F , a d-dimensional solid S , and a d − 1-dimensional interface I dividing the solid and fluid regions. For these, we have Ω = F ̇ I ̇ S .

1. Fluid model

In the time dependent fluid region F ( t ) , we consider the incompressible Navier–Stokes equations. Find a velocity u and a pressure p such that
ρ f t u + ( u ) u div σ ( u , p ) = 0 ,
(1a)
div ( u ) = 0
(1b)
with the non-symmetric stress tensor
σ ( u , p ) = μ f u p Id ,
where μf = ρfνf is the fluid's dynamic viscosity, νf is the kinematic viscosity, ρf is the fluid's density, and Id is the d-dimensional identity operator. Appropriate boundary conditions to complete this system will be discussed later. Note that we use bold face letters to denote vector/matrix valued quantities, while regular faced letters denote scalar objects.

2. Elastic solid and fluid–structure interaction

We consider a linear elastic solid model in S for the solid displacement d and the solid velocity d ̇ , given by
ρ s t d ̇ div σ s d = ( ρ s ρ f ) g , d ̇ = t d ,
(2)
with the acceleration due to gravity g = −9.807 m s−2 and the Cauchy stress tensor σs defined by
σ s ( d ) = 2 μ s E ( d ) + λ s tr ( E ( d ) ) I , E ( d ) = 1 2 d + d T ,
where μs and λs denote the Lamé parameters. Note that we have subtracted the fluid gravitational force ρfg on the right-hand side of the solid equation to be consistent with the equations for a rigid solid presented in the following paragraph. Alternatively, this could be added to the right-hand side of the fluid equations.
Solid and fluid are coupled by means of no-slip coupling conditions on I ,
u = d ̇ , σ ( u , p ) n = σ s ( d ) n .
(3)
The current position of the interface I ( t ) is determined by the displacement variable d.

3. Rigid solid

As the solid materials that we consider are relatively hard, the consideration of a rigid solid yields a good approximation of the FSI dynamics, at least up to the moment when the solid comes close to the lower wall. The movement of the solid is governed by Newton's second law of motion. Let c S ( t ) be the center of mass of the solid S . Since we will consider spherical particles, this is then governed by
d 2 d t 2 c S ( t ) m S = f s ,
(4)
where m S is the mass of the solid and fs are the forces acting on the solid in the horizontal and vertical directions. For simplicity, we assume that the horizontal forces are negligible. The vertical forces are then the gravitational pull, buoyancy, and viscous drag,
f s = 0 0 m S g vol ( S ) ρ f g + F 3 ,
(5)
where vol ( S ) is the volume of the solid and F3 is the viscous drag force in the vertical direction. This is the third component of
F = I σ ( u , p ) n d s .
(6)
Note that we added the effects of buoyancy in (5). This would be naturally included in F if we added body force ρfgxd to the right-hand side of the fluid equation (1a). However, since this would only affect the pressure, it is sufficient to consider the homogeneous equation (1a) and add the effect of buoyancy to (5). Furthermore, this approach is more accurate in our case since only pressure-robust methods are able to reflect this effect of gradient contributions in the forcing term on the pressure exactly on the numerical level.43 
As a result of us neglecting the horizontal movement of the solid, (4) becomes a scalar ordinary differential equation (ODE). We can also simplify the terms in (5) so that, in total, we come to the equation
d d t v S , 3 ( t ) = ρ s ρ f ρ s g + F 3 vol ( S ) ρ s ,
(7)
where v S ( t ) = d d t c S ( t ) is the solid's velocity.
The solid's motion couples back to the fluid equations through the boundary condition at the interface I by requiring the continuity of the velocity, i.e.,
u I = v S = d d t c S ( t ) .
(8)
We note that this model neglects rotational effects. As will be shown below, this does not have a major impact on the quality of the resulting approximations.

Since we consider balls of different sizes, we shall keep the domain description general. The background domain is a cylinder Ω = { x = ( x , y , z ) T R 3 | x 2 + y 2 < R 2 , 0 < z < H } for a given radius R and a given height H. At t = 0, the solid domain is described by S ( 0 ) = { x = ( x , y , z ) T R 3 | x 2 + y 2 + ( z ( h 0 + r S ) ) 2 < r S 2 } for a given ball radius r S and an initial height of the bottom of the ball h0. Accordingly, the volume of the solid is given by vol ( S ) = 4 π r 3 / 3 . An illustration of this can be seen on the left of Fig. 1.

FIG. 1.

The initial spatial configuration: (a) three-dimensional domain and (b) rotationally reduced domain.

FIG. 1.

The initial spatial configuration: (a) three-dimensional domain and (b) rotationally reduced domain.

Close modal

1. Boundary conditions

We denote the top boundary (z = H) of the cylinder as Γtop, the bottom boundary (z = 0) as Γbottom, and the side of the cylinder (x2 + y2 = R2) as Γwall.

On the interface I between the solid and the fluid, the Dirichlet boundary condition is given by the continuity of the velocity; see (3) and (8). On the wall and bottom boundaries Γbottom ∪ Γwall, we shall impose homogeneous Dirichlet boundary conditions u = 0. In order to approximate the free surface at the top of the water tank Ω, we impose a free-slip boundary condition u3 = 0 at Γtop.

1. Spatial parameters

We consider two different particles: one with diameter d S = 6 mm consisting of polytetrafluoroethylene/teflon (PFTE) and one with diameter d S = 22 mm consisting of rubber. We shall refer to these two cases as PTFE6 and Rubber22, respectively.

Both particles are considered inside the same cylindrical fluid tank with radius R = 0.055 m and a height of H = 0.2 m. The spatial setup for both cases is summarized in Table I.

TABLE I.

Geometrical parameters of the test cases.

Geometry Boundary conditions
Experiment R (m) H (m) r S (m) h0 (m) Γwall ∪ Γbottom Γtop I
PTFE6  0.055  0.2  0.003  0.161 661 6  u = 0  ud = 0  u = v S  
Rubber22  0.011  0.146 120 3 
Geometry Boundary conditions
Experiment R (m) H (m) r S (m) h0 (m) Γwall ∪ Γbottom Γtop I
PTFE6  0.055  0.2  0.003  0.161 661 6  u = 0  ud = 0  u = v S  
Rubber22  0.011  0.146 120 3 

2. Material parameters

The numerical computation of the above-mentioned scenarios requires the following material parameters. For the fluid, these are the density and viscosity. For the solid balls, we require the density and the Lamé parameters. Both fluid and solid parameters are summarized in Table II.

TABLE II.

Summary of the benchmark setup in standard units.

Material parameters
Experiment g (m s−2) μf (kg m−1 s−1) ρf (kg m−3) ρs (kg m−3) λs (kg m−1 s−2) μs (kg m−1 s−2)
PTFE6  −9.807  0.008  1141  2122  2.638 70 × 109  2.294 52 × 108 
Rubber22  1361  3.332 89 × 109  6.667 11 × 105 
Material parameters
Experiment g (m s−2) μf (kg m−1 s−1) ρf (kg m−3) ρs (kg m−3) λs (kg m−1 s−2) μs (kg m−1 s−2)
PTFE6  −9.807  0.008  1141  2122  2.638 70 × 109  2.294 52 × 108 
Rubber22  1361  3.332 89 × 109  6.667 11 × 105 

Remark

(source of material parameters). The fluid parameters and the solid densities are given in the original experimental paper.42 We derived the Lamé parameters from the Young's modulus and Poisson ratio of the used materials. In the case of the PTFE6 ball, these are about 670 MPa = 670 000 000 kg m−1 s−2 (https://www.kugelpompel.at/upload/2312783_Datenblatt%20Kunststoffkugel%20PTFE%20V1.01.pdf) and νs = 0.46 (http://www.matweb.com/search/datasheet_print.aspx?matguid=4e0b2e88eeba4aaeb18e8820f1444cdb), respectively.

The chemical analysis of the Rubber22 material suggests this to be hydrogenated nitrile rubber (HNBR). The Young's modulus of this is ∼1.7 MPa to 20.7 MPa = 1 700 000 kg m−1 s2–20 700 000 kg m−1 s2 (https://eriks.de/content/dam/de/pdf/downloads/dichtungen/o-ringe/ERIKS_Technisches-Handbuch-O-Ringe_de.pdf, p. 20), and the Poisson ratio is νs = 0.4999.44 

We will use the following quantities to compare our numerical results with each other and with the experimental data:

  • t* let t 0 = t c S = h 0 be the time at which the center of mass is at h0, i.e., when the ball has traveled r S vertically. We define t* as the time after release when dist ( I , Γ bottom ) = d S relative to t0, i.e., for PTFE6 t * = t c S = ( 0 , 0 , 0.009 ) t 0 and for Rubber22 t * = t c S = ( 0 , 0 , 0.033 ) t 0 .

  • v* the velocity of the ball in the z-direction at t = t* + t0.

  • f* the vertical component of the force F acting on the ball at time t* + t0.

  • tcont the time of the first solid contact relative to t0.

  • tjump the time between contact and the second change in direction is realized, i.e., the amount of time the balls travels upward after the first contact.

  • djump the maximum of dist ( I , Γ bottom ) after contact, i.e., the size of the bounce.

An illustration of how these quantities are defined can be seen in Fig. 2.

FIG. 2.

Illustration for the definitions of the quantities of interest.

FIG. 2.

Illustration for the definitions of the quantities of interest.

Close modal

Remark

(choice of reference values). In the experiments, we observed that the balls do not immediately start to fall after they are released. The settling process starts with some distance to the liquid surface but in the closest vicinity of the vacuum cup. In general, particles experience an increased drag force when they are moving toward or away from a solid wall and free fluid-surfaces. Accordingly, the early stage of the settling process here is dominated by an increased drag force coming from the vacuum cup and the liquid surface.

Due to the rather slow motion of the sphere at the beginning, the moment of release cannot be defined well. To be able to compare the numerical results with the experimental data, we therefore defined the above-mentioned quantities relative to the point in time at which the ball has traveled the distance of the ball's radius. Furthermore, since we have no measurement of the drag force in in the experiment, f* will only be used to compare the computations directly.

To establish v*, tcont, dcont, and djump from the experimental data, we interpolate the height data using a spline of order 3. This spline is then evaluated to establish the time at which c S = ( 0 , 0 , 3 r S / 2 ) , the time of contact, as well as the time and height of the jump. The velocity is then taken as the first derivative of the spline representing the height. The resulting reference values for the quantities of interest are shown together with our numerical results in Tables III and IV, respectively.

TABLE III.

Results for the PTFE6 setup. The experimental values have been reproduced with permission from T. Hagemeier, Particle settling-transitional regime, version 1, Mendeley Data, September 2020. Copyright 2020 Author(s), licensed under a Creative Commons Attribution 4.0 License.

Discretization Results
Method [hmin, hmax] Δt dof nze t* v* f* tcont tjump djump
ALE  [0.000 08, 0.004]  1/200  7.54  0.34  0.542 134  −0.312 207 7  1.113 93 · 10−3  …  …  … 
  [0.000 04, 0.002]  1/800  23.88  1.10  0.539 231  −0.313 924 5  1.119 82 · 10−3  …  …  … 
  [0.000 02, 0.001]  1/3 200  83.07  3.90  0.539 026  −0.313 993 1  1.120 08 · 10−3  …  …  … 
      Extrapolate  0.539 010  −0.313 996 0  1.120 21 · 10−3  …  …  … 
      order (in h 3.8  4.6  4.0  …  …  … 
CutFEM  [0.001 14, 0.008]  1/2 000  14.50  0.43  0.579 181  −0.288 446 2  1.118 97 · 10−3  0.600 328  0.017 827  1.339 74 · 10−3 
  [0.000 57, 0.004]  1/2 000  45.80  1.35  0.536 698  −0.316 408 6  1.127 23 · 10−3  0.556 037  0.030 627  2.706 81 · 10−3 
  [0.000 29, 0.002]  1/2 000  152.06  4.50  0.530 537  −0.321 376 2  1.121 22 · 10−3  0.549 543  0.030 749  2.692 50 · 10−3 
FSI  [0.000 32, 0.00562]  [1/128 000, 1/1 000]  59.5  1.24  0.539 918  −0.314 713 2  …  0.558 503  0.027 515  1.817 16 · 10−3 
  [0.000 16, 0.002 81]  [1/128 000, 1/1 000]  236.6  5.40  0.535 625  −0.315 960 0  …  0.554 358  0.028 015  2.061 63 · 10−3 
Experiment          0.516 403  −0.330 987  …  0.534 503  0.027 92  2.211 70 · 10−3 
Discretization Results
Method [hmin, hmax] Δt dof nze t* v* f* tcont tjump djump
ALE  [0.000 08, 0.004]  1/200  7.54  0.34  0.542 134  −0.312 207 7  1.113 93 · 10−3  …  …  … 
  [0.000 04, 0.002]  1/800  23.88  1.10  0.539 231  −0.313 924 5  1.119 82 · 10−3  …  …  … 
  [0.000 02, 0.001]  1/3 200  83.07  3.90  0.539 026  −0.313 993 1  1.120 08 · 10−3  …  …  … 
      Extrapolate  0.539 010  −0.313 996 0  1.120 21 · 10−3  …  …  … 
      order (in h 3.8  4.6  4.0  …  …  … 
CutFEM  [0.001 14, 0.008]  1/2 000  14.50  0.43  0.579 181  −0.288 446 2  1.118 97 · 10−3  0.600 328  0.017 827  1.339 74 · 10−3 
  [0.000 57, 0.004]  1/2 000  45.80  1.35  0.536 698  −0.316 408 6  1.127 23 · 10−3  0.556 037  0.030 627  2.706 81 · 10−3 
  [0.000 29, 0.002]  1/2 000  152.06  4.50  0.530 537  −0.321 376 2  1.121 22 · 10−3  0.549 543  0.030 749  2.692 50 · 10−3 
FSI  [0.000 32, 0.00562]  [1/128 000, 1/1 000]  59.5  1.24  0.539 918  −0.314 713 2  …  0.558 503  0.027 515  1.817 16 · 10−3 
  [0.000 16, 0.002 81]  [1/128 000, 1/1 000]  236.6  5.40  0.535 625  −0.315 960 0  …  0.554 358  0.028 015  2.061 63 · 10−3 
Experiment          0.516 403  −0.330 987  …  0.534 503  0.027 92  2.211 70 · 10−3 
TABLE IV.

Results for the Rubber22 setup. The experimental values have been reproduced with permission from T. Hagemeier, Particle settling-transitional regime, version 1, Mendeley Data, September 2020. Copyright 2020 Author(s), licensed under a Creative Commons Attribution 4.0 License.

Discretization Results
Method [hmin, hmax] Δt dof nze t* v* f* tcont tjump djump
ALE  [0.000 4, 0.004]  1 / 200   6.9  0.31  0.455 904 5  −0.303 002  1.132 62 · 10−2  …  …  … 
  [0.000 2, 0.002]  1 / 800   21.5  0.99  0.455 355 1  −0.303 552  1.131 07 · 10−2  …  …  … 
  [0.000 1, 0.001]  1 / 3200   73.8  3.45  0.455 333 4  −0.303 616  1.131 17 · 10−2  …  …  … 
      Extrapolate  0.455 332 5  −0.303 625  …  …  …  … 
      order (in h 4.6  3.1  …  …  …  … 
ALE 3D  [0.000 8, 0.032]  1 / 200   18.43  4.05  0.453 695  −0.304 252  1.142 52 · 10−2  …  …  … 
  [0.000 4, 0.016]  1 / 800   68.52  15.37  0.456 145  −0.302 170  1.144 93 · 10−2  …  …  … 
  [0.000 2, 0.008]  1 / 3200   304.4  69.33  0.455 929  −0.302 667  1.135 57 · 10−2  …  …  … 
CutFEM  [0.002 00, 0.008]  1 / 2000   11.46  0.34  0.453 455  −0.309 542 3  1.090 94 · 10−2  0.524 979  0.098 086  7.308 15 · 10−3 
  [0.001 00, 0.004]  1 / 2000   38.98  1.16  0.454 081  −0.308 841 2  1.110 82 · 10−2  0.525 502  0.110 672  1.147 72 · 10−2 
  [0.000 50, 0.002]  1 / 2000   140.96  4.20  0.453 789  −0.306 220 2  1.129 22 · 10−2  0.526 010  0.121 606  1.368 74 · 10−2 
FSI (Es = 5 · 106 [0.001 0, 0.004]  [ 1 / 2000 , 1 / 500 ]   51.4  1.15  0.446 020  −0.320 688 9  …  0.515 197  0.079 14  2.682 98 · 10−3 
  [0.000 5, 0.002]  [ 1 / 2000 , 1 / 500 ]   204.6  4.72  0.449 821  −0.311 385 1  …  0.521 487  0.083 323  3.719 86 · 10−3 
FSI (Es = 2 · 106 [0.001 0, 0.004]  [ 1 / 2000 , 1 / 500 ]   51.4  1.15  0.446 012  −0.320 715 1  …  0.516 197  0.092 5  5.292 47 · 10−3 
  [0.000 5, 0.002]  [ 1 / 2000 , 1 / 500 ]   204.6  4.72  0.449 827  −0.311 375 9  …  0.522 087  0.090 5  5.501 46 · 10−3 
Experiment          0.469 137  −0.309 301  …  0.544 021  0.089 492  4.414 85 · 10−3 
Discretization Results
Method [hmin, hmax] Δt dof nze t* v* f* tcont tjump djump
ALE  [0.000 4, 0.004]  1 / 200   6.9  0.31  0.455 904 5  −0.303 002  1.132 62 · 10−2  …  …  … 
  [0.000 2, 0.002]  1 / 800   21.5  0.99  0.455 355 1  −0.303 552  1.131 07 · 10−2  …  …  … 
  [0.000 1, 0.001]  1 / 3200   73.8  3.45  0.455 333 4  −0.303 616  1.131 17 · 10−2  …  …  … 
      Extrapolate  0.455 332 5  −0.303 625  …  …  …  … 
      order (in h 4.6  3.1  …  …  …  … 
ALE 3D  [0.000 8, 0.032]  1 / 200   18.43  4.05  0.453 695  −0.304 252  1.142 52 · 10−2  …  …  … 
  [0.000 4, 0.016]  1 / 800   68.52  15.37  0.456 145  −0.302 170  1.144 93 · 10−2  …  …  … 
  [0.000 2, 0.008]  1 / 3200   304.4  69.33  0.455 929  −0.302 667  1.135 57 · 10−2  …  …  … 
CutFEM  [0.002 00, 0.008]  1 / 2000   11.46  0.34  0.453 455  −0.309 542 3  1.090 94 · 10−2  0.524 979  0.098 086  7.308 15 · 10−3 
  [0.001 00, 0.004]  1 / 2000   38.98  1.16  0.454 081  −0.308 841 2  1.110 82 · 10−2  0.525 502  0.110 672  1.147 72 · 10−2 
  [0.000 50, 0.002]  1 / 2000   140.96  4.20  0.453 789  −0.306 220 2  1.129 22 · 10−2  0.526 010  0.121 606  1.368 74 · 10−2 
FSI (Es = 5 · 106 [0.001 0, 0.004]  [ 1 / 2000 , 1 / 500 ]   51.4  1.15  0.446 020  −0.320 688 9  …  0.515 197  0.079 14  2.682 98 · 10−3 
  [0.000 5, 0.002]  [ 1 / 2000 , 1 / 500 ]   204.6  4.72  0.449 821  −0.311 385 1  …  0.521 487  0.083 323  3.719 86 · 10−3 
FSI (Es = 2 · 106 [0.001 0, 0.004]  [ 1 / 2000 , 1 / 500 ]   51.4  1.15  0.446 012  −0.320 715 1  …  0.516 197  0.092 5  5.292 47 · 10−3 
  [0.000 5, 0.002]  [ 1 / 2000 , 1 / 500 ]   204.6  4.72  0.449 827  −0.311 375 9  …  0.522 087  0.090 5  5.501 46 · 10−3 
Experiment          0.469 137  −0.309 301  …  0.544 021  0.089 492  4.414 85 · 10−3 

The experimental study42 was conducted so that the horizontal displacement of the particles was minimal. The dataset that we will compare against shows a maximal horizontal displacement of less than 2 mm and 0.75 mm in the PTFE6 and Rubber22 cases, respectively. This compares with the mean over time of the maximal deviation in the center location between experiment repetitions of 0.192 mm and 0.135 mm for the PTFE6 and Rubber22 cases, respectively. However, the experiment is only able to give the projection of the horizontal displacement onto the xz plane. Thus, it is not possible to detect the true horizontal motion. Note that we have ignored the horizontal motion in the computation of the reference values. However, since overall the horizontal deflection is small, we consider this to be reasonable for the present purpose.

The setup described in Subsection II C is symmetric with respect to rotation if viewed in cylindrical coordinates. The experimental data presented in Refs. 41 and 42 show a rotational component in the motion of the solid, and they also show a small deflection of the center of mass c S from the z axis. However, these effects are small, and since the material parameters are such that the resulting flow is in the low to intermediate Reynolds-number regime,42 we assume that the solution is described sufficiently well by a rotationally symmetric flow in cylindrical coordinates. We use this in order to obtain a two-dimensional reduced formulation that is computationally cheaper compared to full three-dimensional computations. In Subsection IV B 3, we will also present a fully resolved three-dimensional simulation to have a closer look at these rotational effects. To distinguish between the full three-dimensional domains and the reduced two-dimensional domain, we shall denote objects stemming from the three-dimensional setup with a superscript “3d” if there is a potential for ambiguity. For the sake of readability, objects based on the reduced two-dimensional setup will not have a special notation.

With the spaces V 3 d = { v H 1 ( F 3 d ) | v I 3 d = u S 3 d , v Γ wall 3 d Γ bottom 3 d = 0 and v z Γ top 3 d = 0 } [ H 1 ( F 3 d ) ] 3 and Q 3 d = L 0 2 ( F 3 d ) , the weak formulation of the Navier–Stokes equation (1) is as follows:
Find  ( u , p ) V 3 d × Q 3 d  such that for all  ( v , q ) V 3 d × Q 3 d ,  it holds  A f 3 d ( u , p ; v , q ) m 3 d ( t u , v ) + a 3 d ( u , v ) + c 3 d ( u , u , v ) + b 3 d ( v , p ) + b 3 d ( u , q ) = 0
(9)
with the multilinear forms
m 3 d ( u , v ) = ρ f F 3 d u v d x , a 3 d ( u , v ) = μ f F 3 d u : v d x , c 3 d ( u , v , w ) = ρ f F 3 d ( u ) v w d x , b 3 d ( q , v ) = F 3 d q v d x .
(10)

In order to reduce this three-dimensional flow problem into a two-dimensional flow problem, we rewrite the problem into cylindrical coordinates r, ϕ, z. The rotational symmetry of a solution u means that ϕu = 0. We use this to rotate the domain into the r+z-plane and transform the weak formulation (9) accordingly. Now, let = ( r z ) and Ω = Ω 3 d ( R + × R ) be the reduced two-dimensional computational domain. A sketch of the three-dimensional domain transformed into two dimensions can be seen on the right of Fig. 1.

At the symmetry boundary (r = 0), the boundary conditions
u r = 0 , r u z = 0
(11)
are valid. The reduced spaces are then V = { v H 1 ( F ) | v I = 0 v S 3 , v Γ wall Γ bottom = 0 ,  v z Γ top = 0 , and v r r = 0 = 0 } [ H 1 ( F ) ] 2 and Q = L 0 2 ( F ) . The reduced Navier–Stokes problem then reads as follows:
Find  ( u , p ) V × Q  such that for all  ( v , q ) V × Q ,  it holds  A f ( u , p ; v , q ) m ( t u , v ) + a ( u , v ) + c ( u , u , v ) + b ( v , p ) + b ( u , q ) = 0
(12)
with the transformed multilinear forms
m ( u , v ) = 2 π ρ f F r u v d x , a ( u , v ) = 2 π μ f F u : v + 1 r u r v r d x , c ( u , v , w ) = 2 π ρ f F r ( u ) v w d x , b ( q , v ) = F q ( v r + r v ) d x ,
cf. Ref. 45. Since the motion of the solid is partially driven by the fluid forces acting on it, we also need to transform these. Let (u3d, p3d) ∈ V3d × Q3d be a rotationally symmetric solution of (9) and (u, p) ∈ V × Q be the solution of (12). We can then transform the weak boundary integral formulations as
I 3 d σ 3 d ( u 3 d , p 3 d ) n 3 d v 3 d d S 3 d = 2 π I r ( μ f u Id p ) n v d S .
(13)
The forces F can then be computed by inserting the appropriate non-conforming test-functions v into this functional. Using this, we can compute the motion of the solid as before using (7).

We give the details of the different numerical approaches applied to attempt to reproduce the observed data and present the results attained with these methods.

The full result datasets of all methods and the source code to reproduce the rigid-body computations according to the descriptions in Subsections IV A 1 and IV A 2 can be found in the zenodo repository.63 This also includes the preparatory examples presented in  Appendix B.

We provide the details on the formulations of the different discretization approaches used.

1. Fluid–rigid body interaction in Arbitrary Lagrangian–Eulerian coordinates

To compute the spatially reduced coupled problem, we formulate the Navier–Stokes equations in Arbitrary Lagrangian–Eulerian (ALE) coordinates46 by introducing a reference map,
T A L E ( c S ) : F F ( c S ) ,
and by transforming the Navier–Stokes equation onto the reference domain F , which is fixed for all times. This setting allows for a direct finite-element triangulation Ωh of F that resolves the interface between fluid and solid well. ALE approaches are highly efficient and accurate and well established for fluid–solid interaction problems. ALE approaches will however fail if the deformation becomes too large47 or if even contact of the solid with an outer boundary happens,14 such as it is the case in our context. We refer to the literature on ALE12,46 and to Appendix A 1 for details on our implementation.

2. Fluid–rigid body interaction in Eulerian coordinates

As in Subsection IV A 1, we consider the problem as a moving domain problem for the fluid, assume the solid to be a rigid body, and decouple the fluid and solid equations. For the resulting moving domain problem, we use an unfitted Eulerian finite-element method from Ref. 23 using BDF2 time stepping, together with Taylor–Hood elements in space, which are inf–sup stable in the CutFEM setting48 with ghost-penalty stabilization.49 In order to obtain the full convergence order of the Taylor–Hood finite-element pair, we use an isoparametric mapping introduced in Ref. 50 for stationary domains to realize the higher-order geometry approximation.

a. Transformed nitsche terms
In the CutFEM method used here, Dirichlet boundary conditions on unfitted boundaries are enforced using Nitsche's method.51 For a consistent and stable method, we also transform these terms into the rotationally symmetric formulation. Using the standard derivation for Nitsche's method, we find that the consistency and penalty terms for the reduced formulation are
n c ( u , v ) = 2 π μ f Γ r ( u ) n v d S
and
n s ( u , v ) = 2 π μ f σ k 2 h Γ r u v d S ,
respectively, with a penalty parameter σ > 0, the velocity space's polynomial order k, and the local mesh diameter h. Similarly, we find for the pressure-coupling operator that the Nitsche term is
n p ( v , q ) = Γ r q v n d S .
The necessary Nitsche term for a symmetric and consistent formulation to enforce the Dirichlet conditions u = g is therefore
n ( u , p ; v , q ) n c ( u , v ) + n c ( v g , u ) + n s ( u g , v ) + n p ( v , p ) + n p ( u g , q ) .
b. Transformed ghost-penalty operators

The role of the ghost-penalty operator in the unfitted Eulerian time-stepping scheme used here is twofold. As in other CutFEM discretizations, it stabilizes arbitrary element cuts such that the method is stable and the resulting matrices are well conditioned.49 Second, appropriately scaled ghosty-penalties provide the necessary implicit extension for the method-of-lines approach to the discretization of the time-derivative.22–24 

We use the direct-version of the ghost-penalty operator.52 To define this operator, let f h be a set of facets between neighboring elements on which the ghost-penalty operator is to act on, cf. Refs. 23 and 24 for further details. For a facet F f h such that F = T ¯ 1 T ¯ 2 , we define the facet patch ωF = T1T2. The velocity ghost-penalty operator for the rotationally symmetric formulation is then
i h ( u , v ) = γ u F f h ω F r h 2 ( u 1 u 2 ) ( v 1 v 2 ) + 1 r ( u r , 1 u r , 2 ) ( v r , 1 v r , 2 ) d x ,
where v i = E P v T i with the canonical extension of polynomials E P : P T P R d . The pressure ghost-penalty operator is
j h ( p , q ) = γ p F f h ω F r ( p 1 p 2 ) ( q 1 q 2 ) d x .
With these transformed ghost-penalty operators, it is easy to show the standard ghost-penalty results in the appropriately transformed norms.
The variational formulation of the CutFEM discretization then reads as follows:
Find  ( u , p ) V h × Q h  such that for all  ( v , q ) V h × Q h ,  it holds A f ( u , p ; v , q ) + n ( u , p ; v , q ) + μ f + 1 μ f i h ( u , v ) + 1 μ f j h ( p , q ) = 0 .
c. Contact algorithm
We consider a very basic contact avoidance scheme, used widely in the literature.18,32,53,54 The idea is to introduce an artificial (lubrication) force acting on the rigid body in the vicinity of the contact wall, which increases when the ball comes closer to the wall and acts in the direction away from the wall. This force is then added to the forces governing the motion of the rigid solid such that contact does not occur. We define this force as
f c ( S ) = 0 if dist ( I , Γ bottom ) d i s t 0 γ c d i s t 0 dist ( I , Γ bottom ) dist ( I , Γ bottom ) if dist ( I , Γ bottom ) < d i s t 0 ,
where dist0 and γc are parameters to be chosen and dist ( I , Γ bottom ) is the minimal distance between I and Γbottom. We then add this to the right-hand side of the ODE (4) and carry this thought so that the right-hand side of (7) becomes ρ S ρ F ρ S g + F y + f c vol ( S ) ρ S . The solution to this ODE then governs the position on the level set describing the solid and the Dirichlet condition enforced on the interface.
d. Implementation

This discretization is implemented using the finite-element library Netgen/NGSolve (see Refs. 55 and 56 and ngsolve.org) together with the add-on package ngsxfem57 for unfitted finite-element functionality.

The background mesh is constructed by defining a local mesh parameter hinner on the left of the reduced domain where r < 2 d S / 3 and then creating a shape regular mesh with h = hmax in the remainder of the domain. This is to obtain more accurate boundary integrals, i.e., when computing the forces acting on the ball. In the Rubber22 setting, we choose hinner = hmax/4 and in the PTFE6 case as hinner = hmax/7. This is to obtain background meshes with a similar number of elements in each settings. On the active part of the mesh, we consider standard P 2 / P 1 Taylor–Hood elements.

The Nitsche parameter is taken as σ = 100, the extension ghost-penalty parameter is γu,e = 0.1, the cut-stability ghost-penalty parameter γu,s = γp,s = 0.01, and the extension strip width parameter is cδ = 4. See Ref. 23 for details on these parameters.

The contact parameters are tuned with respect to the PTFE6 jump height since the model cannot be expected to resolve the elastic nature of rubber. For PTFE6, we take the contact model parameters dist0 = 2 · 10−5 and γc = 0.38. Since the mass of the Rubber22 ball is ∼31.6 times larger than the PTFE6 ball, we take contact model parameters that are appropriately larger such that the resulting acceleration acting on the balls is comparable. For the Rubber22 computations, we take dist0 = 2 · 10−5 and γc = 12.

Each time step is iterated between the fluid system and the solid ODE until the system is solved. We consider the system as solved when the update of the ball velocity in an iteration is less than 10−8.

3. Fluid–structure interaction in fully Eulerian coordinates

We consider the full fluid–structure interaction problem including contact with Γbottom. We adopt here a fully Eulerian approach for the FSI system in order to enable the transition to contact.31,32,34 To allow for an implicit inclusion of the contact conditions into the system (see below), we use a Nitsche-based method for the FSI coupling as presented in Ref. 39.

a. Solid bilinear form and FSI coupling
The fluid bilinear form has already been detailed in (12). To introduce the solid form, we denote the reduced solid domain by S ( t ) and define the bilinear form corresponding to the linear elasticity equations in Eulerian coordinates (2) by
A s ( d , d ̇ ; w , z ) m s ( t d ̇ d ̇ d ̇ , w ) + a s ( d , w ) + m d ̇ ( t d d ̇ d + d ̇ , z ) ,
where
m s ( d , w ) 2 π ρ s S r d w d x , a s ( d , w ) 2 π S r σ s ( d ) : w + σ s , r w r d x , m d ̇ ( d , w ) 2 π S d w d x ,
and
σ s = 2 μ s E ( d ) + λ s tr ( E ( d ) ) + 1 r d r I , σ s , r = 2 μ s + λ s r d r + λ s tr ( E ( d ) ) .
Moreover, we make use of the Nitsche terms defined in Subsection IV A 2 to impose the FSI coupling conditions (3),
n ( u , p , d ̇ ; v , q , w ) n s ( u d ̇ , v w ) + n c ( u , v w ) + n c ( v , u d ̇ ) + n p ( v w , p ) n p ( u d ̇ , q ) .
Note the negative sign in front of the last term, which is required to ensure the stability of the FSI formulation; see Ref. 25.
b. Discretization and stabilization

In order to resolve the interface I within the discretization, we use the locally modified finite-element method introduced in Ref. 58. This fitted finite-element method is based on a coarse unfitted patch mesh, which is independent of the interface location. The coarse cells are then divided in such a way into sub-triangles and sub-quadrilaterals that the interface is resolved in a linear approximation. In this work, we use equal-order locally modified finite elements of first order, in combination with an anisotropic edge-oriented pressure stabilization term sp(p, q); see Ref. 59 for details. We denote the locally modified finite-element space of first order by X h 1 . The discrete spaces for fluid velocity and pressure Vh and Qh and for solid displacement and velocity Wh and Zh are given by applying the respective Dirichlet conditions to the locally modified finite-element space X h 1 and by restricting the degrees of freedom to the fluid or solid sub-domain, respectively.

In addition, we add a stabilization term s d ( d ̇ , z ) of artificial diffusion type to the solid equations (see Ref. 32) and a consistent stabilization at the boundary Γsym corresponding to the rotational axis to ensure that the second boundary condition (11) is accurately imposed in the discrete formulation,
s d ( d ̇ , z ) = α d h 2 ( r d ̇ , z ) S , s r ( u , v ) α sym ρ f μ f Γ sym h n 2 h τ r u z r d S .
Here, hn and hτ refer to the cell-sizes in the normal and tangential directions, respectively. We summarize the stabilization terms in the bilinear form,
s ( u , p , d ̇ ; v , q , z ) = s p ( p , q ) + s d ( d ̇ , z ) + s r ( u , v ) .
For time discretization, we use the modified dG(0) time discretization presented in Ref. 60, which can be seen as a variant of the dG(0)/backward Euler methods that considers the movement of the interface in each space–time slab. The interface positions and the domain affiliations are updated explicitly based on the displacement d(tn−1) of the previous time step and using the initial point set/backward characteristics method.31,32 This means that we set
I n I ( d ( t n 1 ) ) , S n S ( d ( t n 1 ) ) , F n F ( d ( t n 1 ) ) .
c. Contact treatment

When a part I c I Γ bottom of I enters into contact with Γbottom, the FSI conditions need to be substituted with appropriate contact conditions. It has been noted in Refs. 39 and 40 that although the fluid layer between the ball and the lower wall vanishes (from a macroscopical perspective), an extension of the fluid forces to the contact surface I c has to be considered to obtain a physically relevant contact formulation. Here, we use the simplest possible numerical approach, which is to relax the no-penetration condition by a small ϵ = ϵ(h) > 0 such that a very thin mesh-dependent fluid layer remains at all times.

The distance to Γbottom depends on the (Eulerian) displacement d(tn) in the following way:
dist I ( t n ) , Γ bottom = dist I ( t n 1 ) , Γ bottom + d n ( t n ) d n ( t n 1 ) .
The contact conditions, relaxed by a small ϵ = ϵ(h) > 0, read as
dist I ( t n ) , Γ bottom ϵ , σ n 0 , σ n dist I ( t n ) , Γ bottom ϵ = 0  on  I n = I ( t n 1 ) ,
(14)
where
σ n n T σ s n n T σ f n .
In other words, the relaxation means that the contact conditions are already applied at an ϵ-distance from Γbottom. The three conditions (14) can be equivalently formulated in equality form with an arbitrary γC > 0,61 
σ n = γ C [ dist I ( t n ) , Γ bottom ϵ γ C 1 σ n = : P γ ( d , σ n ) ] +  on  I n ,
(15)
where [·]+ stands for the positive part of (·). The contact parameter γC will be chosen as γ C = γ C 0 λ s h 1 and the relaxation parameter will be chosen as ϵ = ϵ0hy, where hy denotes the cell size in the vertical direction at the bottom of the cylinder.

Note that (15) includes both the FSI coupling and the contact condition as in the absence of contact, it is exactly the FSI interface condition in the normal direction. For this reason, the transition between FSI coupling and contact conditions can be included easily in a fully implicit fashion in the variational formulation.

The final variational formulation reads as follows:
Find  ( u , p , d , d ̇ ) V h × Q h × W h × Z h  such that for all  ( v , q , w , z ) V h × Q h × W h × Z h ,  it holds A f ( u , p ; v , q ) + A s ( d , d ̇ ; w , z ) + n ( u , p , d ̇ ; v , q , w ) + s ( u , p , d ̇ ; v , q , z ) + r [ P γ ( d , σ n ) ] + , w n I = ( r f s , w ) Ω .
d. Implementation

The described algorithms and equations have been implemented in the finite-element library Gascoigne3D.62 We use a Cartesian finite-element mesh, which is highly refined in the region where contact occurs.

Concerning time discretization, we start with a relatively coarse time step Δt = 2 · 10−3, which captures the essential dynamics of the case process. When the ball gets close to the lower wall, the time step is reduced in order to capture the contact dynamics and, in particular, to resolve the impact time accurately. We do this by reducing the time step by a factor of two each time the distance to the wall drops below certain thresholds di, i = 0, …, m. In the Rubber22 case, we choose, for example, d0 = 10−2 m and d1 = 10−3 m. For the PTFE6 case, the contact interval is much shorter (around 10−4 s compared to 4 · 10−3 s for rubber). For this reason, we specify seven thresholds to reduce the time step in seven steps until Δt < 10−5 s.

In addition, we use the following numerical parameters:
γ C 0 = 1 , σ = 1 0 5 , α d = 1 , α sym = 1 0 3 .
The contact relaxation parameter is chosen as ϵ 0 = 1 8 for the Rubber22 test case and ϵ 0 = 1 4 for PTFE6. For a detailed sensitivity study of the influence of the contact parameters, we refer to Ref. 39.

In the following, we shall abbreviate the methods described in Subsection IV A 1 as ALE, in Subsection IV A 2 as CutFEM, and in Subsection IV A 3 as FSI.

Remark.

In order to indicate the computational effort needed to solve the (linearized) systems resulting from our computations, we give the size and sparsity of these with our quantitative results. The number fo free degrees of freedom is abbreviated as dof, and the number of non-zero entries in the system is denoted by nze.

In the CutFEM method, dof and nze can vary between time steps since different elements are active in different time steps. We therefore state the number of unconstrained degrees of freedom and non-zero entries of the linearized system in the first time step.

1. PTFE6

Our quantitative results for the PTFE6 setup are presented in Table III. Figure 3 shows the distance between the bottom of the ball to the bottom of the fluid domain over time from the experimental data and all three numerical methods.

FIG. 3.

The distance between the bottom of the ball and the bottom of the tank: Experimental and numerical results for the PFTE6 setup. The experimental data are reproduced with permission from T. Hagemeier, Particle settling-transitional regime, version 1, Mendeley Data, September 2020. Copyright 2020 Author(s), licensed under a Creative Commons Attribution 4.0 License.

FIG. 3.

The distance between the bottom of the ball and the bottom of the tank: Experimental and numerical results for the PFTE6 setup. The experimental data are reproduced with permission from T. Hagemeier, Particle settling-transitional regime, version 1, Mendeley Data, September 2020. Copyright 2020 Author(s), licensed under a Creative Commons Attribution 4.0 License.

Close modal

If we look at the pre-contact quantities of interest in Table III, we see that all methods give very similar results for the given quantities. Looking at the velocity of the PTFE6 ball, we see that the numerical values are within a relative error of 5.1% of the experiment on the finest discretizations. Taking into account that these results ignore the 2 mm deflection from the z axis observed in the experimental data, we consider this to be acceptable. Since the elastic effects of the particle appear to be negligible in this phase of the problem and due to the known good approximation properties of the ALE method, we consider these to be the most accurate values for future comparison.

Looking at the quantities of interest in the later phase, we see that in both the CutFEM and FSI methods, contact occurs later than in the experiment. This is consistent with the smaller speed of the particle compared to the experiments as observed above. With respect to the jump, we see that both methods capture the rebound dynamics since both the point in time at which the peak of the rebound is realized and the size of the jump are consistent with the experiment. As the CutFEM contact parameters were tuned with respect to the size of this jump, this is unsurprising. However, since the contact force only acts for a very small number of time steps, the fact that the time at which the rebound is maximal is also captured well shows that even after contact, the system is still approximated well by the fluid–rigid body system. Nevertheless, it is clear that the FSI system captures the dynamics much more accurately and without the need of essentially unknowable and artificial parameters dist0 and γc in the contact model.

2. Rubber22

Our results for the Rubber22 problem are presented in Table IV. The height of the ball over time can be seen in Fig. 4. To illustrate the applicability of our spatially reduced formulation, we solve the fluid–rigid body system in full three spatial dimensions using the ALE approach. The ALE approach has been chosen for this due to its significantly higher computational performance compared to the CutFEM and FSI methods. These results are also given in Table IV. To illustrate the flow solution, we show the result of the CutFEM simulation on the coarsest mesh, rotated into the xz plane at t = 0.575 in Fig. 5.

FIG. 4.

The distance between the bottom of the ball and the bottom of the tank: experimental and numerical results for the Rubber22 setup. The experimental data are reproduced with permission from T. Hagemeier, Particle settling-transitional regime, version 1, Mendeley Data, September 2020. Copyright 2020 Author(s), licensed under a Creative Commons Attribution 4.0 License.

FIG. 4.

The distance between the bottom of the ball and the bottom of the tank: experimental and numerical results for the Rubber22 setup. The experimental data are reproduced with permission from T. Hagemeier, Particle settling-transitional regime, version 1, Mendeley Data, September 2020. Copyright 2020 Author(s), licensed under a Creative Commons Attribution 4.0 License.

Close modal
FIG. 5.

Velocity solution (left) and pressure solution (center) at t = 0.575 s rotated into the xz plane and computational mesh (right) for the Rubber22 case, resulting from the CutFEM computation with hmax = 0.008 and Δt = 0.005 s.

FIG. 5.

Velocity solution (left) and pressure solution (center) at t = 0.575 s rotated into the xz plane and computational mesh (right) for the Rubber22 case, resulting from the CutFEM computation with hmax = 0.008 and Δt = 0.005 s.

Close modal

We again start by inspecting the pre-contact results in Table IV. We observe that t*, v*, and f* are again very similar for all methods. However, here, the discrepancy between the numeric and experimental velocity values is significantly smaller, with the relative difference being between 2.1% and 0.6% on the finest discretizations. We note at this point that the deviation from the z axis in this experiment was less than 0.75 mm. This shows that even for a more elastic material, a fluid–rigid body system can capture the pre-contact dynamics as well as a full FSI model. Furthermore, we see that with perfect initial data and spatial symmetry, our spatially reduced model, derived under the assumption of a rotationally symmetric solution, captures the dynamics as well as the significantly more computationally expensive full three-dimensional computation.

Looking at the numbers for the contact and rebound dynamics, we again see that the time of contact is similar for both the CutFEM and FSI methods, and this matches the experimental time with an error of less than 5%. For the CutFEM method, we clearly see that the PTFE tuned parameters are not able to capture the rebound dynamics well. In fact, the size of the rebound is approximately three times larger than the physical rebound. This shows that the use of an artificial lubrication force, as considered in a variety of other literature, can lead to physically meaningful results but is heavily dependent on the “correct” choice of parameters for which there is no a priori knowledge. For the FSI model, we see that the overall dynamics are captured well. We also observe that—while the pre-contact dynamics are essentially independent of the elasticity parameters—a variation of the elasticity modulus Es changes the rebound height djump significantly. Here, a softer material (Es = 2 · 106 Pa) leads to a larger rebound as more elastic energy is taken up through the deformation during the impact.

3. Three-dimensional computation including rotational effects

One of the challenges in performing the experimental study42 was to limit the horizontal deflection of the falling objects, i.e., to keep them close to the centerline. Identifying the source of these three-dimensional effects is one of the intriguing questions for future research. The cause may be found in a complex solution pattern of the Navier–Stokes equations, in material inaccuracies such as non-uniform distribution of the mass or the surface roughness but also in experimental inaccuracies, e.g., during the release process.

Since the reduced formulation presented in Sec. III cannot depict 3D effects, we add further studies based on the full 3D formulation. Here, we choose the ALE representation again due to its much better computational performance, even if no contact can be modeled. However, this is not a significant limitation since effects of rotation and horizontal deflection happen during the phase of free fall. The three-dimensional simulations based on the ALE formulation presented above in Subsection IV B 2 did not show any three-dimensional effect if the configuration is fully symmetric and the ball is released at the centerline, i.e., c S = 0 , 0 , d ( 0 ) . To investigate the stability of the Navier–Stokes rigid-body system, we consider further numerical simulations based on the Rubber22 case with distorted initial values. We start the simulation with the initial data
c S = 1 0 3 χ x , 1 0 3 χ y , 0.1461203 m , v S ( 0 ) = 4 1 0 3 χ x , 4 1 0 3 χ y , 0 ms 1 , ω S ( 0 ) = χ 1 , χ 2 , χ 2 2 π 1 0 2 s 1 ,
where χ x , χ y , χ 1 , χ 2 , χ 3 iid N ( 0 , 1 ) are normally distributed random numbers with mean zero and standard deviation 1.

Figure 6 shows the results for multiple experiments based on randomly chosen initial data. The left plot shows the projection of the center of mass onto the xy plane. These results show that numerical simulations cannot predict a substantial deflection from the centerline if an initial deflection is prescribed. While the rigid solids are indeed further removed from the center, the effect is small and the objects remain within 3 mm off the center. On the right, we show the velocity of the particles. The upper figure shows the horizontal velocity component, while the lower plot gives the dominant vertical velocity. Here, we indeed see a substantial impact of the initial disturbances. When the solid comes close to the lower boundary, a deflection to the sides becomes visible. We note that these ALE simulations crash before contact is established. At final time t ≈ 0.64 s, this distance between the lower boundaries is still slightly larger than radius r S = 0.011 m . Since the horizontal velocity is beginning to increase significantly here, further simulations, in which the particle is able to get closer to the bottom boundary, are of interest. We deduce from these results that small fluctuations during the release process can indeed explain the small horizontal displacement observed during the experiments for the PTFE and rubber particles and the discrepancy between the experimental and numerical realizations of our quantities of interest. However, these computations also show that very large horizontal displacements, as observed in Refs. 41 and 42, for example, for the POM particle, cannot be solely explained through this.

FIG. 6.

Left: view from below. Deflection of the center of mass from the centerline (0, 0) for 15 experiments starting with random initial deviations each. Right: velocity of the particles for six experiments. The upper figure shows the horizontal velocity component v S , 1 2 + v S , 2 2 , and the lower plot gives the vertical velocity close to the bottom for t ≥ 0.5 s.

FIG. 6.

Left: view from below. Deflection of the center of mass from the centerline (0, 0) for 15 experiments starting with random initial deviations each. Right: velocity of the particles for six experiments. The upper figure shows the horizontal velocity component v S , 1 2 + v S , 2 2 , and the lower plot gives the vertical velocity close to the bottom for t ≥ 0.5 s.

Close modal

We have presented two setups for a fluid–structure interaction problem with solid contact based on the physical experiments described in Ref. 42 and the data available in Ref. 41. We computed these setups using a spatially reduced model under the assumption of rotational symmetry in cylindrical coordinates. For the discretization, we used a fitted ALE and unfitted CutFEM approach within a fluid–rigid body model and a fully Eulerian FSI approach in a fluid–elastic structure model capable of resolving the solid contact.

We showed how each of these discretizations is able to capture the pre-contact dynamics observed in the physical experiment within a margin of 5.1%–0.6%, even though this ignored any horizontal motion observed in the experimental data. Using a full three-dimensional ALE discretization, we saw that the spatially reduced approach did indeed result in meaningful results under the assumption of perfect initial conditions at a fraction of the computational cost. Furthermore, we presented computation with disturbed initial data. From these results, we deduced that the observed horizontal motion in the PTFE and rubber experiments is within the scope explainable by imperfect starting conditions. This shows that a fluid–rigid body model is suitable for this type of problem before solid contact occurs.

With respect to the contact dynamics, we were able to show that the Eulerian FSI discretization with contact treatment is able to reproduce the spatial and temporal dynamics observed in the experiments very well. This is even though the theoretical modeling of such contact dynamics is not yet fully understood. The moving domain CutFEM approach together with the artificial contact treatment showed that this type of contact treatment can result in physically meaningful results when the artificial parameters are chosen “correctly” and the extent to which the artificial parameters are material dependent.

The resulting datasets, the source code for the fluid–rigid body discretizations, as well as some simplified examples are available in the zenodo repository.63 

We conclude that the two discussed setups are well suited for the validation of fluid–structure interaction models in the moderate Reynolds-number fluid regime both before and after contact. To the best of our knowledge, this is the first example of a computational FSI setup with rebound contact dynamics, which is validated by experimental data. We note that the PTFE6 setup is better suited for validation of contact and rebound models since the model parameters are known precisely. On the other hand, the Rubber22 scenario is well suited to validate models before or without contact since there is less deviation from the centerline in the data from the experiment.

For future research, it remains to be investigated that to what extent the free surface at the top of the fluid domain plays a role in system dynamics. Furthermore, an open question is that of the role of imperfections in the mass distribution within the solid, of the surface roughness of the solid, and whether these can explain larger horizontal displacements and rotation observed, for example, with the POM particle in Ref. 42.

H.v.W. and T.R. acknowledge support from the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Grant No. 314838170, GRK 2297 MathCoRe. T.R. further acknowledges support from the Federal Ministry of Education and Research of Germany (Project No. 05M16NMA).

The numerical data supporting the finding of this study are openly available in Zenodo, Ref. 63. The experimental data, on which this study is based, is openly available in Mendeley Data, Ref. 41.

1. Fluid–rigid body interaction in arbitrary Lagrangian–Eulerian coordinates
The ALE formulation is based on the reference map
T A L E ( c S ) : F F ( c S ) ,
where c S is the solid's center of mass relative to the bottom boundary and where
F ( c S ) = [ 0 , 0.055 m ] × [ 0 , 0.2 m ] \ B r S c S ,
B r S ( c S ) being the open ball of radius r with midpoint cS. As the reference domain, we set F F ( 0.05 m ) where the ball is centered at c S = ( 0 , 0 , 0.5 ) such that the mesh distortion is not too extreme when the solid comes close to the lower boundary. The mapping T A L E ( c S ) is given analytically by
T A L E ( c S ; r , z ) = r , z + ( c S 0.05 m ) f A L E ( z ) , f A L E ( z ) = z 0.05 m 2 r S , z < 0.05 m 2 r S 1 , 0.05 m 2 r S z 0.05 m + 2 r S 0.2 m z 0.15 m + 2 r S , z > 0.05 m + 2 r S
(A1)
such that T A L E ( c S ; r , z ) is a pure translation for all z [ 0.05 m 2 r S , 0.05 m + 2 r S ] , where the reference ball is located. We make sure that the lines at z = 0.05 m ± 2 r S are resolved by the computational mesh such that TALE is differentiable within all elements. This construction of the domain map does not allow us to reduce the distance of the ball from the lower boundary to less than r S .
The reference domain F is the basis for the finite-element discretization. In cylindrical coordinates, the ALE version of the variational formulation takes the form
m A L E ( c S ; u , v ) = 2 π ρ f F r J ( c S ) u v d x ,
a A L E ( c S ; u , v ) = 2 π μ f F J ( c S ) u 1 0 0 J ( c S ) 1 : v × 1 0 0 J ( c S ) 1 + J ( c S ) r u r v r d x ,
c A L E ( c S ; u , v , w ) = 2 π ρ f F r J ( c S ) 1 0 0 J ( c S ) 1 × u t T A L E ( c S ) w d x ,
b A L E ( c S ; q , v ) = 2 π F J ( c S ) q v r + r r v 1 + z v 2 / J ( c S ) d x ,
where J ( c S ) = det ( T A L E ) is the determinant of the deformation gradient.
The discretization is realized by means of quadratic equal-order finite elements for velocity and pressure on a quadrilateral mesh of the reference domain F ; we refer to Ref. 12, Sec. 4.4 for details on the realization. To stabilize the inf–sup condition, we use the local projection scheme as introduced by Becker and Braack,64 given, in the ALE formulation, as
s A L E ( c S ; p , q ) = F r δ r 0 0 δ z κ h ( p ) J ( c S ) 0 0 J ( c S ) 1 × κ h ( q ) + J ( c S ) r δ r κ h ( p ) κ h ( q ) d x .
Here, κh: = id − i2h is the coarse mesh fluctuation operator that subtracts the interpolation to the mesh with double spacing and where δr and δz are local stabilization parameters that depend on the element diameter h and the time step size k,
δ r = 0.1 μ f ρ f h 2 + 1 k 1 , δ z = 0.1 μ f ρ f h 2 J ( c S ) 2 + 1 k 1 .
The different scaling in r- and z-directions reflects the anisotropy induced by the ALE transformation, cf. Ref. 65 or (Ref. 12, Sec. 5.3.3).

In time, we discretize the Navier–Stokes equations and the rigid-body problem with BDF2 time stepping in a decoupled iteration taking 10−8 as tolerance for the solid velocity and deformation update. The forces acting on the solid are evaluated by means of the Babuška–Miller trick66,67 such that quadratic finite elements let us expect fourth order convergence. Hence, each step of spatial refinement will be accompanied by two refinements of the time step.

a. Fluid–rigid body interaction in arbitrary Lagrangian–Eulerian coordinates in three dimensions
To clarify the role of fixed body rotations of the solid and deflections of the center of mass from the z axis, both of which were observed in the experimental analysis,41,42 we perform full three-dimensional simulations. The Navier–Stokes equations are formulated in Cartesian coordinates (1), and the motion of the rigid body is described by
d d t c S ( t ) = v S ( t ) , d d t v S ( t ) m S = 0 0 m S g vol ( S ) ρ f g + I σ ( u , p ) n d s , d d t ω S ( t ) I S + ω S ( t ) × I S ω S ( t ) = I x c S ( t ) × σ ( u , p ) n d s ,
(A2)
where v S ( t ) R 3 is the full velocity vector, c S ( t ) R 3 is the solid's center of mass, and ω S ( t ) R 3 is the rotational velocity vector; cf. Eqs. (4)–(6). Assuming a homogeneous distribution of the density, the moment of inertia is given by
I S = 8 15 π ρ S r S 5 Id .
It follows that the nonlinear term in the rotational ODE (A2) vanishes since ω S ( t ) × I S ω S ( t ) = 8 15 π ρ S r S 5 ω S ( t ) × ω S ( t ) = 0 . On the surface of the sphere, the Navier–Stokes Dirichlet conditions are used to describe the velocity
v ( x , t ) = v S ( t ) + ω S ( t ) × x c S  on  I .
To evaluate the torque correctly through the variational formulation, i.e., by using the Babuška–Miller trick,66,67 the symmetric gradient is used instead in (10); hence,
a 3 d ( u , v ) = μ f Ω 3 d ( u + u T ) : v d x .
We cast the problem in standard ALE formulation and refer to Ref. 12, Chap. 5 for further details. The ALE map TALE(t) is chosen similar to the reduced case (A1), but we must incorporate motion in the xy plane and define
T A L E ( c S ; x , y , z ) = x + c S , 1 g A L E r ( x , y ) y + c S , 2 g A L E r ( x , y ) z + ( c S , 3 r S 0.05 m ) f A L E ( z ) ,
where fALE(·) is defined in (A1), while r(·, ·) and gALE(·) are given by
r ( x , y ) = max 0 , min 1 , x r S R r S , g A L E ( r ) = 1 1 1 + exp 1 2 r r r 2 .
The function r(·, ·) maps the x–y plane to [0, 1] with r(x, y) = 0 for x 2 + y 2 r S and r(x, y) = 1 for x 2 + y 2 R . Furthermore, gALE is a smooth transition function mapping (0, 1) to (0, 1), with all derivatives being zero at 0 and 1. On the other hand, the derivative of fALE(z) is not defined at z = 0.05 m ± 2 r S . These two surfaces are however resolved by the finite-element mesh to give good approximation characteristics of the ALE formulation.
b. Implementation
Both formulations, the reduced two-dimensional ALE formulation in cylindrical coordinates and the full three-dimensional ALE formulation, are implemented in the finite-element library Gascoigne3D.62 The coupling between the Navier–Stokes equation and the rigid-body motion is resolved in a simple iteration until the discrepancy in velocity reached a threshold
v I v S < 1 0 8 .
The nonlinearity of the Navier–Stokes equation is solved by a Newton iteration, and the resulting linear systems of equations are approximated with a parallel generalized minimal residual method (GMRES) iteration, preconditioned by a geometric multigrid solver; see Ref. 68. The meshes are graded with a higher resolution close to the solid.

We define two simplified test cases. This is intended to make it easier for others to reproduce the presented results using different methods and/or implementations.

1. Stationary flow test

For this stationary test, we modify the Rubber22 setup. The sphere is fixed at c S = ( 0 , 0 , 0.1 ) , i.e., the center of the cylinder. We impose an inflow boundary condition u = 0.01 ( 1 ( x 1 2 + x 2 2 ) / R 2 ) on Γtop, no-slip u = 0 on Γ wall I , and a homogeneous Neumann boundary condition σ(u, p)n = 0 on Γbottom.

We consider the stationary Navier–Stokes problem on this domain. As reference quantities, we take the vertical stress acting on the sphere, i.e., testing the reduced formulation in (13) with the non-conforming, continuous test-functions w ̂ = ( 0 , 1 ) T on I and 0 on Γ, respectively.

We compute the problem based on the discretizations discussed in Subsections IV A 1 and IV A 2. The results can be seen in Table V.

TABLE V.

Resulting reference quantities for the stationary test scenario.

Discretization Results
Method [hmin, hmax] dof nze Fz
ALE  [0.000 30, 0.004 9]  6.9  0.310  −4.432 34 · 10−5 
  [0.000 15, 0.002 4]  26.7  1.238  −4.429 92 · 10−5 
  [0.000 08, 0.001 2]  104.9  4.946  −4.429 75 · 10−5 
    Extrapolate  −4.429 74 · 10−5 
    order (in h)  3.8 
CutFEM  [0.002 0, 0.008]  11.9  0.361  −4.402 54 · 10−5 
  [0.001 0, 0.004]  41.3  1.241  −4.423 55 · 10−5 
  [0.000 5, 0.002]  151.4  4.525  −4.429 19 · 10−5 
Discretization Results
Method [hmin, hmax] dof nze Fz
ALE  [0.000 30, 0.004 9]  6.9  0.310  −4.432 34 · 10−5 
  [0.000 15, 0.002 4]  26.7  1.238  −4.429 92 · 10−5 
  [0.000 08, 0.001 2]  104.9  4.946  −4.429 75 · 10−5 
    Extrapolate  −4.429 74 · 10−5 
    order (in h)  3.8 
CutFEM  [0.002 0, 0.008]  11.9  0.361  −4.402 54 · 10−5 
  [0.001 0, 0.004]  41.3  1.241  −4.423 55 · 10−5 
  [0.000 5, 0.002]  151.4  4.525  −4.429 19 · 10−5 
2. Non-stationary flow test with prescribed motion
As a second test case, we keep the Rubber22 setup as the basis. The material parameters and the cylinder boundary conditions are as before, i.e., we consider no-slip on Γwall ∪ Γbottom and a free-slip condition on Γtop. We prescribe the motion of the sphere as follows: Over the time interval I = [0, 20], the sphere is located at
c S = 0 , 0 , d ( t ) with d ( t ) = 0.1 + 0.05 cos ( 0.1 π t ) .
The boundary condition on the interface is accordingly set to u = 0 , 0 , t d ( t ) . Quantities of interest are the maximal value over time of the z-component of the force functional in the reduced formulation.

We compute this again using the reduced formulation with the rigid-body discretizations. The quantitative results can be seen in Table VI, while the force functional is shown over time in Fig. 7.

TABLE VI.

Resulting reference quantities for the non-stationary moving domain test scenario.

Discretization Results
Method [hmin, hmax] Δt dof nze Fz,max tz,max
ALE  [0.000 30, 0.009 9]  1 / 5   11.1  0.49  1.018 38 · 10−4  4.107 346 91 
  [0.000 15, 0.005 0]  1 / 20   25.0  1.10  1.017 48 · 10−4  4.113 563 25 
  [0.000 08, 0.002 5]  1 / 80   58.7  2.61  1.017 20 · 10−4  4.106 704 36 
CutFEM  [0.002 0, 0.008]  1 / 25   11.9  0.361  1.017 00 · 10−4  4.153 
  [0.001 0, 0.004]  1 / 50   41.2  1.240  1.013 95 · 10−4  4.138 
  [0.000 5, 0.002]  1 / 100   151.3  4.526  1.016 71 · 10−4  4.111 
Discretization Results
Method [hmin, hmax] Δt dof nze Fz,max tz,max
ALE  [0.000 30, 0.009 9]  1 / 5   11.1  0.49  1.018 38 · 10−4  4.107 346 91 
  [0.000 15, 0.005 0]  1 / 20   25.0  1.10  1.017 48 · 10−4  4.113 563 25 
  [0.000 08, 0.002 5]  1 / 80   58.7  2.61  1.017 20 · 10−4  4.106 704 36 
CutFEM  [0.002 0, 0.008]  1 / 25   11.9  0.361  1.017 00 · 10−4  4.153 
  [0.001 0, 0.004]  1 / 50   41.2  1.240  1.013 95 · 10−4  4.138 
  [0.000 5, 0.002]  1 / 100   151.3  4.526  1.016 71 · 10−4  4.111 
FIG. 7.

Force functionals acting on the sphere with prescribed motion over time.

FIG. 7.

Force functionals acting on the sphere with prescribed motion over time.

Close modal
1.
G.
Závodszky
,
B.
van Rooij
,
B.
Czaja
,
V.
Azizi
,
D.
de Kanter
, and
A. G.
Hoekstra
, “
Red blood cell and platelet diffusivity and margination in the presence of cross-stream gradients in blood flows
,”
Phys. Fluids
31
,
031903
(
2019
).
2.
J. O.
Marston
,
J. P. K.
Seville
,
Y.-V.
Cheun
,
A.
Ingram
,
S. P.
Decent
, and
M. J. H.
Simmons
, “
Effect of packing fraction on granular jetting from solid sphere entry into aerated and fluidized beds
,”
Phys. Fluids
20
,
023301
(
2008
).
3.
H.
Brenner
, “
The slow motion of a sphere through a viscous fluid towards a plane surface
,”
Chem. Eng. Sci.
16
,
242
251
(
1961
).
4.
E.
Feireisl
, “
On the motion of rigid bodies in a viscous incompressible fluid
,” in
Nonlinear Evolution Equations and Related Topics
(
Birkhäuser
,
Basel
,
2003
), pp.
419
441
.
5.
D.
Gérard-Varet
,
M.
Hillairet
, and
C.
Wang
, “
The influence of boundary conditions on the contact problem in a 3D Navier–Stokes flow
,”
J. Math. Pures Appl.
103
,
1
38
(
2015
).
6.
D.
Gérard-Varet
and
M.
Hillairet
, “
Existence of weak solutions up to collision for viscous fluid-solid systems with slip
,”
Commun. Pure Appl. Math.
67
,
2022
2076
(
2014
).
7.
E.
Feireisl
, “
On the motion of rigid bodies in a viscous compressible fluid
,”
Arch. Ration. Mech. Anal.
167
,
281
308
(
2003
).
8.
R. H.
Davis
,
J.-M.
Serayssol
, and
E. J.
Hinch
, “
The elastohydrodynamic collision of two spheres
,”
J. Fluid Mech.
163
,
479
497
(
1986
).
9.
C.
Grandmont
and
M.
Hillairet
, “
Existence of global strong solutions to a beam-fluid interaction system
,”
Arch. Ration. Mech. Anal.
220
,
1283
1333
(
2016
).
10.
M.
Hillairet
and
T.
Takahashi
, “
Collisions in three-dimensional fluid structure interaction problems
,”
SIAM J. Math. Anal.
40
,
2451
2477
(
2009
).
11.
G.
Gravina
,
S.
Schwarzacher
,
O.
Souček
, and
K.
Tůma
, “
Contactless rebound of elastic bodies in a viscous incompressible fluid
,” arXiv:2011.01932v1 (
2020
).
12.
T.
Richter
,
Fluid-Structure Interactions: Models, Analysis and Finite Elements
, Lecture Notes in Computational Science and Engineering 118 (
Springer
,
2017
).
13.
J.
Donea
,
S.
Giuliani
, and
J. P.
Halleux
, “
An arbitrary Lagrangian-Eulerian finite element method for transient dynamic fluid-structure interactions
,”
Comput. Methods Appl. Mech. Eng.
33
,
689
723
(
1982
).
14.
S.
Frei
,
T.
Richter
, and
T.
Wick
, “
Long-term simulation of large deformation, mechano-chemical fluid-structure interactions in ALE and fully Eulerian coordinates
,”
J. Comput. Phys.
321
,
874
891
(
2016
).
15.
T.
Rezaee
and
K.
Sadeghy
, “
Effect of porosity on the settling behavior of a 2D elliptic particle in a narrow vessel: A lattice-Boltzmann simulation
,”
Phys. Fluids
31
,
123301
(
2019
).
16.
A.
Johansson
,
M. G.
Larson
, and
A.
Logg
, “
High order cut finite element methods for the Stokes problem
,”
Adv. Model. Simul. Eng. Sci.
2
,
24
(
2015
).
17.
B.
Schott
,
C.
Ager
, and
W. A.
Wall
, “
A monolithic approach to fluid‐structure interaction based on a hybrid Eulerian‐ALE fluid domain decomposition involving cut elements
,”
Int. J. Numer. Methods Eng.
119
,
208
237
(
2019
).
18.
R.
Glowinski
,
T.-W.
Pan
,
T. I.
Hesla
, and
D. D.
Joseph
, “
A distributed Lagrange multiplier/fictitious domain method for particulate flows
,”
Int. J. Multiphase Flow
25
,
755
794
(
1999
).
19.
S.
Court
and
M.
Fournié
, “
A fictitious domain finite element method for simulations of fluid–structure interactions: The Navier–Stokes equations coupled with a moving solid
,”
J. Fluids Struct.
55
,
398
408
(
2015
).
20.
S.
Court
,
M.
Fournié
, and
A.
Lozinski
, “
A fictitious domain approach for the Stokes problem based on the extended finite element method
,”
Int. J. Numer. Methods Fluids
74
,
73
99
(
2013
).
21.
E.
Burman
,
S.
Claus
,
P.
Hansbo
,
M. G.
Larson
, and
A.
Massing
, “
CutFEM: Discretizing geometry and partial differential equations
,”
Int. J. Numer. Methods Eng.
104
,
472
501
(
2014
).
22.
E.
Burman
,
S.
Frei
, and
A.
Massing
, “
Eulerian time-stepping schemes for the non-stationary Stokes equations on time-dependent domains
,” arXiv:1910.03054v1 (
2019
).
23.
H.
von Wahl
,
T.
Richter
, and
C.
Lehrenfeld
, “
An unfitted Eulerian finite element method for the time-dependent Stokes problem on moving domains
,” arXiv:2002.02352v1 (
2020
).
24.
C.
Lehrenfeld
and
M.
Olshanskii
, “
An Eulerian finite element method for PDEs in time-dependent domains
,”
ESAIM: Math. Model. Numer. Anal.
53
,
585
614
(
2019
).
25.
E.
Burman
and
M. A.
Fernández
, “
An unfitted Nitsche method for incompressible fluid-structure interaction using overlapping meshes
,”
Comput. Methods Appl. Mech. Eng.
279
,
497
514
(
2014
).
26.
A.
Gerstenberger
and
W. A.
Wall
, “
An eXtended finite element method/Lagrange multiplier based approach for fluid–structure interaction
,”
Comput. Methods Appl. Mech. Eng.
197
,
1699
1714
(
2008
).
27.
P.
Hansbo
and
J.
Hermansson
, “
Nitsche's method for coupling non-matching meshes in fluid-structure vibration problems
,”
Comput. Mech.
32
,
134
139
(
2003
).
28.
A.
Legay
,
J.
Chessa
, and
T.
Belytschko
, “
An Eulerian–Lagrangian method for fluid–structure interaction based on level sets
,”
Comput. Methods Appl. Mech. Eng.
195
,
2070
2087
(
2006
).
29.
F. P. T.
Baaijens
, “
A fictitious domain/mortar element method for fluid–structure interaction
,”
Int. J. Numer. Methods Fluids
35
,
743
761
(
2001
).
30.
G.-H.
Cottet
,
E.
Maitre
, and
T.
Milcent
, “
Eulerian formulation and level set models for incompressible fluid-structure interaction
,”
ESAIM: Math. Model. Numer. Anal.
42
,
471
492
(
2008
).
31.
T.
Dunne
, “
An Eulerian approach to fluid–structure interaction and goal-oriented mesh adaptation
,”
Int. J. Numer. Methods Fluids
51
,
1017
1039
(
2006
).
32.
S.
Frei
, “
Eulerian finite element methods for interface problems and fluid-structure interactions
,” Ph.D. thesis,
Ruprecht-Karls-Universität Heidelberg
,
2016
.
33.
F.
Hecht
and
O.
Pironneau
, “
An energy stable monolithic Eulerian fluid–structure finite element method
,”
Int. J. Numer. Methods Fluids
85
,
430
446
(
2017
).
34.
T.
Richter
, “
A fully Eulerian formulation for fluid–structure-interaction problems
,”
J. Comput. Phys.
233
,
227
240
(
2013
).
35.
S.
Turek
and
J.
Hron
, “
Proposal for numerical benchmarking of fluid–structure interaction between an elastic object and laminar incompressible flow
,” in Fluid-Structure Interaction, Lecture Notes in Computational Science and Engineering Vol. 53 (
Springer Berlin Heidelberg
,
2006
), pp.
371
385
.
36.
C.
Ager
,
B.
Schott
,
A. T.
Vuong
,
A.
Popp
, and
W. A.
Wall
, “
A consistent approach for fluid‐structure‐contact interaction based on a porous flow model for rough surface contact
,”
Int. J. Numer. Methods Eng.
119
,
1345
1378
(
2019
).
37.
C.
Ager
,
A.
Seitz
, and
W. A.
Wall
, “A consistent and versatile computational approach for general fluid-structure-contact interaction problems,”
Int. J. Numer. Methods Eng.
(published online) (
2020
).
38.
E.
Burman
,
M. A.
Fernández
,
S.
Frei
, and
F. M.
Gerosa
, “
3D-2D Stokes-Darcy coupling for the modelling of seepage with an application to fluid-structure interaction with contact
,” arXiv:1912.08503 (
2019
).
39.
E.
Burman
,
M. A.
Fernández
, and
S.
Frei
, “
A Nitsche-based formulation for fluid-structure interactions with contact
,”
ESAIM: Math. Model. Numer. Anal.
54
,
531
564
(
2020
).
40.
S.
Zonca
,
P. F.
Antonietti
, and
C.
Vergara
, “
A polygonal discontinuous Galerkin formulation for contact mechanics in fluid-structure interaction problems
,” MOX-Report No. 26/2020,
2020
.
41.
T.
Hagemeier
, Particle settling-transitional regime, version 1,
Mendeley Data
, September
2020
.
42.
T.
Hagemeier
,
D.
Thévenin
, and
T.
Richter
, “Settling of spherical particles in the transitional regime,”
Int. J. Multiph. Flow
138
,
103589
(
2021
).
43.
V.
John
,
A.
Linke
,
C.
Merdon
,
M.
Neilan
, and
L. G.
Rebholz
, “
On the divergence constraint in mixed finite element methods for incompressible flows
,”
SIAM Rev.
59
,
492
544
(
2017
).
44.
M. L.
Anderson
,
P. H.
Mott
, and
C. M.
Roland
, “
The compression of bonded rubber disks
,”
Rubber Chem. Technol.
77
,
293
302
(
2004
).
45.
C.
Bernardi
,
M.
Dauge
, and
Y.
Maday
,
Spectral Methods for Axisymetric Domain
, Series in Applied Mathematics (
Gauthier-Villars
,
Paris
,
1999
).
46.
T. J. R.
Hughes
,
W. K.
Liu
, and
T. K.
Zimmermann
, “
Lagrangian-Eulerian finite element formulation for incompressible viscous flows
,”
Comput. Methods Appl. Mech. Eng.
29
,
329
349
(
1981
).
47.
T.
Richter
and
T.
Wick
, “
Finite elements for fluid–structure interaction in ALE and fully Eulerian coordinates
,”
Comput. Methods Appl. Mech. Eng.
199
,
2633
2642
(
2010
).
48.
J.
Guzmán
and
M.
Olshanskii
, “
Inf-sup stability of geometrically unfitted Stokes finite elements
,”
Math. Comput.
87
,
2091
2112
(
2017
).
49.
E.
Burman
, “
Ghost penalty
,”
C. R. Math.
348
,
1217
1220
(
2010
).
50.
C.
Lehrenfeld
, “
High order unfitted finite element methods on level set domains using isoparametric mappings
,”
Comput. Methods Appl. Mech. Eng.
300
,
716
733
(
2016
).
51.
J.
Nitsche
, “
Über ein variationsprinzip zur lösung von Dirichlet-problemen bei verwendung von teilräumen, die keinen randbedingungen unterworfen sind
,”
Abh. Math. Semin. Univ. Hambg.
36
,
9
15
(
1971
).
52.
J.
Preuss
, “
Higher order unfitted isoparametric space-time FEM on moving domains
,” MA thesis,
Georg-August Universität Göttingen
,
2018
.
53.
S.
Sathe
and
T. E.
Tezduyar
, “
Modeling of fluid-structure interactions with the space-time finite elements: Contact problems
,”
Comput. Mech.
43
,
51
60
(
2008
).
54.
D.
Wan
and
S.
Turek
, “
Direct numerical simulation of particulate flow via multigrid FEM techniques and the fictitious boundary method
,”
Int. J. Numer. Methods Fluids
51
,
531
566
(
2006
).
55.
J.
Schöberl
, “
NETGEN an advancing front 2D/3D-mesh generator based on abstract rules
,”
Comput. Vis. Sci.
1
,
41
52
(
1997
).
56.
J.
Schöberl
, “
C++11 implementation of finite elements in NGSolve
,” ASC Technical Report Report No. 30/2014,
Institute for Analysis and Scientific Computing, TU Wien
, September
2014
.
57.
Ngsxfem: An add-on to NGSolve for unfitted finite element discretizations, version 1.3.2008,
2020
.
58.
S.
Frei
and
T.
Richter
, “
A locally modified parametric finite element method for interface problems
,”
SIAM J. Numer. Anal.
52
,
2315
2334
(
2014
).
59.
S.
Frei
, “
An edge‐based pressure stabilization technique for finite elements on arbitrarily anisotropic meshes
,”
Int. J. Numer. Methods Fluids
89
,
407
429
(
2019
).
60.
S.
Frei
and
T.
Richter
, “
A second order time-stepping scheme for parabolic interface problems with moving interfaces
,”
ESAIM: Math. Model. Numer. Anal.
51
,
1539
1560
(
2017
).
61.
P.
Alart
and
A.
Curnier
, “
A mixed formulation for frictional contact problems prone to Newton like solution methods
,”
Comput. Methods Appl. Mech. Eng.
92
,
353
375
(
1991
).
62.
R.
Becker
,
M.
Braack
,
D.
Meidner
, and
T.
Richter
,
The Finite Element Toolkit Gascoigne 3D
, www.gascoigne.de (accessed February 24, 2021).
63.
H.
vonWahl
,
T.
Richter
,
S.
Frei
, and
T.
Hagemeier
(
2020
), “
Falling balls in a viscous fluid with contact: Comparing numerical simulations with experimental data
,” Zenodo. https://doi.org/10.5281/zenodo.3989604.
64.
R.
Becker
and
M.
Braack
, “
A finite element pressure gradient stabilization for the Stokes equations based on local projections
,”
Calcolo
38
,
173
199
(
2001
).
65.
M.
Braack
and
T.
Richter
, “
Local projection stabilization for the Stokes system on anisotropic quadrilateral meshes
,” in
Proceedings of ENUMATH 2005 the 6th European Conference on Numerical Mathematics and Advanced Applications, Santiago de Compostela, Spain, July 2005
, edited by
A.
de Castro
,
D.
Gómez
,
P.
Quintela
, and
P.
Salgado
,
2006
.
66.
I.
Babuška
and
A.
Miller
, “
The post-processing approach in the finite element method. I: Calculations of displacements, stresses and other higher derivatives of the displacements
,”
Int. J. Numer. Methods Eng.
20
,
1085
1109
(
1984
).
67.
H.
von Wahl
,
T.
Richter
,
C.
Lehrenfeld
,
J.
Heiland
, and
P.
Minakowski
, “
Numerical benchmarking of fluid-rigid body interactions
,”
Comput. Fluids
193
,
104290
(
2019
).
68.
L.
Failer
and
T.
Richter
, “
A parallel Newton multigrid framework for monolithic fluid-structure interactions
,”
J. Sci. Comput.
82
,
28
(
2020
).
Published open access through an agreement with Otto-von-Guericke-Universitat Magdeburg Fakultat fur Mathematik163504