A DIMENSION SPLITTING AND CHARACTERISTIC PROJECTION METHOD FOR THREE-DIMENSIONAL INCOMPRESSIBLE FLOW

. A dimension splitting and characteristic projection method is proposed for three-dimensional incompressible ﬂow. First, the characteristics method is adopted to obtain temporal semi-discretization scheme. For the remaining Stokes equations we present a projection method to deal with the in- compressibility constraint. In conclusion only independent linear elliptic equations need to be calculated at each step. Furthermore on account of splitting property of dimension splitting method, all the computations are carried out on two-dimensional manifolds, which greatly reduces the diﬃculty and the compu- tational cost in the mesh generation. And a coarse-grained parallel algorithm can be also constructed, in which the two-dimensional manifold is considered as the computation unit.


1.
Introduction. Incompressible flow takes place in many real situations of industrial and natural systems, such as pipe flow, flow around airfoils, atmospheric flow, blood flow, and convective heat transfer inside industrial furnaces. Typically, such flows are governed by the incompressible Navier-Stokes (NS) equations. Although significant progress has been achieved in numerical methods for solving incompressible NS equations in the last decades [28,29], as typically encountered in solving many nonlinear problems, these numerical methods often fail. Several primary reasons have been recognized in the literature and are reviewed briefly as follows.
First, it is well-known that for the convection-dominated problem, the NS equations often present serious numerical difficulties. At a high Reynolds number, the transient NS equations possess some hyperbolic property, which render standard finite element and finite difference methods usually exhibiting some combination of nonphysical oscillation and excessive numerical dispersion [17,5]. Many numerical methods have been developed to overcome such difficulties. For examples, the characteristics methods have demonstrated their efficiency for diffusion problems [11,25,1,8]. The idea in these methods is to rewrite the governing equations in terms of Lagrangian coordinates defined by the particle trajectories associated with the problem under consideration. Then it consists of combining a Galerkin finite element procedure with a discretization of the Lagrangian material derivative along the characteristics. The Lagrangian treatment in these methods greatly reduces the time truncation error as otherwise shown in Eulerian methods. The reason is that the approximation of ∂u/∂t by the standard backward-difference leads to errors in the form of C ∂ 2 u/∂ 2 t ∆t in suitable norms, whereas the characteristics method gives errors in the form of C ∂ 2 u/∂ 2 τ ∆t. For problems with significant convection, the solution changes much less rapidly in the characteristic τ direction than in the t direction [11]. Simultaneously, the characteristics method has been shown to possess remarkable stability [22]. In addition, it combines the advantages of the characteristics method, both being stable and explicit in treating the convection, in which the linear systems to be solved at each time step involve only the diffusion. And the second half-step is a saddle point problem where the pressure is a Lagrange multiplier associated with the incompressibility constraint enforced on the velocity.
Second, another difficulty for the numerical simulation of incompressible flows is that velocity and pressure are coupled by the incompressibility constraint, which is the same as in the second half-step of characteristics method. The interest in using projection method to overcome this difficulty started in the late 1960s in the ground breaking work of Chorin and Temam [9,10,26]. The most attractive feature of projection method is that, at each time step, only a sequence of decoupled elliptic equations need to be solved, which makes it very efficient for large-scale numerical computation and is very suitable for the dimension splitting method proposed in this article. Roughly speaking, the projection method is much more efficient than solving a coupled equations for velocity and pressure which would arise from the characteristic scheme of the NS equations. Moreover, an improved projection method [7] will be adopted in this paper, which is more accurate in dealing with fractional step equations and boundary conditions when compared with the traditional projection method.
Third, the high-cost of mesh generation and computation is a very significant problem in three-dimensional numerical analysis. In the 1950s, a very popular method was originally proposed by Wu [27] for the three-dimensional (3D) NS equations in turbo-machines. It suggested that the 3D solutions can be approximated through computing the flow on two classes of stream surface iteratively, which are referred as blade-to-blade (S1) stream surface and hub-to-tip (S2) stream surface. However, the algorithm is still complex because generation of two classes of stream surface is a difficult issue. Therefore, significant emphasis has been placed on developing simplified approach to formulation of stream surfaces. Especially, K. Li and coworkers introduced a dimension splitting method (DSM) [21] for threedimensional flow. This method is very promising and can be regarded as a renovation to Wu's S1/S2 approach. As the same as in Wu's S1/S2 theory, 3D flow computation can be translated into a series of two-dimensional computations in DSM, which avoids the need to generate 3D meshes. It is noted, nonetheless, unlike Wu's theory, the DSM only introduces one kind of two-dimensional manifold and does not make any assumption about the NS equations. Since the same kind of twodimensional manifold is used in the DSM [21], it is further feasible to construct a coarse-grained parallel algorithm. In our recent work [6], the finite difference is used to discretize the splitting direction and the finite element is used to approximate the equations on the two-dimensional manifold, which is a hybrid scheme, simplifies the designing procedure of original DSM and referred as the dimension-splitting method with difference form (or DSM-D).
The primary aim of this paper is to combine the DSM with the characteristic projection method for simulating the incompressible NS equations in three dimensions. Furthermore, in this newly proposed DSM, the finite element discretization is used in both splitting directions and the two-dimensional manifold. Hereinafter, it is termed as the DSM with characteristic projection (or DSM-C). Compared with DSM-D, mass conservation and nonconforming-pressure solution among manifolds are no longer needed for special treatment. In addition, higher calculation precision is expected. In this paper, first of all, temporal discretization is established along the characteristics according to the Lagrangian material derivative for the governing equations, and semi-discrete momentum equations are Stokes problem and then the projection method is used to deal with the incompressibility constraint. In conclusion, only independent linear elliptic equations are needed to calculate at each time step in the characteristic projection method for incompressible flow, then the DSM-C is constructed to compute these elliptic equations respectively. Because of the splitting property in DSM-C, all computations are carried out on these twodimensional manifolds. Namely, only a series of two-dimensional linear problems need to be solved. Therefore, it is asserted that the proposed algorithm is very effective for large-scale numerical computation.
2. Mathematical formulation. The governing equations of incompressible flow are the conservation of mass and momentum. The density-invariant incompressible model is adopted in this study. Then the dimensionless form of the governing equations defined in a three-dimensional domain Ω are expressed as follows: Continuity equation: Momentum equation: where u = (u, v, w) T is the dimensionless velocity, p the dimensionless pressure. Because of characteristics method the time term is considered. In the above equations the following dimensionalization was employed, These variables with the superscript are dimensional. U ∞ is a suitable characteristic velocity, L is a suitable characteristic length. Last, we consider the same initial-boundary value problem as in [7].
3.1. Characteristic temporal discretization. In this section we present a characteristics scheme for temporal semi-discretization of incompressible flow. First, the material derivative of a physical quantity is introduced, then we propose an Euler approximation of the characteristics curves. For a field φ we denote the material derivative byφ, x(c,t) =x. ( It is the characteristic curve that is simply the trajectory of the motion associated with the velocity field u. Hence, the material point c could be labeled as (x,t) according to the above equation. The position x(c, t) and physical quantity φ(c, t) are rewritten as x(x,t; t) and φ(x,t; t); the material derivative is also rewritten as, For the time variable, we introduce the number of time steps N and the time step t = t f /N , obtaining the uniform mesh of (0, t f ).
Moreover let x n ( x) be defined by which indicates the position at time t n of the material point through the space-time point ( x, t n+1 ). Settingx = x,t = t n+1 and applying Eq.(4) we obtain, Namely, ∂φ ∂t + u · ∇φ at space-time point (x, t n+1 ) is the material derivative at time t n+1 of φ of the material point through (x, t n+1 ).
In order to discretize the material time derivative in Eq.(5) we propose the Euler backward approximation, Then the two terms in the right-hand of the above equation is expressed in Euler point of view, So the temporal discretization of material derivative is further obtained, If φ = u, according to Eq.(4) the right hands of momentum equation (2) is written in the form of material derivative, that is, φ(x, t n ) is approximated by φ n and applying Eq.(6) and Eqs.(7) the characteristic temporal discretization of momentum equation is as follows, where u( x n ( x), t n ) will be given from the numerical approximations of x n ( x) in the following. From the above equation, there is a principal advantage in the characteristic temporal discretization of incompressible flow. First, the material derivative includes the nonlinear convective term; therefore, the nonlinear problem of NS equations has been dealt after the characteristic temporal discretization is applied. In this way, Eqs.(8) are linear elliptic equations. Simultaneously, an effective upwind method is also obtained for the convection dominated problem.
For the numerical approximations of Obviously, x n ( x) =x(t n ) and it is the position at time t n along the characteristic curve through ( x, t n+1 ). If u n is given, the initial-value problem could be approximated by Euler explicit difference scheme, When x n ( x) is solved, it is easy to compute u( x n ( x), t n ) in the finite element approximation. But it should be noted that x n ( x) is possibly out of the domain Ω, where u( x n ( x), t n ) could not be directly calculated. However t is very small, so there is little difference between x n ( x) and x. In the way u( x n ( x), t n ) could be estimated by the Taylor series,

Projection method.
It is emphasized again that through characteristic temporal discretization Eqs. (1) and (2) are linearized as follows: where Eqs. (11) are the Stokes problem for velocity u n+1 and pressure p n+1 , which are interconnected by the divergence constraint. Therefore, the projection method can be used to approximate the solution of the coupled system by first solving an analog to Eq.(8)(without regard to the divergence constraint) for an intermediate This is termed as the evolution step, in which ∇p n+1 is replaced by ∇p n , so it is the equation only for the intermediate quantity u * . And its boundary condition would be discussed later. Then we introduce well-known Helmholtz decomposition theorem [20]; it plays a key role in designing numerical procedure of projection methods. From Helmholtz decomposition theorem the intermediate velocity u * could be uniquely decomposed in the following form, In above decomposition the divergence-free is just the u n+1 that has been validated [7]. Then the decomposition (13) is substituted into the evolution step (12), Applying Eq.(11)−Eq.(14) the formula for pressure update is obtained, The last term t Re ∇(∆ψ n+1 ) of this formula plays an important role in computing the correct pressure gradient, which is different from Chorin's projection method and is one of approaches to overcome the numerical boundary layer of pressure [7]. Due to ∇ · u * = t∆ψ n+1 , the formula for the pressure update is further rewritten as, Compared with the first formula for pressure update the second derivative is not needed to calculate in Eq. (15). And the formula for velocity update is formed according to decomposition (13), Obviously, in Eq. (15) and Eq.(16) there is ψ n+1 , which is obtained by divergence operation on the two sides of Eq.(13), This is termed as the projection step. We would refer to methods of this type generically as incremental pressure-correction projection methods [15], since the projection step serves to compute an incremental-pressure gradient correction.
The equation for u * in the evolution step has been presented, now we discuss the corresponding boundary condition. In the traditional projection method it is always assumed that u * may not differ significantly from the fluid velocity because ∇p n is considered as a good approximation to ∇p n+1 . And thus a reasonable choice for the boundary condition may be u * | Γ = u D in the traditional method, such as Chorin's projection method, for ease of DSM-C algorithm design we also employ the boundary condition in the paper. Up to now the characteristic projection method for three-dimensional incompressible flow has been constructed and then we would describe the complete algorithm. If u n , p n are given, u n+1 , p n+1 could be solved by means of the steps (12), (17), (15) and (16). The characteristic projection method is of first order accuracy in time only, but it suffices for the steady problem considered in this article and high computational efficiency is also obtained.
From the above steps, the equations for u, p of incompressible flow have been completely decoupled. Furthermore the only independent linear elliptic equations are needed to compute at each time step. But it should be noted that there are no initial conditions of p, so at n = 0 a reasonable choice [15] for p 0 is like that the divergence operation is carried out on the two sides of and then the equation for p 0 is obtained. So far we present the semi-discrete approximation for three-dimensional incompressible flow, and subsequently dimension splitting method would be considered for the remaining three-dimensional elliptic equations.
where d = Z/K, K is the given integer, I = {0, 1, 2, · · · , K} is the index set and the function space on two-dimensional manifold is proposed, It is worth noting that the subscript D indicates not the two-dimensional manifold but the Dirichlet boundary condition in Γ. Then let T h = {κ} be a triangulation of two-dimensional manifold D, the P2-P1 finite element space on D is presented, Likewise X d = span{φ k , k ∈ I, φ k+ 1 2 , k ∈ {0, 1, · · · , K − 1}} would be P2 finite element space on T d and basis function φ k (z) at the integral node is defined by, basis function φ k+ 1 2 at the fractional node is described by, Therefore is the three-dimensional P2 finite element space, that serves as the base of DSM-C for adopting finite element discretization in both splitting direction and two-dimensional manifold simultaneously.

4.2.
Dimension splitting method. In characteristic projection algorithm only independent three-dimensional elliptic equations are needed to solve, therefore we introduce the dimension splitting method for three-dimensional Poisson equation.
The following Poisson equation with homogeneous Dirichlet boundary condition is taken as an example, which is defined in the three-dimensional domain, where Ω = D × (0, Z) is a cylindrical bounded domain in R 3 and D is a bounded domain in R 2 , shown in Figure 1.
Firstly introducing the following notations, According to the three-dimensional P 1 finite element space the approximation to u could be expressed by, The Laplace operator ∆ is divided into two parts, including of splitting direction and two-dimensional manifold, then because of the splitting property in M hd , the variational form about z coordinate direction could be derived firstly, ∀ϕ ∈ M d and ϕ| z=0,z=Z = 0, And then the variational form on two-dimensional manifold D is considered, ∀q ∈ Q h,0 , exchanging the integration order of first item in above equation, Therefore, the variational problem with dimension splitting property is obtained, According to Eq.(23) u i,h are the unknowns for the variational problem and are arguments on two-dimensional manifold D. Furthermore u i,h is only the function of (x, y) and could be regarded as a constant for z, so the scalar product (·, ·) z can be computed analytically and finally the above problem is just the variational problem on two-dimensional manifold D.
When the basis function in M d is used as test function ϕ, finite element discretization of Poisson equation with dimension splitting form is described as to find where, From Eq.(24) some two-dimensional variational problems are only needed to obtain the approximation u hd of three-dimensional Poisson equation, that is the basic idea of dimension splitting method. Taking (u hd , ϕ k ) z for example, it would be explained how to calculate the scalar product (·, ·) z and on account of compactly supported property of basis function ϕ k , In the same way the scalar product ( ∂u hd ∂z , ∂ϕ k ∂z ) z can be obtained. Although variational problems Eq.(24) are linear form, these unknowns u i,h on two-dimensional manifolds still couple. If these unknowns are solved simultaneously on two-dimensional manifolds, it would result in large scale computation of algebraic equation, even more important, it could not come true that two-dimensional problem on each manifold is computed independently, which is the main character of DSM. Therefore a natural way is that the argument on current manifold is regarded as the unknown and the arguments on other manifolds are regarded as the known. Then Jaccobi iterative scheme [23] is proposed, The finite element method could be employed to compute the variational problem on each manifold independently, hence a coarse-grained parallel algorithm would be constructed and the two-dimensional manifold is considered as the computation unit. At last dimension splitting method for Poisson equation is constructed and the same method could be applied to form the dimension splitting method for Eqs. (12) and (17). Since the derivation process is all the minor details and analogous to the one given above, the dimension splitting schemes for Eqs. (12) and (17) are directly presented in the following.
In the same way the following notations are firstly introduced, u k = u(x, y, z)| z=z k , p k = p(x, y, z)| z=z k , u k+ 1 2 = u(x, y, z)| z=z k+ 1 2 , ψ k = ψ(x, y, z)| z=z k Due to the three-dimensional finite element space in previous section, these approximations of unknowns in incompressible flow could be expressed by, The above formulas u h,k , p h,k , ψ h,k are finite element approximations to the unknowns u k , p k , ψ k on two-dimensional manifold. Then the finite element approximations to elliptic equations in characteristic projection algorithm are presented in the frame of DSM-C.

1) the evolution equation for
In Dirichlet boundary condition the intermediate velocity u * 0 , u * K on manifold D 0 and D K are known and test function should be picked as φ k , k = 1, · · · , K − 1. Hence there are 2K-1 two-dimensional problems needed to be computed, instead of computing original three-dimensional problem. Similarly the construction principle of scalar product (·, ·) z has been proposed in Eq. (25).
2) the projection equation for ψ n+1 In order to keep the uniqueness of the solution, a penalty term αψ n+1 hd is considered in the above variational problem. And in Neumann boundary condition the intermediate arguments ψ 0 , ψ K on manifold D 0 and D K are unknown and test function should be picked as φ k , k = 0, 2, · · · , K. For ψ n+1 hd there are K+1 two-dimensional problems needed to be computed, instead of computing original three-dimensional problem. And for the velocity and pressure update equation (15), (16) and pressure initialization equation, the analogous method could be adopted to obtain relevant formulation.
Although only two-dimensional problems are computed on manifold, these unknowns on two-dimensional manifolds are also coupled. For the purpose of decoupling the computation between manifolds and constructing a coarse-grained parallel algorithm, a natural idea is that the argument on current manifold is still regarded as the unknown, while the arguments on other manifolds are regarded as the known, that is just Jaccobi iteration. Nowadays two-dimensional problem on manifold could be computed independently, thus the parallelization strategy of DSM-C could be constructed easily. And if np processors are used, the index of manifold computed by the processor p is, when mod(K − 1, np) = 0, The idea of above strategy is that amount of calculation is equally loaded on each processor, when the number of manifold K-1 is divisible by the number of processor np. When it is not divisible, the divisible part is equally loaded on each processor and the rest is loaded on processor orderly.

