三、 Matrix Factorization Method
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.