Published: 30 September 2013

Extension of piecewise exact solution method for two- and three-dimensional fluid flows

Hyoseob Kim1
Jongsup Ahn2
Changhwan Jang3
Jaeyoull Jin4
Byungsoon Jung5
1, 2Kookmin University, Seoul, Korea
3Korean Intellectual Property Office, Daejeon, Korea
4Korea Institute of Ocean Science & Technology (KIOST), Ansan, Korea
5Sea & River Consultant Inc., Seoul, Korea
Corresponding Author:
Changhwan Jang
Views 32
Reads 15
Downloads 1211

Abstract

Extended forms of a pseudo-numerical scheme for advection terms in fluid momentum equations are proposed here. The fact that analytic solution exists for the Burgers equation, if velocity distribution in space is straight for one-dimensional flow, was shown by Jang et al. Analytic solution also exists for two- or three-dimensional fluid flows, if the velocity components in two- or three-direction are linearly distributed in space, and the existing piecewise exact solution method is extended for two- and three-dimensions here. The analytic solution is adopted for computation of the advecting property of fluid momentum in two- or three-dimensional directions. This method produces zero numerical error during one time increment so that it is distinguished from any other numerical scheme which produces small or large numerical error within one time increment. The behavior of the new scheme is demonstrated for two- and three-dimensional examples. The nonlinear modifications of velocity profiles towards singularity with time progress are well simulated for three test cases. The computed maximum relative errors for a given condition for one-, two-, and three-dimensions become larger as the number of dimension increases. The scheme is believed to work well for two- and three-dimensional flows.

1. Introduction

Fluid motion is governed by conservation equations including momentum conservation equation. Not only the advection equation but also the advection-diffusion equation has been intensively used to describe the flow of fluids or other material in engineering problems [1, 2, 3].

Fundamentally, the one-dimensional advection-diffusion equation has a form:

1
ut+aut=v2ux2,

which is reduced to the advection-only equation:

2
ut+aux=0,

where a is the transporting celerity, and u is a physical property. The inviscid Burgers equation is a specific form of the above advection equation, when a=u:

3
ut+uux=0.

From mathematical point of view the momentum equations in three dimensions are composed of several terms, i.e, advection terms, diffusion terms, and the others. Fractional step method [4, 5] or operator splitting method [6] is often adopted in numerical solution procedure due to its convenience and simplicty. The one-dimensional invisicid Burgers equation can be written in conservative form as:

4
ut+F(u)x=0,
5
Fu=u22,

where x is the spatial axis, t is the time, and u is the velocity in the x axis. The velocity field u(x,t) at the initial time (t=0) is u(x,0) and is called f(x0). The analytic solution of Eq. (4) for the given initial condition is:

6
uLx0,t=ux0,0=fx0,

where uL is the velocity for Langrangean position. The Eulerian description of the velocity becomes:

7
x=x0+uLt=x0+fx0t,
8
ux,t=ux0+fx0t,t=uLx0,t=uLx0,0=fx0.

The above advection equation is called the inviscid Burgers equation [7, 8]. The two- or three-dimensional Burgers equation in a Lagrangean form is:

9
DV/Dt=0,

where V is the velocity vector of fluid flow, V=(u,v,w), (u,v), or (u) for three-dimensional, two-dimensional, or one-dimensional flow, respectively, where u, v and w are the velocity components in the x, y and z directions, respectively. Eulerian forms of the Burgers equations are then:

10
ut+uux+vuy+wuz=0,
vt+uvx+vvy+wvz=0,
wt+uwx+vwy+wwz=0.

The problem is that it is not always possible to get the solution in an explict form for an arbitrary initial function f. We need adquate nemerical methods which are good for fixed grid system [9, 10]. One of frequently used numerical schemes to solve the advection equation is the upwind scheme [9]:

11
uin+1+uinΔt=-u-inΔxui+1n-uin or uin-ui-1n,