5.
Numerical examples. In order to assess the main characteristics of the proposed DSM-C and to test its effectiveness in the context of high-performance computing, we develop two numerical examples herein in this paper. The first is a classical problem with analytical solution; the second is the well-known 3D lid-driven cavity flow problem. The numerical implementation was programmed through using OpenMP (Open specifications for Multi-Processing) for parallel computing [16].

Analytical solution problem.
It is well known that the analytical solution to the incompressible flow equation exists in a few special cases. In order to avoid unnecessary influence, the source term f is constructed specifically to lead the known analytical solution problem. Particularly, the same steady problem as used in Ref [6] is adopted, for which the prescribed solution is given in the domain of u(x, y, z) = (u 1 (x, y, z), u 2 (x, y, z), u 3 (x, y, z)), p(x, y, z) = 3 − (x 3 + y 3 + z 3 ), u 1 (x, y, z) = sin(πy(y − 1)) sin(πz(z − 1)), A simple calculation could verify that ∇ · u = 0. The source term f is obtained such that (u 1 , u 2 , u 3 , p) is the solution of the Navier-Stokes equations (1) and (2), then the velocity boundary condition is set equal to the exact solution. The uniform mesh partition of the two-dimensional manifold D into triangular elements is obtained by dividing D into sub-squares of equal size and then drawing the diagonal in each sub-square (Figure 7), where h is the side of sub-square and partition interval d = h is used in the z direction. In order to verify the performance of the proposed method quantitatively, we define some parameters, where κ div is used to denote that numerical solution meets the continuity equation quantitatively. Rigorous error estimation for the DSM-C is not in the scope of this paper. Instead, it is obtained through the known analytical solution example. The L 2 norms of the error of the numerical solution u hd , p hd are defined as, where u, p are the analytical solutions, u hd , p hd are the numerical solutions. Let the error estimate for velocity be, where C and α are constants independent of the partition. The rate of convergence α could be estimated as follows. Based on the two different partitions h 1 and h 2 . we have, If we only consider the case of equality in the formulations above, the logarithm of quotient obtained by dividing can be treated as the convergence rate.
log(h 1 /h 2 ) As such, the convergence rate of pressure is obtained similarly. In order to minimize the impact of temporal discretization scheme and estimate more accurate numerical error of spatial discretization scheme, the size of time step is set to be ∆t = 0.0001. The resulting L 2 norms of numerical solution, κ div , and the rate of convergence α are summarized in Table 1. From κ div in the table, the incompressibility constraint problem is well solved by the projection method. Moreover, there is no special treatment in the DSM-C and it is more convenient than the DSM-D [6]. We further compare the DSM-C with the DSM-D results, and Table 2 presents the convergence rates of numerical solutions and CPU time for using the two methods. It is obviously observed that the accuracy of DSM-C with only finite element scheme is better than DSM-D with hybrid scheme. It is more important to note that the characteristic projection method is acceptable for dealing with the convective term and for decoupling the velocity and pressure in DSM-C. Therefore, the CPU cost is much less than DSM-D, in which the coupled computation of velocity and pressure is adopted. However, the algorithm construction of DSM-C is more complicated; and the CPU time in the generation of stiffness matrix is more than DSM-D, although the total CPU time of DSM-C and DSM-D is relatively close on a coarse mesh. It is further observed that the superiority of decoupled calculation becomes more outstanding with mesh refinement; this is because that solving linear equations takes up most of CPU time on fine mesh.
Last, the parallel performance of DSM-C is tested (31) and (32). In the following, the parallel speedup ratio S p and the efficiency E p are used to analyze the parallel performance, where T 1 is the CPU time of serial computation, T p the CPU time of parallel computation using p processors. In order to make fair comparison, 200 time steps and one Jaccobi iteration are applied. Table 3 and 4 tabulate the performance analysis results from using the parallel algorithm in DSM-C with various partitions h.
The results show that the parallel efficiency is lower in the coarse mesh. This is expected because the proportion of parallel computation relative to the total  computation is less. When the proportion rises as the mesh is refined, the parallel efficiency is observed to rise. Scalability is another important issue in parallel computing. The issue becomes significant when solving large-scale problems with many processors, because the scalability represents the parallel efficiency E p , which should not decrease greatly as p increases. For DSM-C, as p increases, the parallel efficiency E p is expected to decrease when the number of processors is greater than the number of two-dimensional manifolds, because the extra processors do not save CPU time and increase the communication time among threads. Hence the scalability would become stronger as the mesh scale h gets smaller. On the other hand, when compared with the DSM-D, the parallel performance of DSM-C is a little higher at coarse mesh. This is because that DSM-C is more complicated than DSM-D and the extra computation could be applied on two-dimensional manifold in parallel, which is the proportion of serial computation amount in DSM-C is smaller. Lastly, it is noted that the parallel performance of these two methods approach to be closer as the mesh is more refined.

