Extension of piecewise exact solution method for two- and three-dimensional fluid 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:

which is reduced to the advection-only equation:

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$:

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:

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\left(x_0\right)$. The analytic solution of Eq. (4) for the given initial condition is:

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

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

where $\mathbf{V}$ is the velocity vector of fluid flow, $\mathbf{V}=(u,v,w)$, $(u,v)$, or $\left(u\right)$ 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:

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]:

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

This virtually takes a representative velocity as $(u_i^n+u_{i+1}^n)/2$ or $(u_{i-1}^n+u_i^n)/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:

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:

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$, $t_0$ 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 $x_i$ and $x_{i-1}$ at time $t_0$ are:

We can remove $c$ from Eq. (15), and get $u_i^{n+1}$ as:

If the velocity at $x_i$ is negative, the velocity at $x_i$ at new time step is:

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:

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 $t_0$:

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:

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

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 $u_i^{n+1}$ which is identical to Eq. (16). Now we expand this approach to the two-dimensional advection equations:

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

The coefficients $a_1$, $a_2$, $a_3$, $b_1$, $b_2$, $b_3$ 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 $t_0$. For example, if $u_{i,j}^n\ge 0$ and $v_{i,j}^n\ge 0$, then the three corner points will move to other positions ${\left(i,j\right)}^{*}$, ${\left(i-1,j\right)}^{*}$ and ${\left(i,j-1\right)}^{*}$ at $t_0+\Delta t$:

where:

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

Then, at time $t_0+\Delta t$ the six unknowns $c_1$, $c_2$, $c_3$, $d_1$, $d_2$, $d_3$ 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 $u_{i,j}$ and positive $v_{i,j}$, and $(i,j)$, $(i+1,j)$, $(i,j+1)$ for positive $u_{i,j}$ and negative $v_{i,j}$. Eqs. (27) to (29) can be expressed as:

where $M_1$, $M_2$ are the coefficient matrices, $c$ and $d$ are the unknown vectors $\left(c_1,c_2,c_3\right)$ and $\left(d_1,d_2,d_3\right)$, respectively; $u$ and $v$ are known boundary value vector $\left(u_{{\left(i,j\right)}^{*}}^{n+1},u_{{\left(i-1,j\right)}^{*}}^{n+1},u_{{\left(i,j-1\right)}^{*}}^{n+1}\right)$ and $\left(v_{{\left(i,j\right)}^{*}}^{n+1},v_{{\left(i-1,j\right)}^{*}}^{n+1},v_{{\left(i,j-1\right)}^{*}}^{n+1}\right)$, respectively. Then:

Thus $u_{i,j}^{n+1}=c_3$ and $v_{i,j}^{n+1}=d_3$ are new velocity components at $(i,j)$, i.e. $x=y=0$. Next we expand this approach to the three-dimensional advection equations:

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

The coefficients $a_1$, $a_2$, $a_3$, $a_4$, $b_1$, $b_2$, $b_3$, $b_4$, $c_1$, $c_2$, $c_3$, $c_4$ 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 $t_0$. For example, if $u_{i,j}^n\ge 0$, $v_{i,j}^n\ge 0$ and $w_{i,j}^n\ge 0$ then the four corner points will take other positions at $t_0+\Delta t$:

where:

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

Then, at time $t_0+\Delta t$ the twelve unknown coefficients $d_1$, $d_2$, $d_3$, $d_4$, $e_1$, $e_2$, $e_3$, $e_4$, $f_1$, $f_2$, $f_3$, $f_4$ 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:

where $M_1\text{,}$$M_2\text{,}$$M_3$ are the coefficient matrices, $d\text{,}$$e$ and $f$ are the unknown vectors $\left(d_1,d_2,d_3,d_4\right)$, $\left(e_1,e_2,e_3,e_4\right)$, $\left(f_1,f_2,f_3,f_4\right)$ respectively, and $u$, $v$ and $w$ are known boundary value vectors $\left(u_{{\left(i,j,k\right)}^{*}}^{n+1},u_{{\left(i-1,j,k\right)}^{*}}^{n+1},u_{{\left(i,j-1,k\right)}^{*}}^{n+1},u_{{\left(i,j,k-1\right)}^{*}}^{n+1}\right)$, $\left(v_{{\left(i,j,k\right)}^{*}}^{n+1},v_{{\left(i-1,j,k\right)}^{*}}^{n+1},v_{{\left(i,j-1,k\right)}^{*}}^{n+1},v_{{\left(i,j,k-1\right)}^{*}}^{n+1}\right)$ and $\left(w_{{\left(i,j,k\right)}^{*}}^{n+1},\mathrm{}{w}_{{\left(i-1,j,k\right)}^{*}}^{n+1},w_{{\left(i,j-1,k\right)}^{*}}^{n+1},w_{{\left(i,j,k-1\right)}^{*}}^{n+1}\right)$ respectively. Then:

Thus, $u_{i,j,k}^{n+1}=d_3\text{,}$$v_{i,j,k}^{n+1}=e_3$ and $w_{i,j,k}^{n+1}=f_3$ are new velocity components at $(i,j,k)\text{,}$ i.e. $x=y=z=0$.

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:

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

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

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

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

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

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

a) Profile

b) Zoomed view

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

a) Profile

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:

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 $\Delta 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

a) Profile

b) Zoomed view

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 $t_0$ preserved at time $t_0+\Delta 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.