We evaluate a number of different finiteelement 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 finiteelement 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 wellestablished 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 FSIcontact 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 threedimensional 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 opensource finite element libraries, and the results are made freely available to aid reproducibility.
I. INTRODUCTION
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 flows^{1} 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 nonlinear incompressible Navier–Stokes together with noslip 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 freeslip condition,^{5} the rough nature of the surface is taken into account^{6} 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 wellestablished 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 latticeBoltzmann 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 approaches^{21} 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 subdomains are meshed separately. The two meshes are then coupled by means of Nitsche's method^{25,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 fluidstructure interaction problems such as those in Ref. 35 completely avoid contact since the methods that handle contact remain relatively new. For rigidbody 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 rigidbody ALE, a rigidbody 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 threedimensional computations. To the best of our knowledge, there is currently no comparable benchmark that considers such a multiphase flow/fluidstructure 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.
II. DESCRIPTION
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.
A. Physical experiment
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 highspeed CMOScamera. 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 inplane 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.
B. Mathematical model
We consider a bounded domain $\Omega \u2208 R d $, with d ∈ {2, 3}, over a finite, nonempty, time interval [0, T_{end}]. This is divided into a ddimensional fluid region $F$, a ddimensional solid $S$, and a d − 1dimensional interface $I$ dividing the solid and fluid regions. For these, we have $\Omega =F \u222a \u0307 I \u222a \u0307 S$.
1. Fluid model
2. Elastic solid and fluid–structure interaction
3. Rigid solid
C. Domain description
Since we consider balls of different sizes, we shall keep the domain description general. The background domain is a cylinder $\Omega = { x = ( x , y , z ) T \u2208 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 \u2208 R 3  x 2 + y 2 + ( z \u2212 ( 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 h_{0}. Accordingly, the volume of the solid is given by $vol ( S ) =4\pi r 3 /3$. An illustration of this can be seen on the left of Fig. 1.
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 (x^{2} + y^{2} = R^{2}) 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 freeslip boundary condition u_{3} = 0 at Γ_{top}.
D. Setup
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.
.  Geometry .  Boundary conditions .  

Experiment .  R (m) .  H (m) .  $ r S $ (m) .  h_{0} (m) .  Γ_{wall} ∪ Γ_{bottom} .  Γ_{top} .  $I$ . 
PTFE6  0.055  0.2  0.003  0.161 661 6  u = 0  u_{d} = 0  $u= v S $ 
Rubber22  0.011  0.146 120 3 
.  Geometry .  Boundary conditions .  

Experiment .  R (m) .  H (m) .  $ r S $ (m) .  h_{0} (m) .  Γ_{wall} ∪ Γ_{bottom} .  Γ_{top} .  $I$ . 
PTFE6  0.055  0.2  0.003  0.161 661 6  u = 0  u_{d} = 0  $u= v S $ 
Rubber22  0.011  0.146 120 3 
2. Material parameters
The numerical computation of the abovementioned 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.
.  .  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 × 10^{9}  2.294 52 × 10^{8} 
Rubber22  1361  3.332 89 × 10^{9}  6.667 11 × 10^{5} 
.  .  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 × 10^{9}  2.294 52 × 10^{8} 
Rubber22  1361  3.332 89 × 10^{9}  6.667 11 × 10^{5} 
(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} s^{2}–20 700 000 kg m^{−1} s^{2} (https://eriks.de/content/dam/de/pdf/downloads/dichtungen/oringe/ERIKS_TechnischesHandbuchORinge_de.pdf, p. 20), and the Poisson ratio is ν_{s} = 0.4999.^{44}
E. Quantities of interest
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 h_{0}, i.e., when the ball has traveled $ r S $ vertically. We define t_{*} as the time after release when $dist ( I , \Gamma bottom ) = d S $ relative to t_{0}, i.e., for PTFE6 $ t * = t c S = ( 0 , 0 , 0.009 ) \u2212 t 0 $ and for Rubber22 $ t * = t c S = ( 0 , 0 , 0.033 ) \u2212 t 0 $.

v_{*} the velocity of the ball in the zdirection at t = t_{*} + t_{0}.

f_{*} the vertical component of the force F acting on the ball at time t_{*} + t_{0}.

t_{cont} the time of the first solid contact relative to t_{0}.

t_{jump} 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.

d_{jump} the maximum of $dist ( I , \Gamma bottom ) $ after contact, i.e., the size of the bounce.
An illustration of how these quantities are defined can be seen in Fig. 2.
(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 fluidsurfaces. 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 abovementioned 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_{*}, t_{cont}, d_{cont}, and d_{jump} 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.
Discretization .  Results .  

Method .  [h_{min}, h_{max}] .  Δt .  dof .  nze .  t_{*} .  v_{*} .  f_{*} .  t_{cont} .  t_{jump} .  d_{jump} . 
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 .  [h_{min}, h_{max}] .  Δt .  dof .  nze .  t_{*} .  v_{*} .  f_{*} .  t_{cont} .  t_{jump} .  d_{jump} . 
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 .  [h_{min}, h_{max}] .  Δt .  dof .  nze .  t_{*} .  v_{*} .  f_{*} .  t_{cont} .  t_{jump} .  d_{jump} . 
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 (E_{s} = 5 · 10^{6})  [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 (E_{s} = 2 · 10^{6})  [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 .  [h_{min}, h_{max}] .  Δt .  dof .  nze .  t_{*} .  v_{*} .  f_{*} .  t_{cont} .  t_{jump} .  d_{jump} . 
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 (E_{s} = 5 · 10^{6})  [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 (E_{s} = 2 · 10^{6})  [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 study^{42} 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 x–z 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.
III. REDUCED MODEL
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 Reynoldsnumber 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 twodimensional reduced formulation that is computationally cheaper compared to full threedimensional computations. In Subsection IV B 3, we will also present a fully resolved threedimensional simulation to have a closer look at these rotational effects. To distinguish between the full threedimensional domains and the reduced twodimensional domain, we shall denote objects stemming from the threedimensional setup with a superscript “3d” if there is a potential for ambiguity. For the sake of readability, objects based on the reduced twodimensional setup will not have a special notation.
In order to reduce this threedimensional flow problem into a twodimensional 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^{+} − zplane and transform the weak formulation (9) accordingly. Now, let $\u2207=( \u2202 r \u2202 z )$ and $\Omega = \Omega 3 d \u2229 ( R + \xd7 R ) $ be the reduced twodimensional computational domain. A sketch of the threedimensional domain transformed into two dimensions can be seen on the right of Fig. 1.
IV. NUMERICAL COMPUTATIONS
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 rigidbody 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.
A. Discretizations
We provide the details on the formulations of the different discretization approaches used.
1. Fluid–rigid body interaction in Arbitrary Lagrangian–Eulerian coordinates
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 finiteelement method from Ref. 23 using BDF2 time stepping, together with Taylor–Hood elements in space, which are inf–sup stable in the CutFEM setting^{48} with ghostpenalty stabilization.^{49} In order to obtain the full convergence order of the Taylor–Hood finiteelement pair, we use an isoparametric mapping introduced in Ref. 50 for stationary domains to realize the higherorder geometry approximation.
a. Transformed nitsche terms
b. Transformed ghostpenalty operators
The role of the ghostpenalty operator in the unfitted Eulerian timestepping 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 ghostypenalties provide the necessary implicit extension for the methodoflines approach to the discretization of the timederivative.^{22–24}
c. Contact algorithm
d. Implementation
This discretization is implemented using the finiteelement library Netgen/NGSolve (see Refs. 55 and 56 and ngsolve.org) together with the addon package ngsxfem^{57} for unfitted finiteelement functionality.
The background mesh is constructed by defining a local mesh parameter h_{inner} on the left of the reduced domain where $r<2 d S /3$ and then creating a shape regular mesh with h = h_{max} 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 h_{inner} = h_{max}/4 and in the PTFE6 case as h_{inner} = h_{max}/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 ghostpenalty parameter is γ_{u,e} = 0.1, the cutstability ghostpenalty 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 dist_{0} = 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 dist_{0} = 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 Nitschebased method for the FSI coupling as presented in Ref. 39.
a. Solid bilinear form and FSI coupling
b. Discretization and stabilization
In order to resolve the interface $I$ within the discretization, we use the locally modified finiteelement method introduced in Ref. 58. This fitted finiteelement 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 subtriangles and subquadrilaterals that the interface is resolved in a linear approximation. In this work, we use equalorder locally modified finite elements of first order, in combination with an anisotropic edgeoriented pressure stabilization term s_{p}(p, q); see Ref. 59 for details. We denote the locally modified finiteelement space of first order by $ X h 1 $. The discrete spaces for fluid velocity and pressure V_{h} and Q_{h} and for solid displacement and velocity W_{h} and Z_{h} are given by applying the respective Dirichlet conditions to the locally modified finiteelement space $ X h 1 $ and by restricting the degrees of freedom to the fluid or solid subdomain, respectively.
c. Contact treatment
When a part $ I c \u2254I\u2229 \Gamma 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 nopenetration condition by a small ϵ = ϵ(h) > 0 such that a very thin meshdependent fluid layer remains at all times.
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.
d. Implementation
The described algorithms and equations have been implemented in the finiteelement library Gascoigne3D.^{62} We use a Cartesian finiteelement 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 d_{i}, i = 0, …, m. In the Rubber22 case, we choose, for example, d_{0} = 10^{−2} m and d_{1} = 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.
B. Results
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.
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 nonzero 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 nonzero 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.
If we look at the precontact 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 dist_{0} 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 x–z plane at t = 0.575 in Fig. 5.
We again start by inspecting the precontact 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 precontact 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 threedimensional 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 precontact dynamics are essentially independent of the elasticity parameters—a variation of the elasticity modulus E_{s} changes the rebound height d_{jump} significantly. Here, a softer material (E_{s} = 2 · 10^{6} Pa) leads to a larger rebound as more elastic energy is taken up through the deformation during the impact.
3. Threedimensional computation including rotational effects
One of the challenges in performing the experimental study^{42} was to limit the horizontal deflection of the falling objects, i.e., to keep them close to the centerline. Identifying the source of these threedimensional 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 nonuniform distribution of the mass or the surface roughness but also in experimental inaccuracies, e.g., during the release process.
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 x–y 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.011m$. 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.
V. SUMMARY AND CONCLUSION
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 precontact 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 threedimensional 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 Reynoldsnumber 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.
ACKNOWLEDGMENTS
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).
DATA AVAILABILITY
APPENDIX A: DETAILS ON THE NUMERICAL REALIZATION
1. Fluid–rigid body interaction in arbitrary Lagrangian–Eulerian coordinates
In time, we discretize the Navier–Stokes equations and the rigidbody 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 trick^{66,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
b. Implementation
APPENDIX B: COMPUTATIONAL TEST CASES
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=\u22120.01 ( 1 \u2212 ( x 1 2 + x 2 2 ) / R 2 ) $ on Γ_{top}, noslip u = 0 on $ \Gamma wall \u222aI$, 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 nonconforming, continuous testfunctions $ w \u0302 = ( 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.
Discretization .  Results .  

Method .  [h_{min}, h_{max}] .  dof .  nze .  F_{z} . 
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 .  [h_{min}, h_{max}] .  dof .  nze .  F_{z} . 
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. Nonstationary flow test with prescribed motion
We compute this again using the reduced formulation with the rigidbody discretizations. The quantitative results can be seen in Table VI, while the force functional is shown over time in Fig. 7.
Discretization .  Results .  

Method .  [h_{min}, h_{max}] .  Δt .  dof .  nze .  F_{z,max} .  t_{z,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 .  [h_{min}, h_{max}] .  Δt .  dof .  nze .  F_{z,max} .  t_{z,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 