where superscript n or n+1 means the time step n or n+1, subscript i means the spatial grid number i, u-in is a representative velocity at the index number i at time step n, which could either be uin, ui+1n, ui-1n, (uin+ui+1n)/2 or (ui-1n+uin)/2. If we use the difference equation for the Burgers equation in the conservative form, the difference equation becomes:

12
uinΔt=-12Δxui+1n2-uin2 or uin2-ui-1n2.

This virtually takes a representative velocity as (uin+ui+1n)/2 or (ui-1n+uin)/2 depending on the flow direction. Jang et al. [12] proposed a useful piecewise exact solution method (PESM) to solve the one-dimensional advection equation, Eq. (3). If the velocity is linearly distributed in space as:

13
u=x+ct.

Then the above equation itself satisfies the advection equation, Eq. (3), if conditions are well-posed. Eq. (13) may have a singular solution, if the velocity gradient becomes infinity as time goes by. This corresponds to shock wave in physical domain [11]. The solution at a fixed point of index i becomes:

14
uin+1=Δx(uin-ui-1n)Δt+Δxuin.

However, if the governing momentum equation is to be solved with the splitting method, the analytical solution for some long period may not be very useful, because not only invisid Burgers equation but also the diffusion equation or other equation with other terms should be alternately solved every time step. Therefore the analytical solution for the invisid Burgers equation may be useful for each time increment only.

Interesting feature of this solution is that the orgin of t, t0 is not a priori, but should be found, and the solution is not temporarilly linear from the Eulerion point of view. If a fixed grid is used, and two velocities, at two points xi and xi-1 at time t0 are:

15
uin=xi+ct0, ui-1n=xi-Δx+ct0.

We can remove c from Eq. (15), and get uin+1 as:

16
uin+1=xi+ct0+Δt=Δxuin(uin-ui-1n)Δt+Δx.

If the velocity at xi is negative, the velocity at xi at new time step is:

17
uin+1=Δxuin(ui+1n-uin)Δt+Δx.

The above PESM for one-dimensional advection equation was compared with Godunov's method [11], and demonstrated higher accuracy relative to Godunov's method. Since the above solution is the analytic solution, it does not produce any numerical error during a time increment, but errors develop during the spatial interpolation. The PESM was compared with a few existing numerical schemes, and showed unconditional stability for CFL number larger than 1.0 by using the following extended solution:

18
uin+1=Δx{(m+1)ui-mn-mui-m-1n}(ui-mn-ui-m-1n)Δt+Δx, uin0 and ui-mnΔtmΔx,uin+1=Δx{(m+1)ui-mn-mui-m+1n}(ui-m+1n-ui-mn)Δt+Δx, uin<0 and ui-mnΔt<mΔx,

where m is an integer. Jang et al. suggested the above solutions could be used for the first-order portion of any other existing higher order schemes, e.g. Fromm's or Lax-Wenddroff scheme. However, Jang et al.’s method has been confined to one-dimensional form.

2. Extension of PESM for two- and three-dimensional cases

We are extending Jang et al.’s PESM for higher dimension here, i.e. for two-dimensional and three-dimensional advection equations [12]. Starting from one-dimensional problem, we adopt coefficients for spacial distribution of velocity, not involving t0:

19
uin=axi+b,
ui-1n=axi-Δx+b,

where b implicity include the time effect. We can also use the invariant Lagrangean advection property of the velocity component. Furthermore we also know that the solution must be spatially linear:

20
uin+1=dxi+e.

We have two boundary values at two points at the advanced time:

21
ui*n+1=uin at xi*=xi+uinΔt,
u(i-1)*n+1=ui-1n at x(i-1)*=xi-0-Δx+ui-1nΔt,

where superscript * means the new position of the particle which was at grid point i at the previous time. Inserting Eq. (21) into Eq. (19), we obtain the solution for uin+1 which is identical to Eq. (16). Now we expand this approach to the two-dimensional advection equations:

