III. Applications 137
11.2. Continuum mechanical model
In this case, the stationarity condition ˜E0=0 implies H=− σw
ρgR=σcosθ
ρgR . (11.9)
Non-dimensional formulation: It is convenient to rewrite eq. (11.7) to arrive at a dimensionless form of the functional, i.e.
E˜=πR2σ 2H R
σw
σ +ρgR2 2σ
H R
2
+1
!
=πR2σ
2 ˜Hσ˜w+Bo ˜H2 2 +1
.
(11.10)
Here ˜H=H/Ris the dimensionless rise height, ˜σw=σw/σ is the dimensionless wetting energy satisfying the dimensionless form of the Young–Dupr´e equation
cosθ+σ˜w=0. (11.11)
Hence, the dimensionless form of equation (11.1) reads as H˜ =−2 ˜σw
Bo =2 cosθ
Bo . (11.12)
In the two-dimensional case, the simplified functional reads as
E˜=Rσ 2 ˜Hσ˜w+Bo ˜H2+1
(11.13) leading to the equilibrium rise height
H˜ =−σ˜w
Bo =cosθ
Bo . (11.14)
Improved stationary rise height estimate: It turns out that the approximations made in the derivation of (11.1) and (11.9) lead to systematic error which is clearly visible in the numerical simulations and which should also be measurable in an experiment. The main simplification is that the details of the interface shape are ignored approx-imatingΣas a flat interface in the simplified energy functionals (11.7) and (11.8). It is known that for small Bond numbers, the interface can be well approximated by a spherical section [Con68]. Hence, an improved description can be achieved by introducing an additional degree of freedom and minimizing the energy over liquid columns with a spherical interface. Note that this is also an approximation since the true liquid-gas interface must have a non-constant curvature (see [GFB19]).
The problem can be simplified further by assuming that the interface is a spherical cap with a contact angle de-termined by the Young-Dupr´e equation. Then the additional mass in the interface region can be computed and a correction to the rise height is obtained. The correction derived in [Gr¨u20b] for the two-dimensional case reads as
H=σcosθ
ρgR −∆H with ∆H= R 2 cosθ
2−sinθ−arcsin cosθ cosθ
>0. (11.15) See also [LLL18], where a similar correction term is found.
Chapter 11. The capillary rise problem
Figure 11.2.: Computational domain for the continuum mechanical model.
11.2. Continuum mechanical model
11.2.1. Governing equations and boundary conditions
We consider the standard model (3.22) with a constant slip length and a fixed contact angle, i.e.
ρDv
Dt −η∆v+∇p=0, ∇·v=0 inΩ\Σ(t), JvK=0, Jp1−2ηDKnΣ=σ κnΣ onΣ(t),
hv,n∂Ωi=0,vk+2L(Dn∂Ω)k=0 on∂Ωs\Γ(t), VΣ=hv,nΣi onΣ(t), VΓ=hv,nΓi, hnΣ,n∂Ωi=−cosθ onΓ(t),
(11.16)
to model the dynamic capillary rise problem. In order to simplify the computational setup, we only solve for the flow inside the capillary and prescribe effective boundary conditions at the (artificial) in- and outflow boundaries
∂Ωband∂Ωt (see Fig. 11.2). Neglecting the hydrostatic pressure in the gas phase, we force the pressure to be equal at the top and bottom boundaries, i.e.
p=0 at∂Ωt∪∂Ωb. (11.17)
In order to allow the liquid to flow across the top and bottom boundaries, we set3
h(∇v)n∂Ω,n∂Ωi=0 at∂Ωt∪∂Ωb, (11.18) wheren∂Ωis the unit outer normal field to
∂Ω=∂Ωs∪∂Ωt∪∂Ωb.
3In general, the definition of appropriate inflow and outflow boundary conditions at artificial boundaries is a non-trivial subject [SG94,HRT96].
However, in the present case, the choice of the velocity boundary conditions at the artificial boundaries showed only little influence on the results since the flow profile is essentially a Poiseuille Flow (see, e.g., [Gr¨u20b, p.70]).
142
11.2. Continuum mechanical model
Moreover, to specify the type of fluid entering at the top and bottom boundaries, we require
∂ χ
∂n∂Ω =0 at∂Ωt∪∂Ωb. (11.19)
11.2.2. Choice of physical parameters
Motivated by the work by Fries and Dryer [FD09], who carried out a dimensional analysis of the ODE model [Bos23], we choose a set of physical parameters (see Table 11.1) to vary the non-dimensional group
Ω= s
9σcosθ η2 ρ3g2R5 =3√
cosθ Oh
Bo, (11.20)
where Bo and Oh denote the Bond and Ohnesorge number defined as Bo=ρgR2
σ and Oh= η
√σ ρR,
in the range ofΩ=0.1,0.5,1,10,100. Owing to the increasing influence of viscosity, varyingΩleads to different rise regimes ranging from a highly oscillatory (smallΩ) to a purely monotone rise to the stationary rise height (largeΩ). Note that
Ωc=2 (11.21)
has been shown to be the critical damping for the classical ODE model [Bos23] (see [QRO99,Gr¨u20b,GSA+20a]).
It is important to note that the parameterΩ, which has been introduced for the three-dimensional case by Qu´er´e et. al [QRO99], does not take into account the slip length. In fact, the slip length shows a strong impact on the rise dynamics in the full continuum mechanical model. We show below that the slip length does also affect the critical valueΩcfor the occurrence of oscillations in the dynamic rise.
In order to simplify the setup for the parameter study, we fix the contact angle to beθ=30◦for all cases. More-over, the radius and the height of the capillary are fixed and the parameters are chosen such that the Jurins height is fixed as four times the radius, i.e.
Rcosθ
Bo =4R ⇔ Bo=cos(π/6)
4 ≈0.217.
We choose an initial rise heighth0=2Rfor all cases. The height of the complete computational domain is chosen to be 8Rwhich is sufficient to capture any rise height oscillations. The density and viscosity in the gas phase is chosen small enough, i.e.
ρ ρg
=1000 and η ηg
=1000,
to make sure that the influence of the gas phase is negligible. This is necessary since we aim at a direct comparison with a free surface formulation of the problem, see below. The physical parameters used for the study of the dynamics in Section 11.3 and [GSA+20a] are summarized in Table 11.1. The reported Camax is the maximum observed value for the capillary number. Note that the test cases serve the purpose of anumericalbenchmark and do not necessarily correspond to the physical parameters of some real fluid-solid system in the lab. In particular, we chose the slip lengthmacroscopically large such that it can be easily resolved by the computational mesh.
Even though this choice is physically unrealistic, it still allows to establish a numerical reference or “numerical benchmark” which contains the wetting as a driving force.
Chapter 11. The capillary rise problem
Table 11.1.: Physical parameters for theΩ-study.
Ω R ρ η g σ θe Camax Eo
- m kg m−3 Pa s m s−2 N m−1 ◦ -
-0.1 0.005 1663.8 0.01 1.04 0.2 30 0.003 0.217
0.5 0.005 133.0 0.01 6.51 0.1 30 0.015 0.217
1 0.005 83.1 0.01 4.17 0.04 30 0.029 0.217
10 0.005 3.3255 0.01 26.042 0.01 30 0.106 0.217 100 0.005 0.33255 0.01 26.042 0.001 30 0.110 0.217
Viscous vs. capillary timestep limit: To estimate the computational costs, we consider the numerical timestep limits already discussed in Section 7.1. The timestep limit due to surface tension is given by(∆t)σ=p
ρ(∆x)3/(4π σ).
Here∆xdenotes the mesh size which is assumed to be constant for the purpose of estimating the computational costs. Since the viscous terms are explicitly discretized in FS3D, there is also a viscous timestep limit for the latter method given by(∆t)η =ρ(∆x)2/(6η). In order to identify which one is the limiting value for a given set of parameters, we take the ratio of the two quantities yielding
(∆t)σ (∆t)η = 3
√π
√σ ρRη R
∆x 1/2
= 3
√πOhN
1 2
cells, (11.22)
whereNcells=R/∆xis the number of computational cells per radius. Hence, the two timestep limits are equal for Ncells∗ = π
9 Oh2=πcosθ Ω2Eo2.
The capillary timestep limit(∆x)σ is dominant forNcells≤Ncells∗ , while the viscous timestep limit(∆x)η is dom-inant forNcells ≥Ncells∗ . Hence the viscous timestep limit dominates over the capillary timestep limit for large Ω.
Estimated number of timesteps: The necessary number of timesteps due to surface tension and viscosity is given by
Nstepsσ = T (∆t)σ =T
s4π σ
ρlR3Ncells3/2 and Nstepsη = T
(∆t)η =T 6ηl ρlR2Ncells2 ,
respectively. HereT denotes the physical time to be simulated. The values forT are roughly estimated from Figure 11.6 in Section 11.3 to obtain an estimate for the number of timesteps, see Table 11.2. Note that the number of timesteps for FS3D is limited by the viscous timestep limit forΩ=10 andΩ=100. In particular, the latter case requires more than 108timesteps if 32 cells per radius are used. Since this would lead to impractical computational costs, we used a coarser grid with 12 cells per radius forΩ=100 (reducing the number of timesteps by a factor of 7). Note the viscous timestep limit can be avoided by an implicit discretization of the viscous terms. The latter approach is not available in the current version of FS3D but will be subject of future work.
Table 11.2.: An estimate of the number of timesteps forNcells=32 and the physical parameters defined in Ta-ble 11.1.
Ω Ncells∗ T [s] Nstepsσ Nstepsη
0.1 5804 10 1.99·105 1.48·104 0.5 232 1 4.98·104 1.85·104 1 58 0.5 1.99·104 1.48·104 10 0.58 2 1.99·105 1.48·106 100 0.006 20 1.99·106 1.48·108
144