5.2.
Three-dimensional lid-driven cavity problem. Herein we start with numerically solving the classical lid-driven flow problem. This test has been widely used as a benchmark for numerical methods and has been analyzed by a number of authors [14,3]. In the following, it is used to examine the computational capability of DSM-C for simulating the real flow. The geometric configuration and the prescribed boundary conditions are depicted in Figure 8. The geometry describes a cubic domain with the upper wall in linear motion. The fluid is assumed to be fully enclosed with no free surfaces, and gravity is neglected. Therefore, setting f = 0 forms the incompressible flow in Eqs. (2). The fluid motion is generated by the top lid that moves in the x-direction with a constant velocity U ∞ = 1, in fact in order to avoid discontinuity in the Dirichlet boundary conditions, we resort to the following velocity profile as reported in [4], The cavity height L = 1 is used as the characteristic length, and the velocity of the moving wall U ∞ as the characteristic velocity. The viscosity is adjusted in order to obtain different Reynolds number of 100, 400 and 1000. The two-dimensional mesh made of 1 h × 1 h = 64 × 64 elements, the z-direction partition d = h and the size of time step ∆t = 0.01 have been used in the calculations.  The experimental visualization of three-dimensional lid-driven cavity flow was originally accomplished at Stanford University by Koseff et al. [19,18], and the numerical prediction of the flow shall be credited to the work of Freitas et al. [12,13]. Their results of experimental and numerical prediction show that the steady state does not exist but the transient behavior becomes periodic for Re ≥ 3200 and the steady problem is considered in the example. Therefore, the Reynolds number is not chosen too high herein.
Due to the main interest in comparing the performance of the implemented method described in the previous sections and the behavior of the reported solution. Therefore, a set of numerical experiments has been performed. Figures 4,5 and 6 plot the velocity profiles for u x and u y along y = 0.5 and x = 0.5 for different Reynolds numbers on the middle plane z=0.5. For comparison, the results of Babu&Korpela [2] and Shu [24] are also included. First of all, one may clearly see the agreement of the solutions. Figures 7 and 8 show that streamline profiles on the three middle planes x=0.5, y=0.5 and z=0.5 for various Reynolds numbers, which depicts that Taylor-Görter-Like (TGL) [19] vortices change with the Reynolds number. When the Reynolds number is low, the primary eddy on the middle planes z=0.5 locates at the upper center. With the rising of the Reynolds number the primary eddy approaches to the center and the corner eddies begin to appear. More flow structure is illustrated by 3D streamline displayed in Figures 9 and 10. It is noticed that in the steady case, the particles are limited to a single side of the cavity, which make the flow fairly symmetric. All qualitative and quantitative results demonstrate that DSM with characteristic projection method is effective for computing the real flow. Conclusions. In this paper we propose a dimension splitting method for the three-dimensional incompressible flow, in which the characteristic projection method is employed to deal with the nonlinear convective term and the incompressibility constraint such that velocity and pressure could be computed independently. In using this novel dimension splitting method, termed the DSM-C, the threedimensional solution can be approximated by the solutions on two-dimensional manifolds, which is the most significant feature of this new method. Therefore, there is no need to generate three-dimensional grids, which greatly reduces the difficulty and computational cost in mesh generation. A coarse-grained parallel algorithm can be constructed, where the two-dimensional manifold is considered as the computation unit. Compared with the dimension splitting method with difference form (DSM-D) [6], the DSM-C presents two main improvements. Velocity and pressure have been decoupled and only independent linear elliptic equations are needed to calculate at each time step, while the velocity and pressure are solved simultaneously for DSM-D. Thus the computational cost of DSM-C is much less than DSM-D, which is validated in the numerical example. On the other hand, the finite element discretization is used in both splitting direction and two-dimensional manifold for DSM-C; whereas for DSM-D, the finite difference is used to discretize the splitting direction and the finite element is used to approximate the equations on the two-dimensional manifold, leading to a hybrid scheme. Additionally, incompressibility constraint among manifolds is no longer needed for special treatment and higher calculation precision is achieved in DSM-C, which is also validated in the numerical examples. A comprehensive comparison between the present results and ones appearing in the literature is reported in the numerical examples, and the data obtained by DSM-C have relatively good agreement with the given reference results. Therefore, the dimension splitting method in this paper is a highly efficient algorithm for simulating three-dimensional incompressible flow.