22
ut=uux+vuy,
vt=uvx+vvy.

A rectangular grid is considered at this stage. Planar distribution of the two velocity components at an instance are assumed as:

23
u=a1x+a2y+a3t, v=b1x+b2y+b3t,
24
ux,y,t0=a1x+a2y+a3t0=c1x+c2y+c3, vx,y,t0=d1x+d2y+d3.

The coefficients a1, a2, a3, b1, b2, b3 are assumed to satisfy Eq. (22). We have six boundary values at three triangle corner points A, B and C of a triangle at time t0. For example, if ui,jn0 and vi,jn0, then the three corner points will move to other positions i, j*, i-1, j* and i, j-1* at t0+Δt:

25
ui,j*n+1=c1xi,j*+c2yi,j*+c3,
ui-1,j*n+1=c1xi-1,j*+c2yi-1,j*+c3,
u(i,j-1)*n+1=c1x(i,j-1)*+c2y(i,j-1)*+c3,
26
vi,j*n+1=d1xi,j*+d2yi,j*+d3,
vi-1,j*n+1=d1xi-1,j*+d2yi-1,j*+d3,
v(i,j-1)*n+1=d1x(i,j-1)*+d2y(i,j-1)*+d3,

where:

27
xi,j*n+1=ui,jnΔt,
y(i,j)*n+1=v(i,j)nΔt,
28
xi-1,j*n+1=-Δx+ui-1,jnΔt,
y(i-1,j)*n+1=v(i-1,j)nΔt,
29
xi,j-1*n+1=ui,j-1nΔt,
y(i,j-1)*n+1=-Δy+vi,j-1nΔt.

Fig. 1Position movement of three corner points of triangle for 2-dimensional case

Position movement of three corner points of triangle for 2-dimensional case

Then, at time t0+Δt the six unknowns c1, c2, c3, d1, d2, d3 are obtained from the above six equations. For different flow directions other three corner points are chosen, i.e. (i, j), (i-1, j), (i, j+1) for negative ui,j and positive vi,j, and (i, j), (i+1, j), (i, j+1) for positive ui,j and negative vi,j. Eqs. (27) to (29) can be expressed as:

30
M1c=u,
M2d=v,

where M1, M2 are the coefficient matrices, c and d are the unknown vectors c1, c2, c3 and d1, d2, d3, respectively; u and v are known boundary value vector ui,j*n+1,ui-1,j*n+1,ui,j-1*n+1 and vi,j*n+1,vi-1,j*n+1,vi,j-1*n+1, respectively. Then:

31
c=M1-1u
d=M2-1v.

Thus ui,jn+1=c3 and vi,jn+1=d3 are new velocity components at (i, j), i.e. x=y=0. Next we expand this approach to the three-dimensional advection equations:

32
ut=uux+uuy+uuz,
vt=vvx+vvy+vvz,
wt=wwx+wwy+wwz.

A rectangular parallelepiped grid is considered at this stage. Planar distribution of the three velocity components at an instance are assumed as:

33
u=a1x+a2y+a3z+a4t,v=b1x+b2y+b3z+b4t,w=c1x+c2y+c3z+c4t,
34
ux,y,z,t0=d1x+d2y+d3z+d4,
vx,y,z,t0=e1x+e2y+e3z+e4,
wx,y,z,t0=f1x+f2y+f3z+f4.

The coefficients a1, a2, a3, a4, b1, b2, b3, b4, c1, c2, c3, c4 are assumed to satisfy Eq. (32). We have twelve boundary values at four corner points A, B, C and D of the quadrangular pyramid at time t0. For example, if ui,jn0, vi,jn0 and wi,jn0 then the four corner points will take other positions at t0+Δt:

