4.3 Incompressible Navier Stokes
5.1.1 Interface Thickness
Γ (x) is not unique, meaning that the point xis equidistant to at least two points on Γ. An example is the center of the circle in Figure 5.2,x= (0,0), which is equidistant to the entire interface.
The unit normal of the interface is determined by nΓ= ∇φ
|∇φ|
φ=0 (5.4)
and we can theoretically find the closest point on the interface from xmin
Γ (x) =x−φ(x)w(x), w(x) = ∇φ(x)
|∇φ(x)| (5.5)
since the shortest distance fromxto the interface is the straight line perpendicular to the interface at xmin
Γ (x), and the gradient at x will take the direction of the interface normal at xmin
Γ (x), a side effect of φ being a signed distance. This is again not possible everywhere, i.e. whenxis equidistant to at least two points on Γ and∇φ(x) is not defined.
The mean curvature of the interface can be defined as κ=∇ · ∇φ Since for a distance function|∇φ|= 1, above could be simplified. However we may in certain cases allow the level set function to only be approximately a signed distance function, hence e.g. Equation (5.5) may not be exact, and it is appropriate to keep the normalization constant in the presentation here.
5.1.1 Interface Thickness
To avoid handling discontinuities in e.g. the viscosity µ(φ) and densityρ(φ), the interface is given a thickness of size, i.e., we use a smoothed Heaviside function
H(φ, ) =
The smoothed Heaviside function is depicted in Figure 5.3. This is only one of many possible definitions of a smoothed Heaviside function [69]. Using a high order method to enable high precision, it is important that the function is as smooth as possible. The above Heaviside function has continuous first and second derivative, while the third is discontinuous at φ= ±. It is possible to design higher order polynomials or combinations of polynomials and trigonometric functions which are as smooth as required.
The interface thicknessshould be chosen as small as possible for accuracy but large enough to stabilize the system: If chosen too small, the Heaviside function
x
Figure 5.3: Smoothed Heaviside function
will be too steep, hence it becomes difficult to represent and oscillations due to Runge phenomenon may occur, making the system stiff and difficult to solve. A typical choice is of the size of an element. Note that does not have to be a global constant, but may vary throughout the domain.
Example 5.1: Smoothed Heaviside and δ-functions
For the Heaviside function, determine the polynomialp(x) of degreeN= 2k+ 1 p(x) =
Inserting (5.8) we get 2k+ 2 equations and equally many unknowns. Then define
H(φ, ) =
which will havekcontinuous derivatives everywhere. Figure 5.4 shows 5 different smoothed Heaviside functions of different polynomial degrees. The polynomial with highest degree is the most smooth but has also the steepest gradient atφ= 0. Hence the coefficients and the size of the derivatives grow larger with the degree of smoothness, which also makes the differentiation less accurate.
When a fully smooth function is needed, and the thickness is less important, Hs(φ, ) = 1
2+ φ
2p
φ2+2 (5.11)
can be used. Figure 5.5 shows that this function actually never reaches zero and one, but
5.1.1. Interface Thickness 79
−1 0 1
0 0.5
1 P1
P3 P5 P7 P9
Figure 5.4: Smooth Heaviside function,= 1, for different polynomial degrees
−2 −1 0 1 2
0 0.5 1
Figure 5.5: Very smooth Heaviside function,= 1
only approaches asymptotically. It is, however, infinitely smooth. Theδ-function can be defined as:
δ(φ) = ∂
∂φH(φ). (5.12)
Hence smoothed versions of theδ-function can be defined from the Heaviside function.
Figure 5.6 shows similar plot for the smoothedδ-functions, where theδ-function denoted
−2 −1.5 −1 −0.5 0 0.5 1 1.5 2
0 0.5 1 1.5
P0 P2 P4 P6 P8 s
Figure 5.6: Smooth delta functions,= 1
Piis based on a polynomial of degree i, andsdenotes the very smooth version.
Note that all the smoothed δ-functions have the property:
Z ∞
−∞
δ(φ, ) dx= 1. (5.13)
Most of them actually have the property that Z
−
δ(φ, ) dx= 1, (5.14)
since most of the smoothedδ-function are zero outside±.
End of Example 5.1.
The smoothed δ-function must give the same contribution as the exact δ-function, for example when used for applying the surface tension force, i.e.
Z
We will therefore require that our smoothedδ-function fulfill l(Γ) =
wherel(Γ) is the length of the interface.
Example 5.2: Error in smoothedδ-functions
Consider the three different interfaces in Figure 5.7. We will test the accuracy of
integra-Γ
tion ofδ(φ, ) in Equation (5.16) for different values ofand the degree of the polynomial used when creating the δ-function, Pi, from Figure 5.6. The test is carried out in a [−12;12]×[−12;12] domain. The circle has radius 0.3.
\Pi 0 2 4 6 8
0.1 1.7e-03 8.4e-04 1.5e-04 5.8e-05 1.7e-04 0.2 1.9e-03 9.2e-05 3.1e-06 6.2e-07 3.3e-06 0.3 4.1e-03 4.7e-06 5.0e-06 3.3e-07 2.9e-07 0.4 3.7e-03 2.1e-05 4.4e-06 1.1e-07 8.7e-08
Table 5.1: Error in smoothedδ-function for linear interface, Figure 5.7(a)
5.1.1. Interface Thickness 81
\Pi 0 2 4 6 8
0.05 1.4e-02 3.0e-03 2.5e-03 2.5e-04 1.4e-03 0.10 1.6e-03 5.8e-04 2.6e-04 5.2e-05 6.4e-04 0.15 5.1e-03 4.8e-05 7.1e-05 5.2e-05 5.0e-05 0.20 8.5e-03 1.1e-04 2.1e-05 4.5e-06 5.0e-06
Table 5.2: Error in smoothedδ-function for circular interface Figure 5.7(b)
\Pi 0 2 4 6 8
0.1 7.6e-02 7.5e-02 6.3e-02 5.5e-02 4.9e-02 0.2 2.0e-01 1.5e-01 1.3e-01 1.1e-01 9.8e-02 0.3 3.0e-01 2.3e-01 1.9e-01 1.6e-01 1.5e-01 0.4 4.0e-01 3.0e-01 2.5e-01 2.2e-01 2.0e-01
Table 5.3: Error in smoothedδ-function for linear interface Figure 5.7(c) Results in Table 5.1, 5.2, and 5.1 show the error when integrating theδ(φ, )-function in Equation (5.16), using a spatial discretization having 226 elements of polynomial order 5. Table 5.1 and 5.2 shows that the error decay for increasing polynomial degreePi, and to some degree also with the thickness of the interface. However, in Table 5.3 the error is smallest for smallest, more or less independent of the polynomial degreePi. This is
Γ
(a)
∂Ω
Γ
(c)
∂Ω
Figure 5.8: Integrating close to the boundary
due to the fact that integrating up to the boundary will be inaccurate. Figure 5.8 shows the problem, where the gray part is where the smoothedδ-function should be integrated, however the gray-white part outside the domain are integrated. In case (a) the white boxed area is integrated and cancels out the gray-white part outside the domain, while in case (c) it does not. The error in case (c) will be smaller for smaller , since the part outside the domain decreases with.
The comforting news is to see that for the circular interface, case (b), the smoothed δ-function works well.
End of Example 5.2.
It is important for the results in Example 5.2 that the level set functionφis a true distance function. If φis not a distance function, the interface thickness will change and the property (5.13) will not be true: Assumeφis a distance function, thenψ(x) =aφ(x),a6= 1, is not a distance function,∇ψ=a∇φ=a. Figure 5.9
−2 −1 0 1 2 0
0.5 1
ψ
δ(ψ,1)
ψ=φ ψ=0.75φ
Figure 5.9: δ-function when ψis not a distance function
shows the smoothed δ-function for a = 1 and a = 34, and the area under the smoothed δ-function is bigger fora= 34 than fora= 1. This means that one of the consequences of the level set not being a distance function is that we do not apply the correct surface tension force.