二、 Fractional Step Method
2.4 Boundary conditions
First, we discuss the case of Dirichlet boundary conditions where the domain Ω 0, X 0, Y . That is, the boundary conditions 0, , , X, , , are not defined at the boundary. We approximate these boundary conditions as follows:
After we know the discrete boundary conditions, the discrete operators can be written as matrices. Consider the laplacian term , it contains the laplacian operators to the horizontal velocity and vertical velocity, called and .
0 0
The discrete Δ , and Δ , have been introduced already. Since we suppose ∆ ∆ , and can be written as :
0 let’s see a simple example. Consider the two cells case. There are only four horizontal velocity terms and three vertical velocity terms, called , , , ,
, , . See fig. 2.
Fig. 2. The two cells case
Consider the discrete divergence term, we know that
∆ ∆ ∆ ∆ 0
∆ ∆ ∆ ∆ 0
If we rewrite it to the matrix form, then
11
12 simple factorization .
Similarly, we also can get the matrix form of and factorize . In the two
The relation that T still exists in general situation.
The boundary condition term comes from . can be written as:
, , 0, ,0, , , , , 0, ,0, , , , , , 0, ,0, , T
, , , , , , , 0, ,0 , , , , , , , T If there are not Dirichlet boundary conditions at all. That is, some boundary conditions maybe are Neumann boundary conditions. Let the boundary
∂Ω ∂Ω ∂Ω ∂Ω ∂Ω , and the boundary conditions , , | Ω , , , | Ω , ∂ , , | Ω , ∂ , , | Ω are known, where is the outer normal vector of the boundary.
For example, we suppose that ∂Ω , Ω| 0 , ∂Ω ∂Ω\ ∂Ω ,
∂Ω , Ω| X , and ∂Ω ∂Ω\ ∂Ω . In figure 3, it shows the relationship between the boundary condition type and the place. The boundary conditions that we know are 0, , , X, , , , 0, , , Y, ,
0, , , X, , , , 0, and , Y, .
But the thing we really care about is the discrete boundary condition. We have known how to handle the discrete boundary condition when the boundary condition is Dirichlet. Now we discuss the discrete boundary condition when the boundary condition is Neumann.
Fig. 3. The left figure represents the boundary condition type of u. The right figure represents the boundary condition type of v.
13
14
We can use the central difference to approximate , 0, .
1 , 0, 1 ∆ , ,
∆ .
The discrete boundary condition , can be approximated by
, , ∆ . · 1 , 0, 1 ∆
But we do not know , , so this part can be replace by , , or use the velocity at the time step n and n-1 to approximate , . That is,
, 2 , ,
This approximation also can be applied when the boundary condition is Dirichlet. So the discrete boundary condition , can be approximated by
, 2 , , ∆ . · 1 , 0, 1 ∆
Similarly, the other discrete boundary conditions can be written as
, 2 , , ∆ . · 1 , Y, 1 ∆
, 2 , , ∆ . · X, 1 , 1 ∆
, 2 , , ∆ · X, 1.5 , 1 ∆
Notice that we are not use the central difference to approximate , . Because these boundary conditions are approximated by and . After solving , we need to update the discrete boundary conditions. At the next time step, they will be used to compute N .
Although we use the Neumann boundary conditions approximate the boundary velocity, we would rather use the Neumann boundary conditions than use the approximated boundary velocity.
Consider the laplacian operator . For example, we would rather use 0.5 , 0, 1 ∆ than , , ⁄∆ when we compute the approximation of ∆ , . Suppose ∆ ∆ , the discretization of ∆ , is
15
∆ , , 2 , , , , 0.5 , 0, 1 ∆
, 3 , , , 0.5 , 0, 1 ∆
Follow this idea, the matrix and the boundary condition also have some change. In this case, the matrix changes to
0
The corresponding boundary condition also has some change. First, e denote some values
16
The corresponding boundary condition can be written as:
, T/2
where
, , 0, ,0, , , , 0, ,0, , , , , 0, ,0, T h , , , , 0, ,0 , , , , T
, , 0, ,0, , , , 0, ,0, , , , , 0, ,0, T
, , , , , , , 0, ,0 , , , , , , , T Neumann boundary condition also can be used to compute , but we still use the discrete boundary condition and original . If we use the Neumann boundary condition, the matrix would be changed, the symmetrization that T would be destroyed.
3 Matrix Factorization Method
3.1 Mathematical formulation
The immersed boundary method was first introduced by Peskin in [7, 8].
Consider the two-dimensional incompressible Navier-Stokes equation with an immersed massless boundary. Let the two-dimensional domain Ω 0, X
0, Y . Let the immersed boundary ∂ , which is a closed curve, as shown in Figure 4. The immersed boundary is tracked by the parametric function s , 0 s , where is the length of the immersed boundary . The mathematical formulation is
∂ 1
17
∆ , δ d
0
, , δ d
Ω
where , , , is the Lagrangian boundary force defined on the immersed boundary, and , is the velocity at the immersed boundary
. δ is the two-dimensional Dirac delta function.
Fig. 4. The domain Ω , immersed object , immersed boundary , and the Lagrangian points
18
3.2 Discretization
Similar to the discretization of the incompressible Navier-Stokes equations, we can discretize the mathematical formulation. Refer to [9] , the discretization can be written as:
where and are the discrete regularization operator and the discrete interpolation operator, respectively.
The discrete boundary force f and the discrete boundary velocity are defined on the Lagrangian markers . We can regard as , where is a monotone sequence between 0 and , 1 , as shown in figure 4. Let
, , we define the distance between each Lagrangian markers ∆ as follows:
∆ | | for 1 1
∆
If the arc length of each part is known, we would rather define the arc length
∆ then use distance. Then, we use ∆ define ∆ as follows:
∆ ∆ ∆ ⁄ 2
∆ ∆ ∆ ⁄ for 22
Before we introduce and , the discrete delta function should be introduced first. Here, we use the discrete delta function which was introduced by Roma et al [10].
δ , ·
where the function is
5 3| |
To make the scheme simple,
suppose that ∆ ∆
when near the Lagrangian points.
Consider the two cells case, as shown in fig. 5. Let
, , where
and are the horizontal and vertical direction of discrete boun- dary velocity fields, respectively.
Fig. 5. The two cells case with the Lagrangian points
20
Similarly, consider the discrete regularization operator . Let , , where and are the horizontal and vertical direction discrete boundary force, respectively. We can discretize the term , δ d as of the definition of the discrete delta function if we let ∆ .
21
3.3 Matrix factorization method
Put the unknown terms , , and on the left side of the equal sign, and the other terms on the right side. The discretization can be written as:
2 2
The matrix factorization method is a projection approach of the immersed boundary method introduced by Taira and Colonius [9]. Similar to the fractional step method, we have known the facts that and . The discretization can be rewritten to
1 ̂ 1
1 0 2
Let 1, , , ̂ , and
. It can be represented by the matrix form as follow:
1
01
1
0 2
We know that T, but and do not have the symmetric relation directly. So a transformed forcing function was introduced by Taira and Colonius [9]. It satisfies T . That is, T . The matrix form becomes:
22
The fractional step method could be applied here.
T gradient method can be applied to solve and λ. Consider a special case that
∆ ∆ , then . T can be regard as:
T 2
So in this case, .
23
4 Numerical Results
4.1 Error estimate of fractional step method
Consider the incompressible Navier-Stokes equations consisting of the momentum and continuity equations. There exists an exact solution that satisfies the Navier-Stokes equations. They are
, , 2 ⁄
, , 2 ⁄
, , 2 2 4 ⁄ ⁄ 4
Use this data, we can know the initial and boundary conditions, and compare the numerical solution with the exact solution. We set the computational domain Ω 0,1 0,1 , the computational time T 1, and the Reynolds number Re 40.
To check the order of accuracy of the fractional step method, let ∆t 1 400⁄ , 1 800⁄ , 1 1600⁄ , 1 3200⁄ , and 1 6400⁄ , respectively. N 3. Suppose that ∆ ∆ , and let 40∆t. We define the maximum error by subtracting the numerical solution from the exact solution and take infinite norms. The order is known by taking on the ratio of the errors. As shown in table. 1, the fractional step method has second-order accuracy in time.
Table 1
The maximum errors from different ∆t
∆t maximum error order 1 400⁄ 1 10⁄ 3.2781e-004 2.47 1 800⁄ 1 20⁄ 5.9040e-005 1.97 1 1600⁄ 1 40⁄ 1.5094e-005 2.01 1 3200⁄ 1 80⁄ 3.7449e-006 1.96 1 6400⁄ 1 160⁄ 9.6555e-007 －
4.2 The flow past a cylinder
Consider a situation that the flow past a cylinder, where the diameter of the cylinder is . The matrix factorization method can be applied here to simulate this situation.
In this test, the computational domain is set by Ω 0,16 0,8 . The Reynolds number Re 100 and 200, respectively. We choose ∆t 1 640⁄ and ∆x ∆y 1 16⁄ . The initial condition is given by | 0 and
| 0. The boundary condition is shown in Fig.6. To construct the cylinder, we put the center of the cylinder at the position 4,4 , and let the diameter
1. Choose the number of the markers 64, then ∆ ⁄ . Notice 64 that ∆ 0.785 .
Fig. 6. The boundary condition and the computational domain
The simulations which the flow past a cylinder at Reynolds number 100 and 200 show the periodic vortex shedding. In figure 7, we can see the periodic vortex shedding in the vorticity contours. The vorticity of the flow is defined as
.
24
Fig. 7. The vorticity contours at Reynolds number 100(left) and 200(right) at the time T 60
There are three quantities, the drag coefficient , lift coefficient , and the Strouhal number . We can use the three quantities to compare this simulation with others. The three quantities are defined as:
⁄2
⁄2
where and are the drag force and lift force, 1 in this simulation, is the frequency of the vortex shedding. and can be obtain in [11], which are
, , δ d d
Ω
and are approximated by
·
·
In figure 8 and 9, we can observe the periodic vortex shedding. In table 2 and 3, we compare the three quantities with the previous numerical results which refer to [9, 11, 12, 13, 14, 15] at Reynolds number Re 100 and 200.
25
Fig. 8. The time evolution of drag(left) and lift(right) coefficients at Reynolds number Re 100
26
Fig. 9. The time evolution of drag(left) and lift(right) coefficients at Reynolds number Re 200
Table 2
The comparisons of lift and drag coefficients and Strouhal number of Re 100 Present Lai and Peskin[11] Kim et al.[12] Silva et al.[13]
1.64 1.44 1.33 1.39 0.177 0.165 - 0.16 0.40 0.33 0.32 -
Table 3
The comparisons of lift and drag coefficients and Strouhal number of Re 200 Present Taira and Colonius [9] Linnick and Fasel[14] Liu et al[15]
1.56 1.36 1.34 1.31 0.206 0.197 0.197 0.192 0.72 0.69 0.69 0.69
4.3 The flow past two cylinders
In this test, we simulate the flow past two cylinders at the Reynolds number Re 200. Let the diameter of the two cylinders are the same, and set to be
1. Similar to the previous case, we set the computation domain Ω 0,16 0,8 , ∆t 1 640⁄ and ∆x ∆y 1 16⁄ . The initial and boundary conditions are as before, too. To construct the cylinder, we put the centers of this two cylinders at the position 4, 2.5 and 4, 5.5 . We also choose ∆
⁄ , so the number of the markers 64 128. The time evolution of drag and lift coefficients and the vorticity contours are shown in figure 10 and 11. We can observe that the drag coefficients of the two cylinders are similar, and the lift coefficients of the two cylinders are symmetric.
Fig. 10. The time evolution of drag(left) and lift(right) coefficients of the upper(up) and lower (down) cylinders.
27
(a) (b)
(c) (d)
(e) (f)
(g) (h)
(i) (j)
(k) (l)
Fig. 11. The vorticity contours when the flow past two cylinders : (a)t=7, (b)t=8, (c)t=9, (d)t=10, (e)t=11, (f)t=12, (g)t=30, (h)t=31, (i)t=32, (j)t=33, (k)t=34, (l)t=35.
28
4.4 The flow past a wing
Here, we simulate the flow past through a wing of the airplane at the Reynolds number Re 200. We use the same computational domain, mesh, initial condition and boundary conditions as before, but the immersed object. See figure 12, the shape of the immersed object is a airfoil, which is a thin winglike structure. We rotate the wing by an angle, which is 30° here, as shown in figure 13. The periodic vortex shedding also can be observed behind the wing in figure 14.
Fig. 12. The shape of the immersed object
Fig. 13. The placement of the wing
29
(a) (b)
(c) (d)
(e) (f)
(g) (h)
(i) (j)
Fig. 14. The vorticity contours when the flow past a wing : (a)t=51, (b)t=52, (c)t=53, (d)t=54, (e)t=55, (f)t=56, (g)t=57, (h)t=58, (i)t=59, (j)t=60.
30
4.5 The flow past an oscillating cylinder
Consider the same simulation as the flow past a cylinder at Reynolds number Re 100, but the cylinder is moving here. We impose the velocity , at the boundary as , 0, 0.14 cos 2 , where is the frequency that the cylinder oscillates. That is, the cylinder oscillates vertical to the stream. Here we choose 2 . The time evolution of drag coefficient and the vorticity contours are shown in figure 15 and 16. Compare with figure 8, we can see that the frequency of the vortex shedding is influenced by the oscillation of the cylinder after the cylinder moves. In table 4, we compare the drag and lift coefficients with the previous numerical results which refer to [17, 18].
Fig. 15. The time evolution of the lift(right) coefficients
Table 4
The comparisons of the lift and drag coefficients
Present Su et al.[17] Hurlbut et al.[18]
1.84 1.70 1.68 1.75 0.97 0.95
31
(a) (b)
(c) (d)
(e) (f)
(g) (h)
Fig. 16. The vorticity contours when the flow past a oscillating cylinder : (a)t=4.21875, (b)t=4.921875, (c)t=5.625, (d)t=6.328125, (e)t=7.03125, (f)t=7.734325, (g)t=8.4375, (h)t=9.140625.
32
33
5 Conclusion
In this thesis, we use the matrix factorization method introduced by Taira and Colonius [9] to simulate the flow past as immersed object. In the numerical result, we see that this method can handle the immersed object with complex shape, and which the immersed object is moving, even the two or more objects. It is useful in the engineering applications. In Section 4.4, we simulate the situation when the flow past a wing. The flow produces vortex shedding behind the wing.
Here we suppose that the mesh widths of each cell are the same, that is, we let ∆ ∆ for all , . In fact, the matrix factorization method can use the different grid size of each cell. If the higher accuracy is desired, the grid size near the immersed object has to be small. We can adapt the code to handle the different mesh widths, to improve the accuracy.
34
References
[1] A.J. Chorin, “Numerical solution of the Navier–Stokes equations.” Math. Comput. Vol. 22, pp. 745–762, 1968.
[2] J. Kim, P. Moin, “Application of a fractional-step method to incompressible Navier–Stokes equations.” J. Comput. Phys.
Vol. 59, pp. 308–323, 1985.
[3] J.B. Perot, “An analysis of the fractional step method.” J. Comput. Phys. 108 51–58, 1993.
[4] R. Codina, “Pressure stability in fractional step finite element methods for incompressible flows.” J. Comput. Phys. Vol.
170, pp. 112–140, 2001.
[5] F.H. Harlow, J.E. Welsh, “Numerical calculation of time-dependent viscous incompressible flow of fluid with a free
surface.” Phys. Fluids. Vol. 8, pp. 2181-2189, 1965.
[6] D.L. Brown, R. Cortez, M.L. Minion, “Accurate projection methods for the incompressible Navier–Stokes equations.” J.
Comput. Phys. Vol. 168, pp. 464–499, 2001.
[7] C.S. Peskin, “Flow patterns around heart valves: a numerical method.” J. Comput. Phys. Vol. 10, pp. 252–271, 1972.
[8] C.S. Peskin, “The immersed boundary method”. Acta Numer. Vol. 11, pp. 479-517, 2002.
[9] K. Taira, T. Colonius, “The immersed boundary method: A projection approach.” J. Comput. Phys. Vol. 225, pp.
2118-2137, 2007.
[10] A.M. Roma, C.S. Peskin, M.J. Berger, “An adaptive version of the immersed boundary method.” J. Comput. Phys. Vol.
153, pp. 509–534, 1999.
[11] M. Lai, C.S. Peskin, “An immersed boundary method with formal second-order accuracy and reduced numerical
viscosity.” J. Comput. Phys. Vol. 160, pp. 705–719, 2000.
[12] J. Kim, D. Kim, H. Choi, “An immersed-boundary finite-volume method for simulations of flow in complex geometries.”
J. Comput. Phys. Vol. 171, pp. 132-150, 2001.
[13] A. L. F. Lima E Silva, A. Silveira-Neto, J. J. R. Damasceno, “Numerical simulation of two-dimensional flows over a
circular cylinder using the immersed boundary method.” J. Comput. Phys. Vol. 189, pp. 351-370, 2003.
[14] M.N. Linnick, H.F. Fasel, “A high-order immersed interface method for simulating unsteady incompressible flows on
irregular domains.” J. Comput. Phys. Vol. 204, pp. 157–192, 2005.
[15] C. Liu, X. Zheng, C.H. Sung, “Preconditioned multigrid methods for unsteady incompressible flows.” J. Comput. Phys.
Vol. 39, pp. 35–57, 1998. 1
[16] R. Témam, “Sur l’approximation de la solution des équations de Navier–Stokes par la méthode des pas fractionnaires.”
Arch. Rat. Mech. Anal. Vol. 32 (2), pp. 135–153, 1969.
[17] S. Su, M. Lai, C. Lin, “An immersed boundary technique or simulating complex flows with rigid boundary.” Comput.
Fluids. Vol. 36, pp. 313-324, 2007.
[18] S.E. Hurlbut, M.L. Spaulding, F.M. White, “Numerical Solution for Laminar Two Dimensional Flow About a Cylinder
Oscillating in a Uniform Stream”, Trans ASME, J. Fluids Eng. Vol. 104, pp. 214-221, 1982.
[19] W. Chang, F. Giraldo, B. Perot, “Analysis of an exact fractional step method.” J. Comput. Phys. Vol. 180, pp. 183–199,
2002.