35
ui,j,k*n+1=c1xi,j,k*+c2yi,j,k*+c3zi,j,k*+c4,
u(i-1,j,k)*n+1=c1x(i-1,j,k)*+c2y(i-1,j,k)*+c3z(i-1,j,k)*+c4,
ui,j-1,k*n+1=c1xi,j-1,k*+c2yi,j-1,k*+c3zi,j-1,k*+c4,
u(i,j,k-1)*n+1=c1x(i,j,k-1)*+c2y(i,j,k-1)*+c3z(i,j,k-1)*+c4,
36
vi,j,k*n+1=d1xi,j,k*+d2yi,j,k*+d3zi,j,k*+d4,
vi-1,j,k*n+1=d1xi-1,j,k*+d2yi-1,j,k*+d3zi-1,j,k*+d4,
vi,j-1,k*n+1=d1xi,j-1,k*+d2yi,j-1,k*+d3zi,j-1,k*+d4,
v(i,j,k-1)*n+1=d1x(i,j,k-1)*+d2y(i,j,k-1)*+d3z(i,j,k-1)*+d4,
37
wi,j,k*n+1=a1xi,j,k*+a2yi,j,k*+a3zi,j,k*+a4,
wi-1,j,k*n+1=a1xi-1,j,k*+a2yi-1,j,k*+a3zi-1,j,k*+a4,
wi,j-1,k*n+1=a1xi,j-1,k*+a2yi,j-1,k*+a3zi,j-1,k*+a4,
w(i,j,k-1)*n+1=a1x(i,j,k-1)*+a2y(i,j,k-1)*+a3z(i,j,k-1)*+a4,

where:

38
xi,j,k*n+1=ui,j,knΔt,
yi,j,k*n+1=vi,j,knΔt,
z(i,j,k)*n+1=w(i,j,k)nΔt,
39
xi-1,j,k*n+1=-Δx+ui-1,j,knΔt,
yi-1,j,k*n+1=vi-1,j,knΔt,
zi-1,j,k*n+1=wi-1,j,knΔt,
40
xi,j-1,k*n+1=ui,j-1,knΔt,
yi,j-1,k*n+1=-Δy+vi,j-1,knΔt,
z(i,j-1,k)*n+1=w(i,j-1,k)nΔt,
41
xi,j,k-1*n+1=ui,j,k-1nΔt,
yi,j,k-1*n+1=vi,j,k-1nΔt,
z(i,j,k-1)*n+1=-Δz+wi,j,k-1nΔt.

Fig. 2Position movement of four corner points of quadangular pyramid for 3-dimensional case

Position movement of four corner points of quadangular pyramid for 3-dimensional case

Then, at time t0+Δt the twelve unknown coefficients d1, d2, d3, d4, e1, e2, e3, e4, f1, f2, f3, f4 are obtained from the above twelve equations. For different flow directions from the above case, other four corner points are chosen, see Table 1. Eqs. (35) to (37) can be expressed as:

42
M1d=u,
M2e=v,
M3f=w,

where M1,M2,M3 are the coefficient matrices, d,e and f are the unknown vectors d1, d2, d3, d4, e1, e2, e3, e4, f1, f2, f3, f4 respectively, and u, v and w are known boundary value vectors ui,j,k*n+1,ui-1,j,k*n+1,ui,j-1,k*n+1,ui,j,k-1*n+1, vi,j,k*n+1,vi-1,j,k*n+1,vi,j-1,k*n+1,vi,j,k-1*n+1 and wi,j,k*n+1, wi-1,j,k*n+1,wi,j-1,k*n+1,wi,j,k-1*n+1 respectively. Then:

43
d=M1-1u
e=M2-1v,
f=M3-1w

Thus, ui,j,kn+1=d3,vi,j,kn+1=e3 and wi,j,kn+1=f3 are new velocity components at (i, j, k), i.e. x=y=z=0.

Table 1Selection of four corner points for different signs of u, v and w

