二、 Poisson Equation with Jump Discontinuities
2.5 Numerical Example
× +
d
d d
d d
ph (1 − p)h qh
(1 − q)h
Figure 2: Interpolation on point ∗ using the four corner points.
2.5 Numerical Example
We consider three modified Poisson equations with exact solutions on Ω = [−1, 1]×[−1, 1]
with jump discontinuities defined on a ellipse Γ : x2
a2 +y2 b2 = 1.
We show two test cases: a = 0.6, b = 0.4, and a = 0.8, b = 0.2.
Example. 1 Example. 2 Example. 3
ψ− 1 x2− y2 sin(x) cos(y)
ψ+ 1 + ln(2px2+ y2) 0 0
ω− 0 0 −2 sin(x) cos(y)
ω+ 0 0 0
Table 1: Three test cases for Eq. (2.1).
N k · k∞ ratio k · kL2 ratio Example. 1 16 3.0347E-02 0.0000 9.1019E-03 0.0000 32 5.4561E-03 5.5621 6.4576E-04 14.0949 64 1.7111E-03 3.1887 2.2393E-04 2.8838 128 5.3475E-04 3.1997 1.6647E-04 1.3452 256 8.7802E-05 6.0904 1.6368E-05 10.1703 Example. 2 16 7.2378E-02 0.0000 1.0009E-02 0.0000 32 1.7418E-02 4.1553 1.3929E-03 7.1855 64 2.7228E-03 6.3972 4.1537E-04 3.3535 128 1.0008E-03 2.7206 1.3772E-04 3.0160 256 1.0402E-04 9.6208 1.6698E-05 8.2479 Example. 3 16 4.2294E-02 0.0000 6.2970E-03 0.0000 32 7.7546E-03 5.4540 7.7777E-04 8.0963 64 1.2982E-03 5.9732 1.0272E-04 7.5721 128 5.6115E-04 2.3135 1.1584E-04 0.8867 256 1.4999E-04 3.7413 2.3262E-05 4.9799 Table 2: Test results for a = 0.6, b = 0.4.
In the first case, the convergent rate fluctuates in a wide range. The accuracy may be affected by how the spline approximates the actual ellipse interface. General speaking, it has second order accuracy on average.
N k · k∞ ratio k · kL2 ratio Example. 1 16 2.1539E+00 0.0000 5.7049E-01 0.0000 32 5.4730E-01 3.9355 9.5492E-02 5.9743 64 1.2800E-02 42.7584 6.9241E-04 137.9125 128 2.5374E-03 5.0446 1.9674E-04 3.5193 256 9.1706E-04 2.7668 4.2801E-05 4.5967 Example. 2 16 3.9062E-01 0.0000 4.0670E-02 0.0000 32 2.3950E-01 1.6310 4.0821E-02 0.9963 64 2.1336E-02 11.2252 1.0262E-03 39.7800 128 4.1392E-03 5.1545 3.9887E-04 2.5727 256 1.4158E-03 2.9237 5.4272E-05 7.3495 Example. 2 16 5.6932E-01 0.0000 7.1246E-02 0.0000 32 1.2497E-01 4.5556 1.3391E-02 5.3204 64 8.8370E-03 14.1420 5.5884E-04 23.9622 128 1.5545E-03 5.6848 9.7281E-05 5.7446 256 5.4722E-04 2.8408 3.2040E-05 3.0363 Table 3: Test results for a = 0.8, b = 0.2.
In the second case, the curvature varies severely. In our test, the cartesian grid cannot fetch the interface behavior as N < 64 in this table. As N exceeds this value, we also observe second order accuracy on average.
Figure 3: The results for the above examples, a = 0.6, b = 0.4.
3 Stokes Flow with Singular Force
3.1 Stokes Equation
The steady Stokes equation is of this form
−∇p + µ∆u = f ,
∇ · u = 0.
where u = u(x), v(x) is the fluid velocity, p = p(x) is the fluid pressure, and the µ is the constant fluid viscosity. The force density f (x) = (f1(x), f2(x)) is applied by the immersed interface Γ into the fluid, which the boundary force is F (s) = (F1(s), F2(s)).
More clearly, the singular force f has support only on Γ and has the form f (x) =
Z
Γ
F (s)δ(x − X(s)) ds
where δ is the Dirac delta function defined in 2-D. The boundary configuration of Γ is given by X(s) = (X(s), Y (s)), 0 ≤ s ≤ L, where the are-length parameter s tracks every material point of the initial configuration.
We write a 2-D vector form of this equation, which yields
−px
One can expect that the velocity and the pressure will not be smooth across the boundary. In fact the boundary force F plays an important role to the jump condition of u, p, and their derivatives. In the computational fluid dynamics, we usually impose the no-slip condition of fluid velocity, i.e., [u] = 0 along the interface. Thus, we have the following theorem.
Theorem. Let u, p be solution of Stokes’ equation. Then the jump across the interface satisfies where n, τ are unit normal and tangent vector respectively on the interface.
We will give a proof in the next section.
3.2 Marker and Cell Method
We use marker and cell (MAC) method and take the computation domain Ω to be square (or rectangular) with staggered grid in Cartesian coordinates. For simplicity, we take the mesh size h = ∆x = ∆y to be equal and let Nx, Ny be the number of mesh in the x-and y- directions respectively. Usually, we pick the same number of grid points in each direction, and denote that N = Nx = Ny.
r
Figure 4: The computational domain Ω using staggered Grid with mesh size h.
3.3 Numerical Scheme
In the following numerical scheme, we assume s be an approximation of arc-length pa-rameter such that |∂X/∂s| is nearly equal to 1. For simplicity, we denote Fn = F · n and Fτ = F · τ . By translating the boundary force F to the jump condition, we have
Consequently, by taking divergent to the above equations we have
∆p = 0,
We need to impose a boundary condition for the above Laplace system. For the test example, we give the Dirchlet boundary condition with given solution p calculated on the boundary. In general, we consider the periodic boundary condition. Although any two solutions of this system may be up to a constant, it is fine since we need only ∇p.
Here we summarize the solution procedures.
1. Solve the Laplace equation
∆hp + C(p) = 0
and obtain p, where the correction tern pi,j)/h respectively. It needs some correction when the p-grid points straddle across the interface. Let us explain the procedure more clearly by taking an example.
Note that the location of ui,j is the middle point between those of pi,j and pi,j. Let pi+1,j ∈ Ω− and pi,j, ui,j ∈ Ω+. The key point is to correct the point on the different side with the middle point. We need to add the correction term if the point is in Ω−, and vice versa. In this case we need to add the correction term Ci+1,j(p) of pi+1,j, which is defined as above. That is, the approximation of px can be expressed by ((pi+1,j + Ci+1,j(p) ) − pi,j)/h. We will discuss the detail in the next section.
3. Compute [∇p] which appears in the correction term C(u) later. Recall that τ = (τ1, τ2), the unit tangent vector, is equal to (−n2, n1) at each point on Γ, where
4. Solve the Poisson equation
µ∆hu + C(u) = ∇hp
3.4 Handling the Derivatives with a Crossing Interface
Consider the derivatives approximated by center difference scheme with a crossing inter-face. In general, there are four cases about locations of grid points and the interinter-face. We use a table to illustrate how to calculate the first order derivative at a certain point.
Case Figure Difference formula
Table 4: The four cases about locations of grid points L, M, R and the interface. h denotes the length of the interval [L, R], and M is the middle point in [L, R] here. CL, CRare the correction terms at grid points L, R respectively.
3.5 Numerical Example
Example (1). Normal forces on a circle. Consider a circular boundary of radius R = 1 parametrized by x(θ) = (cos(θ), sin(θ)) in Ω = [−2, 2] × [−2, 2] and define the force on the boundary by
f (θ) = 2 sin(3θ)x(θ)
so that the forces are normal to the boundary. The exact solution is given by
p(r, θ) =
N kuerrk∞ Ratio kverrk∞ Ratio kperrk∞ Ratio 16 6.31625996E-02 - 5.71511520E-02 - 1.80064548E-01 -32 1.15154249E-02 5.4850 1.04312328E-02 5.4788 6.70713180E-02 2.6846 64 2.46761660E-03 4.6666 3.44416038E-03 3.0286 2.41626231E-02 2.7758 128 2.31427096E-04 10.6626 6.66326129E-04 5.1688 3.79474345E-03 6.3673
Table 5: A grid refinement analysis of u, v, p for Example(1).
Example (2). Tangential forces on a circle. We consider a similar case. The domain Ω and the circular boundary are defined as Example(1). Define the force on the boundary by
f (θ) = 2 sin(3θ)x0(θ)
so that the forces are tangential to the boundary. The exact solution is given by
p(r, θ) =
(−r3cos(3θ), r < 1
−r−3cos(3θ), r > 1 u(r, θ) =
(1
8r2cos(2θ) + 161r4cos(4θ) −14r4cos(2θ), r < 1
−18r−2cos(2θ) + 165r−4cos(4θ) − 14r−2cos(4θ), r ≥ 1 v(r, θ) =
(−18r2sin(2θ) + 161r4sin(4θ) + 14r4sin(2θ), r < 1
1
8r−2sin(2θ) + 165r−4sin(4θ) − 14r−2sin(4θ), r ≥ 1 We note that on the interface R = 1,
[p] = 0, ∂p
∂n
= 6 cos(3θ), [∆p] = 0, ∂2[p]
∂s2 = 0 and
[u] = 0, ∂u
∂n
=
2 sin(3θ) sin(θ)
−2 sin(3θ) cos(θ)
, [∆u] = 3 cos(2θ) + 3 cos(4θ)
−3 sin(2θ) + 3 sin(4θ)
, ∂2[u]
∂s2 = 0.
N kuerrk∞ Ratio kverrk∞ Ratio kperrk∞ Ratio 16 5.56812172E-02 - 1.84446032E-02 - 1.92715009E-01 -32 9.06017465E-03 6.1457 5.85585435E-03 3.1497 2.08500578E-02 9.2429 64 3.34797871E-03 2.7061 1.43026464E-03 4.0942 5.88284456E-03 3.5442 128 6.04375427E-04 5.5395 3.19944420E-04 4.4703 1.32162086E-03 4.4512
Table 6: A grid refinement analysis of u, v, p for Example(2).
Example (3). Combine the normal and tangential forces on a circle. Since the Stokes equation is linear, we combine the above two example by taking summation.
N kuerrk∞ Ratio kverrk∞ Ratio kperrk∞ Ratio 16 1.17188492E-01 - 6.20299174E-02 - 2.02998893E-01 -32 1.22218170E-02 9.5884 1.25574945E-02 4.9396 8.55532629E-02 2.3727 64 3.73073530E-03 3.2759 3.67928831E-03 3.4130 2.61792711E-02 3.2679 128 5.23096557E-04 7.1320 7.33869339E-04 5.0135 4.55700360E-03 5.7448
Table 7: A grid refinement analysis of u, v, p for Example(3).
Example (4). A moving interface problem. This example is given in [8]. We consider a moving interface problem which has uniform density and viscosity. The initial interface is given by
r(θ) = r0(1 + sin(kθ)), 0 ≤ θ ≤ 2π.
The initial velocity and pressure are all set to be zero. The only driven force of the fluid motion is the surface tension which is proportional to the curvature κ with a constant E, that is,
Fn = −Eκ, Fτ = 0.
The velocity is smooth but the pressure has a non-constant jump [p] = Fn, ∂p
∂n
= 0.
We simply take an interface moving formula which says Xn+1= Xn+ Un+1∆t.
In our test, we take r0 = 0.5, k = 5, = 0.2, E = 0.05 and µ = 1.
Figure 5: The numerical evolution of the moving interface using the IIM coupled with the interpolation formula introduced in previous section. Here we pick the number of grid:
Nx = Ny = N = 64, number of interface marker: M = 42, and ∆t = 0.002.
Figure 6: The measure of the bubble area. The geometry of interface eventually converges to a circle with conserving measure of area.
4 Navier-Stokes Equation
4.1 Introduction to Navier-Stokes Equation
The Navier Stokes Equation is of this form ρ ∂u force defined on the computational domain. Note that the divergent-free condition ∇·u = 0 denotes the fluid is incompressible.
Here we have some important theorems for Navier-Stokes equation, which can also be applied to time-independent Stokes equation.
Lemma. If [u] = 0, then
Theorem. Let u, p be solution of Navier-Stokes equation without the external force g.
Then the jump across the interface satisfies [p] = F · n
where n, τ are unit normal and tangent vector respectively on the interface Γ.
Proof. First we choose a banded domain Ω(t) enclosing the interface Γ(t) with the outer and inner boundary Γ+ (t) and Γ− (t). The is a small parameter of the distance from Γ(t) to Γ+ (t) and to Γ−(t), where Γ(t) = {X(s, t), 0 ≤ s ≤ L} be the parametric curve with parameter s and dependence on time t. Let φ(x) be any smooth function with compact support defined on Ω(t). We multiply φ(x) on both sides of Navier-Stokes equation and integrate over Ω(t), which yields
Z
Using the relation of the transport theorem in the first term of the left hand side, and the fact of Ω(t) enclosing Γ(t) and the definition of Dirac delta function in the last term of right hand side, we have
d
Let us consider the first component of the above equation, that is, d
Assume u is smooth enough and bounded, the first term of left hand side tends to zero as → 0. Moreover as → 0. Therefore, since φ is arbitrary, we have
[p] |∂X/∂s| n1 = µ ∂u
∂n
|∂X/∂s| + F1. Similarly, for the second component, we have
[p] |∂X/∂s| n2 = µ ∂u
∂n
|∂X/∂s| + F2.
Multiplying n1, n2 on both sides of the above equations respectively, summing them, and by previous lemma, we obtain the first equality
[p] = F · n
|∂X/∂s|.
Now substituting this pressure jump into the above relations, we get
In order to derive the normal derivative jump of the pressure, let us apply the di-vergence operator to both side of Navier-Stokes equation and use the incompressibility constraint, we have
∇ ·
As before, we multiply φ on both sides and integrate over Ω, which yields Z and, using twice the Green’s identity,
Z
Using the definition of the dipole distribution (the generalized derivation of the Dirac
delta function), we have
Since φ is arbitrary, we obtain the normal derivative jump for the pressure
∂p Thus, the proof is complete.
Remark. The results can be generalized to the case with additional bounded piecewise continuous force g across the interface. The proof is straightforward without much effort needed, and the result of [u] ,∂u
∂n are still valid. However, the normal derivative jump condition for the pressure should be
∂p
We consider the fluid density ρ ≡ 1, and translate the interface force to the jump conditions. We have
∂u
4.2 Numerical Scheme
There are 3 steps to approximate the solution numerically. Here we consider the scheme without interfaces first.
1 In the prediction step, we solve u∗ explicitly where u∗− un
∆t + (u · ∇)un+12 + ∇pn−12 = µ
2∆un+ gn+12, (4.4)
and ˜un+1 by a fast Poisson solver with jump discontinuity,
Note that the advection term (u · ∇)un+12 is approximated by the extrapolation 3
2(u · ∇)un−1
2(u · ∇)un−1.
We explain this step in detail. In the first component of the vector field we have these terms of the space discretization
and the interpolation
vstag = 1
4(vi,j+ vi,j−1+ vi+1,j + vi+1,j−1) from v-grid to u-grid. The second component is similar.
2 We perform the projection step to correction the velocity to be divergent free and update the pressure. According to Helmholtz-Hodge decomposition, we have
˜
un+1 = un+1+ ∆t ∇φn+1, (4.6)
∇ · un+1 = 0.
That is, we want to project the vector field ˜un+1 to a divergent free one, un+1. Taking divergent on both sides, we obtain
∇ · ˜un+1
Note that we should calculate the divergent term ∇ · u∗ on the p-grid, solve the Laplace equation of φ, and then obtain un+1.
3 Finally, the pressure should be corrected. See [6] which says that
∇pn+12 = ∇pn−12 + ∇φn+1− µ∆t
2 ∇(∆φn+1).
Similarly,
pn+12 = pn−12 + φn+1−µ
2∇ · ˜un+1.
It can be obtained by substituting (4.6) into the summation of (4.4) and (4.5).
4.3 Numerical Example
4.3.1 Accuracy Check without Interfaces
Here we introduce some example with exact solutions. In the following examples, we put the time marching step ∆t = h/2 where h is again the mesh size of the staggered grid.
Example (1). Consider the Navier-Stokes equation on a domain Ω = [0, 1] × [0, 1], and give the exact solution
u(x, y, t) = − cos(πx) sin(πy) exp(−2π2t/ Re), v(x, y, t) = sin(πx) cos(πy) exp(−2π2t/ Re),
p(x, y, t) = −0.25 (cos(2πx) + cos(2πy)) exp(−4π2t/ Re).
where Re= 1/µ. We have
∂u
∂x = sin(πx)π sin(πy) exp(−2π2t/ Re), ∂u
∂y = − cos(πx) cos(πy)π exp(−2π2t/ Re),
∂2u
∂x2 = cos(πx)π2sin(πy) exp(−2π2t/ Re), ∂2u
∂y2 = cos(πx)π2sin(πy) exp(−2π2t/ Re),
∂v
∂x = cos(πx) cos(πy)π exp(−2π2t/ Re), ∂v
∂y = − sin(πx)π sin(πy) exp(−2π2t/ Re),
∂2v
∂x2 = − sin(πx)π2cos(πy) exp(−2π2t/ Re), ∂2v
∂y2 = − sin(πx)π2cos(πy) exp(−2π2t/ Re),
∂p
∂x = 1
2sin(2πx)π exp(−4π2t/ Re), ∂p
∂y = 1
2sin(2πy)π exp(−4π2t/ Re).
N kuerrork2 Ratio kverrork2 Ratio kperrork2 Ratio
16 2.5618E-05 - 1.7061E-05 - 6.1794E-03
-32 2.7016E-06 9.4825 1.5956E-06 10.6925 1.2137E-03 5.0913 64 2.8951E-07 9.3316 1.7090E-07 9.3364 2.6025E-04 4.6635 128 3.5376E-08 8.1837 2.0645E-08 8.2780 6.0763E-05 4.2830
Table 8: A grid refinement analysis in L2-Norm for Example(1).
N kuerrork∞ Ratio kverrork∞ Ratio kperrork∞ Ratio
16 8.2522E-05 - 6.6571E-05 - 1.0643E-02
-32 1.5195E-05 5.4308 9.4561E-06 7.0400 2.1667E-03 4.9120 64 2.4185E-06 6.2828 1.4674E-06 6.4441 4.8205E-04 4.4947 128 4.2419E-07 5.7014 2.4059E-07 6.0991 1.1171E-04 4.3151
Table 9: A grid refinement analysis in L∞-Norm for Example(1).
Example (2). Consider the Navier-Stokes equation with µ = 1 on a domain Ω = [0, 1] ×
[0, 1], and give the exact solution u(x, y, t) = 1
4ye−t, v(x, y, t) = −1
4x(1 − x2)e−t, p(x, y, t) = (x2+ y2− 1)2e−t, g1(x, y, t) = −1
4ye−t− 1
16x(1 − x2)e−2t+ 4(x2+ y2− 1)e−tx, g2(x, y, t) = 1
4x(1 − x2)e−t+ 1
4ye−t(−1
4(1 − x2)e−t +1
2x2e−t) − 3
2xe−t+ 4(x2+ y2− 1)e−ty.
Here we have
∂u
∂x = 0, ∂u
∂y = 1 4e−t,
∂2u
∂x2 = 0, ∂2u
∂y2 = 0,
∂v
∂x = 1
4e−t(−1 + 3x2), ∂v
∂y = 0,
∂2v
∂x2 = 3
2xe−t, ∂2v
∂y2 = 0,
∂p
∂x = 4(x2+ y2− 1)e−tx, ∂p
∂y = 4(x2+ y2− 1)e−ty.
N kuerrork2 Ratio kverrork2 Ratio kperrork2 Ratio
16 5.3838E-04 - 8.1989E-04 - 1.3604E-02
-32 1.4644E-04 3.6764 2.2699E-04 3.6120 4.3058E-03 3.1594 64 3.8076E-05 3.8459 5.9219E-05 3.8330 1.4074E-03 3.0594 128 9.6009E-06 3.9658 1.4984E-05 3.9521 3.8226E-04 3.6817
Table 10: A grid refinement analysis in L2-Norm for Example(2).
N kuerrork∞ Ratio kverrork∞ Ratio kperrork∞ Ratio
16 1.3922E-03 - 3.2717E-03 - 8.6817E-02
-32 4.0255E-04 3.4584 9.6117E-04 3.4038 4.8699E-02 1.7827 64 1.1010E-04 3.6562 2.5788E-04 3.7271 2.6233E-02 1.8564 128 2.8170E-05 3.9084 6.6552E-05 3.8748 1.3360E-02 1.9635
Table 11: A grid refinement analysis in L∞-Norm for Example(2).
4.3.2 Cavity Flow
Example (3). Let Ω = [0, 1]×[0, 1] be computational domain. Assume that flow is driven by the upper wall. That is, we have the boundary condition
u = 0 on ∂Ω \ ∂Ωlid
and
u = (u, v) = (1, 0) on Ωlid,
where Ωlid = {(x, y) ∈ ∂Ω, y = 1}. Let N = 128 be the number of partition. We compare our result with [14].
Figure 7: Figures of velocities u and v on x = 0.5 and y = 0.5 respectively. Re = 1000.
Figure 8: A contour plot of a driven cavity with Re= 1000.
4.3.3 Circular Flow with a Line Force
Example (4). Let us assume the curve Γ to be a circle x2+y2 = 14 within a computational domain Ω = [−1, 1] × [−1, 1]. That is, the interface can be represented by
X(s), Y (s) = 0.5 cos θ, 0.5 sin θ
where θ = 2s, 0 ≤ s ≤ π. Let the pressure and velocity components be
u(x, y, t) =
(h(t) yr − 2y
r > 12
0 r ≤ 12
v(x, y, t) =
(h(t) −xr + 2x
r > 12
0 r ≤ 12
p(x, y, t) =
(1 r > 12 0 r ≤ 12
where r = px2+ y2. It is easy to verify the velocity satisfies the incompressibility constraint, and it is continuous but has finite jump across the interface. That is,
∂u
∂n
= −2h(t) sin θ, ∂v
∂n
= 2h(t) cos θ, [p] = 1 and then
Fn = 1, Fτ = −2h(t).
The external force g can be computed to satisfy the Navier-Stokes equation, and yields
g1(x, y, t) =
(h(t)2cos θ 4 − 4r − 1r + h(t) sin θrµ2 + h0(t) sin θ(1 − 2r) r > 12
0 r ≤ 12
g2(x, y, t) =
(h(t)2sin θ 4 − 4r − 1r − h(t) cos θrµ2 − h0(t) cos θ(1 − 2r) r > 12
0 r ≤ 12 .
One can verify that on Γ, (r = 1/2), [g] · n = h(t)2
4 − 4r − 1 r
− 0 = 0.
Some results are given here. In our tests, we assume the boundary force is given. We observe that the accuracy in u, v is nearly second order.
N kuk∞ Ratio kvk∞ Ratio kpk∞ Ratio 16 6.2133E-03 0.0000 6.0783E-03 0.0000 1.3309E-01 0.0000 32 2.0049E-03 3.0990 1.9912E-03 3.0526 6.5489E-02 2.0322 64 5.6914E-04 3.5228 5.6695E-04 3.5121 4.4361E-02 1.4762 128 1.4834E-04 3.8367 1.4805E-04 3.8295 1.5033E-02 2.9509 N kukL2 Ratio kvkL2 Ratio kpkL2 Ratio 16 1.4411E-03 0.0000 1.4221E-03 0.0000 2.5676E-02 0.0000 32 3.8649E-04 3.7287 3.8360E-04 3.7073 1.0012E-02 2.5645 64 1.1245E-04 3.4370 1.1209E-04 3.4223 3.3265E-03 3.0099 128 3.2281E-05 3.4835 3.2235E-05 3.4772 1.2290E-03 2.7066
Table 12: The grid refinement analysis for u, v, and p. h(t) = 1.
N kuk∞ Ratio kvk∞ Ratio kpk∞ Ratio
16 1.8770E-03 0.0000 1.8741E-03 0.0000 3.7830E-02 0.0000 32 6.3666E-04 2.9483 6.4121E-04 2.9227 2.0364E-02 1.8577 64 1.8387E-04 3.4625 1.8356E-04 3.4932 1.4330E-02 1.4211 128 4.8301E-05 3.8068 4.8899E-05 3.7539 4.7389E-03 3.0238 N kukL2 ratio kvkL2 ratio kpkL2 ratio 16 4.3742E-04 0.0000 4.3582E-04 0.0000 7.8480E-03 0.0000 32 1.2249E-04 3.5712 1.2222E-04 3.5660 3.2169E-03 2.4396 64 3.5968E-05 3.4054 3.5932E-05 3.4013 1.0809E-03 2.9762 128 1.0295E-05 3.4937 1.0291E-05 3.4917 4.1351E-04 2.6139 Table 13: The grid refinement analysis for u, v, and p. h(t) = 1 − e−t.
5 Concluding Remarks
We have introduced a numerical technique for the computation of interface problems.
The result shows the behavior near the immersed interface can be captured well, and in general second order accuracy can be achieved. Moreover, the method is efficient since fast solvers can still be applied.
One of the future works is to solve piecewise constant coefficient Poisson equation
∇ · (β±∇ψ±) = ω± on Ω±.
The immersed interface method can also handle this type of equation. Thus, we can extend this solver to solve piecewise constant coefficient Stokes equation or Navier-Stokes equation.
It also need to be considered that how the spline we constructed approximates the actual interface, and how accurate the curvature and the distance to the foot of perpen-dicular we obtain. The error generated by approximation of an interface by an spline may lead to the fluctuation of convergent ratio in grid refinement analysis.
It is a challenge that extending this method to three-dimension. Although there are some fast Poisson solvers in 3-D, constructing the interface of surface and finding foot of perpendicular seem to be a heavy task. It is important that maintaining the efficiency and accuracy when extending this method into 3-D.
References
[1] Alfio Quarteroni, Riccardo Sacco, Fausto Saleri, “Numerical Mathematics.” Texts in Applied Mathematics 37, Springer.
[2] Alexandre J. Chorin, Jerrold E. Marsden, “A Mathematical Introduction to Fluid Mechanics, 3rd. ed.” Texts in Applied Mathematics 4, Springer-Verlag.
[3] J. Stoer, R. Bulirsch: “Introduction to Numerical Analysis.” Springer-Verlag, 1983 [4] Bu-Qing Su, Ding-Yuan Liu: “Computational Geometry: Curve and Surfce
Model-ing.” Academic Press, INC.
[5] M. C. Lai, “Jump conditions and test examples for the 2D Navier-Stokes equations involving an immersed boundary.”
[6] R. Cortez, “The Method of Regularized Stokeslets.” SIAM Journal of Scientific Com-puting, vol. 23, No. 4, pp. 1204-1225.
[7] J. Thomas Beale and Anita T. Layton, “On The Accuracy of Finite Difference Meth-ods for Elliptic Problems with Interfaces.” Aug 19, 2005.
[8] Zhilin Li and M. C. Lai, “The Immersed Interface Method for the Navier-Stokes Equations with Singular Forces.” Journal of Computational Physics, Vol. 171, n. 2, pp. 822-842.
[9] A. McKenny and L. Greengard and A. Mayo, “A Fast Poisson Solver for Complex Geometries.” Journal of Computational Physics, Vol. 118, pp. 348-355.
[10] M. C. Lai and C. S. Peskin, “An immersed boundary method with formal second-order accuracy and reduced numerical viscosity. ” Journal of Computational Physics, Vol. 160, pp. 705-719. (2000)
[11] Boyce E. Griffith, Charles S. Peskin, “On the order of accuracy of the immersed boundary method: Higher order convergence rates for sufficiently smooth problems.”
Journal of Computational Physics, Vol. 208, pp. 75-105.
[12] D. Russell. Z. Jane Wang, “A cartesian grid method for modeling multiple mov-ing objects in 2D incompressible viscous flow.” Journal of Computational Physics Vol. 191 (2003), pp.177-205.
[13] Sheng Xu, Z. Jane Wang, “An immersed interface method for simulating the in-teraction of a fluid with moving boundaries.” Journal of Computational Physics, in press.
[14] U. Ghia, K. N. Ghia, and C. T. Shin, “High-Re Solutions for Incompressible Flow Us-ing the Navier-Stokes Equations and a Multigrid Method.” Journal of Computational Physics Vol. 48, pp.387-411.