ui,j,k
vi,j,k
wi,j,k
Four points
+
+
-
(i, j, k), (i-1, j, k), (i, j-1, k ), (i, j, k+1)
+
-
+
(i, j, k), (i-1, j, k), (i, j+1, k ), (i, j, k-1)
+
-
-
(i, j, k), (i-1, j, k), (i, j+1, k ), (i, j, k+1)
-
+
+
(i, j, k), (i+1, j, k), (i, j-1, k ), (i, j, k-1)
-
+
-
(i, j, k), (i+1, j, k), (i, j-1, k ), (i, j, k+1)
-
-
+
(i, j, k), (i+1, j, k), (i, j+1, k ), (i, j, k-1)
-
-
-
(i, j, k), (i+1, j, k), (i, j+1, k ), (i, j, k+1)

3. Examinations of present extention

Now we apply the present PESM extended for two- and three-dimensional fluid flow to two simple two- and three-dimensional situations. The initial distribution of the velocity field of the the two- dimensional flow is given as:

44
u=221+exp-x-502-y-502500,
v=221+exp-x-502-y-502500.

And the initial velocity field for a simple three-dimensional situation is:

45
u=331+exp-x-502-y-502-z-502500,
v=331+exp-x-502-y-502-z-502500,
w=331+exp-x-502-y-502-z-502500.

Fig. 3Speed contour at two instants for two dimensional test flow: t=0

Speed contour at two instants  for two dimensional test flow: t=0

Fig. 4Speed contour at two instants for two dimensional test flow: t=15

Speed contour at two instants  for two dimensional test flow: t=15

Fig. 5Speed contour at two instants for three dimensional test flow: t= 0

Speed contour at two instants  for three dimensional test flow: t= 0

Fig. 6Speed contour at two instants for three dimensional test flow: t= 15

Speed contour at two instants  for three dimensional test flow: t= 15

Fig. 7Velocity profiles along x axis for one-dimensional flow

Velocity profiles along x axis for one-dimensional flow

a) Profile

Velocity profiles along x axis for one-dimensional flow

b) Zoomed view

Fig. 8Velocity profiles along x axis for two-dimensional flow

Velocity profiles along x axis for two-dimensional flow

a) Profile

Velocity profiles along x axis for two-dimensional flow

b) Zoomed view

One-dimensional case is also used for comparision with the two- and three-dimensional cases. The initial velocity distribution of the one-dimensional case is:

46
u=1+exp-x-502500.

The application results of this method to the two-dimensional and three-dimensional progressive flood cases are shown in Figs. 3, 4, 5 and 6, respectively.

Computed profiles of the speed along central lines at two instants for one-, two-, and three-dimensional cases are shown in Figs. 7 to 9. Figs. 7 to 9 show that even though the PESM is the exact solution for each Δt, the method still produces numerical diffusion due to the process of spacial interpolation every time step, see Jang et al. [12]. The computed distributions of speed along central lines are close to analytical solution for one-, two-, and three-dimensional cases, respectively. The maximum relative errors of the computed speed for one-, two-, and three-dimensional cases are shown in Table 2. The relative error becomes larger as the number of dimensions increases, which means the spacial interpolation produces high diffusion. Nontheless the expanded PESM works well for two-and three-dimensional problems according to the present examinations.

Fig. 9Velocity profiles along x axis for three-dimensional flow

Velocity profiles along x axis for three-dimensional flow

a) Profile

Velocity profiles along x axis for three-dimensional flow

b) Zoomed view

Table 2Comparison of relative errors for different dimensions

Dimension
Analytic solution
PESM result
Relative error (%)
One
2.0000
1.9931
0.3457
Two
2.0000
1.9844
0.7799
Three
2.0000
1.9812
0.9408

4. Conclusions

Previous PESM for one-dimensional flow was expanded here for two-, and three-dimensional fluid flows. The PESM takes exact solution within one time step and therefore it is irrelevant to any numerical error, although it is not free from interpolation error due to the non-moving property of grid computational. The PESM secures unconditional stability, which is important for practical use in numerical modeling workes. Algorithms for two- or three-dimensional advection equations are extensions of that for one-dimensional advection equation. The boundery values at corners of a triangle or a quadrangular pyramid at time t0 preserved at time t0+Δt according to the characteristics of the Lagrangean advection equation. The method was applied to simple situations, and the results are satisfactory. The diffusion due to spacial interpolation become larger as the dimension becomes larger. The PESM is believed to work well for arbitrary number of dimensions.

References

  • Morton K. W., Mayers D. F. Numerical Solution of Partial Differential Equations. Cambridge University Press, Cambridge, UK, 2005, p. 86-150.
  • Ibragimov N. H. Handbook of Lie Groups to Analysis of Differential Equations. Vol. 1, CRC Press, Boca Raton, USA, 1994.
  • Polyanin A. D., Zaitsev V. F. Handbook of Nonlinear Partial Differential Equations. Chapman and Hall / CRC, Boca Raton, USA, 2004.
  • Perot J. B. An analysis of the fractional step method. J. of Comput. Phys., Vol. 108, 1993, p. 51-58.
  • Armfiled S., Street R. Modified fractional-step methods for the Navier-Stokes equations. ANZIAM J., Vol. 45, 2004, p. 364-377.
  • Karlsen K. H., Risebro N. H. An operator splitting method for nonlinear convection-diffusion equation. Number. Math., Vol. 77, No. 3, 1994, p. 365-382.
  • Burgers J. M. The Nonlinear Diffusion Equation. Reidel, Dordreht, Netherlands, 1974.
  • Burgers J. M. A mathematical model illustrating the theory of turbulence. Adv. Appl. Meth., Vol. 1, 1948, p. 171-199.
  • Welch J. E., Harlow F. H., Shannon J. P., Daly B. J. The MAC Method. Los Alamos Scientific Lab. Rept. LA-3425, Los Alamos, New Mexico, USA, 1966.
  • Fromm J. E. A method for reducing dispersion in convective difference schemes. J. Comput. Phys., Vol. 3, 1968, p. 176-189.
  • Godunov S. K. A difference scheme for numerical solution of discontinuous solution of hydrodynamic equations. Math. Sbornik, Vol. 47, 1989, p. 271-306.
  • Jang C., Kim H., Choi C., Kim J. Piecewise exact solution of nonlinear momentum conservation equation with unconditional stability for time increment. Journal of Vibroengineering, Vol. 14, 2012, p. 1041-1051.
  • Lax P. D., Wendroff B. Systems of conservation laws. Comm. Pure Appl. Math., Vol. 13, No. 2, 1960, p. 217-237.
  • Sarra S. A. The method of characteristics with application to conservation laws. J. of Online Math. and Its Appl., Vol. 3, 2003.
  • Roe P. L. Approximate Riemann solvers, parameter vectors, and difference schemes. J. of Comput. Phys., Vol. 43, 1981, p. 357-372.
  • Harten A. One class of high resolution total-variation-stable finite difference schemes. SIAM J. Numer. Anal., Vol. 21, 1984, p. 1-23.
  • Thomee V. A stable difference scheme for the mixed boundary value problem for a hyperbolic, first order system in two dimensions. J. Soc. Indust. Appl. Math., Vol. 10, 1962, p. 229-245.
  • Flather R. A., Heaps N. S. Tidal computations for Morecamble bay. Geophys. J. R. Astr. Soc., Vol. 42, 1975, p. 489-517.
  • Van Leer B. Monotonicity and conservation combined in a second-order scheme. J. Comput. Phys., Vol. 14, 1974, p. 361-370.

About this article

Received
10 May 2013
Accepted
04 September 2013
Published
30 September 2013
Keywords
three-dimensional flows
PESM
advection
Burgers equation
Acknowledgements

This paper is a part of the results of a Research Project titled “Marine Environmental Prediction System” and “Developement of Coastal Erosion Control Technology” funded by the Ministry of Land, Tranport and Maritime Affairs, Korea, and the Kookmin University Grant, 2013.