Construction of Solenoidal Immersed Velocity Vectors Using the Kinematic Velocity–Vorticity Relation
Abstract
The present paper suggests a method for obtaining incompressible solenoidal velocity vectors that satisfy approximately the desired immersed velocity boundary conditions. The method employs merely the mutual kinematic relations between the velocity and vorticity fields (i.e, the curl and Laplacian operators). An initial nonsolenoidal velocity field is extended to a regular domain via a zerovelocity margin, where an extended vorticity is found. Recalculation of the velocities (subjected to appropriate boundary conditions), yields the desired solenoidal velocity vector. The method is applied to the two and threedimensional problems for the homogeneous Dirichlet, as well as periodic boundary conditions. The results show that the solenoidality is satisfied up to the machine accuracy for the periodic boundary conditions (employing the Fourier–spectral solution method), while an improvement in the solenoidality is achievable for the homogenous boundary conditions.
keywords:
solenoidal velocity vectors; immersed boundary conditions; velocity–vorticity relation; harmonic functions; periodic and homogenous boundary conditions1 Introduction
The NavierStokes Equations (NSE) are derived for a continuum. Therefore, many difficulties arise in dealing with any kind of discontinuity, including finitesize solution domains (simplyconnected or multiplyconnected). Particularly in an incompressible flow, the problem usually appears as nonsolenoidal velocity fields. Consequently, satisfying the mass conservation has been one of the most challenging issues of almost all algorithms of incompressible flow simulations, whether in the primitive variables formulations or the other formalisms. As it is revealed, the boundary conditions have a vital role here; and therefore, appropriate definition of the boundary conditions received significant attention in the literature [9, 12, 13].
The literature of the subject is too rich to completely be summarized in this paper (see e.g. [13] and the references in there); however, as the most pertinent works, we shall make reference to the works of Quartapelle and ValzGris [11], and also Quartapelle [12], which studied the issue extensively. Calhun [3] used the immersed boundary technique in the vorticity–stream function formulation, and introduced an appropriate distribution of the vorticity in order to satisfy the overall mass balance. From another viewpoint, the problem was studied in relation to the compatibility conditions (also called the space–time corner singularity). The required compatibility conditions for the NSE was studied by Temam [17], Gresho and Sani [9], and also by Boyd and Flyer [1], and then by Flyer and Swarztrauber [8]).
The main contribution of the present paper is to suggest a method for construction of a nearly solenoidal velocity field, which satisfies approximately the desired immersed velocity boundary conditions, using the kinematic velocity–vorticity relation. The central idea is to use the wellknown property that the divergence of the velocities, obtained from the velocity–vorticity equation is a harmonic function; and therefore, gets its maximum values on the boundaries. Consequently, the divergence can be reduced (or even vanished) all over the field, by restricting (or vanishing) the divergence near the boundaries. In the present method, these restrictions are implemented by considering a zerovelocity margin in the vicinity of the outer regular boundaries.
Embedding nonregular solution domains in a bigger regular domain was shown to be useful in solution of the scalar Poisson’s problems [2, 14]. Now the present work may be seen as an extension of the method to the vector Poisson’s problems that the solenoidal solutions are desired. As it will be shown, a solenoidal velocity field (within the machine accuracy) is achievable for the periodic boundary conditions (that is, the noboundary problems), and also improvements in the solenoidality is possible for the homogeneous Dirichlet boundary conditions. The method can be used in recently popularized onefluid [19], or onecontinuum [16] immersed boundary methods, as it was used successfully our previous work [15].
The paper is continued by description of the suggested algorithm, follows by a section contains our arguments about the reasons that the method works (that is, 3). As our numerical experiments, two and threedimensional fixed and moving boundary problems are examined, with homogeneous as well as periodic boundary conditions.
2 The suggested procedure
Given an arbitrary (may be nonsolenoidal) velocity vector , defined on the flow domain , with arbitrary boundary , where ; one can find the corresponding vorticity vector on the . Now, it is wellknown that solution of the Poisson equations
(1) 
on implies
(2) 
which means although is not necessarily zero, it is a harmonic function, and therefore, gets its maximum values on the boundaries. The above statement has been demonstrated several times by many authors [11, 12], and motivated many attempts in obtaining solenoidal velocities by vanishing the divergence right at the boundaries [3, 6].
Now, we suggest limiting the divergence inside the solution domain via controlling the divergence in the vicinity of the outer boundaries (not merely at the boundaries), via definition of a zerovelocity margin. According to Fig. 1, The regular solution domain contains all the fluid–solid system, including the moving regions of the fluid with given velocity boundary
conditions (that is ), in addition to the fixed or moving immersed solid objects. Now this fluid–solid system is embedded in a regular domain via a zerovelocity margin .
With the above definitions in mind, the following procedure is suggested:

Given the arbitrary (may be nonsolenoidal) velocity vector defined on , the extended velocity , defined on the regular domain , is found via addition of the zerovelocity margin , and other immersed regions ( and others) to the .

The vorticity is obtained on .

The modified extended velocity vector is obtained from solution of the Poisson’s equations
(3) subjected to some boundary conditions on , which will be discussed in the next section.

The desired modified velocity vector is obtained by discarding and the other immersed regions ( and others) .
The above procedure is illustrated in Fig. 2.
Remark 1
In fact for the multiply connected domains, the extremum of the divergence may occur on the boundaries of the holes (e.g. or others). However, in the present work we followed an immersed boundary approach, that is, includes the velocities of all the holes of . In this way, the extremum of the divergence occurs on , and therefore, it is just needed to control the divergence on .
3 On the Solenoidality of
We shall examine the solenoidality of . To this end, we consider two different boundary conditions for Eq. (3), that is, the periodic boundary conditions, which is suitable for the Fourier spectral solvers, and the homogenous Dirichlet boundary conditions, which can be implemented easily in the finitedifference or finitevolume solvers.
3.1 Periodic Boundary Conditions
In the case of periodic boundary conditions (that is, noboundary problems), equation (3) can be transformed into the Fourier space, which leads to
(4) 
in which ; and stands for the wavenumber vector, and . Obviously, is perpendicular to the wavenumbers vector , that is
(5) 
which is a translation of solenoidality of in the Fourier space [4]. Now the desired (solenoidal) velocity vector can be obtained simply by ignoring the solution in the and the other immersed regions. More than the above arguments, our numerical experiments have also confirmed that in this case is solenoidal within the machine accuracy. In fact, achieving perfect conservation of mass for the noboundary problems is not so surprising, because many assumptions that the NSE are based on are satisfied [7].
3.2 Homogenous Dirichlet boundary conditions
For the homogenous boundary conditions on , the modified extended velocity is not strictly solenoidal. However, since , the divergence of gets its extremum on the boundary . In the other words, the extremum of divergence occurs in , and therefore, discarding the solution in reduces the divergence such that
(6) 
In fact, the maximum values of are discarded by cutting the solution in ; and therefore, the divergence of is less than the divergence of .
As it will be seen in our numerical experiments, finitedifference solution of Eqns. (3) with homogeneous boundary conditions, modifies the solenoidality of an initially nonsolenoidal velocity field.
4 Numerical Experiments
In this section, the proposed algorithm is assessed via examination of some fixed, as well as moving boundary problems in two and three dimensions. We implemented the method both in the Fourierspectral, and the finitedifference solutions.
4.1 The modified Chandrashkar–Reid flow
In the study of the orthogonal basis functions that satisfy the noslip conditions, Chandrashkar and Reid proposed a velocity field that satisfies the noslip boundary condition in one direction, and is periodic in the other direction [5]. Then this flow was appeared in the work of Rempfer [13]. In the light of their works, it is convenient here to construct a velocity field, which satisfies the noslip conditions in both directions. Tho this end, we first introduce a onedimensional characteristic function
(7) 
where , and are chosen such that , within the machine accuracy. Using this characteristic function, we construct the desired streamfunction
(8) 
in which the normalization factor is chosen such that the unit maximum velocity occurs in the field, that is
(9) 
Now the velocities are directly obtainable from this streamfunction as
(10)  
(11) 
The left panel of Fig. 3 shows the above velocity filed on a point grid, which formally satisfies the noslip conditions in both directions. In this configuration we defied .
It is worth mentioning that although the above velocity vector is obtained from a streamfunction, it is not solenoidal. In fact, the divergence can be written as
(12) 
which obviously is not necessarily zero on the boundaries. For example, on the left boundary , while .
This fundamental problem is one of the difficulties that the formulation is suffered from. In fact, in solution of , two sets of boundary conditions are needed to be satisfied (that is, , and ), which makes the problem overdetermined [13]. The divergence of the flow filed is shown in Figure. 4. As one can see, the divergences are not zero, especially right at the boundaries.
Now the method is applied to this flow. To this end, the solution domain is extended to a bigger regular domain , and then the corresponding vorticity is found on . The extended velocity field is shown in the right panel of Fig. 3. Now Eqns. (3) are solved for the periodic as well as homogenous boundary conditions using the Fourier spectral, and the finitedifference methods; and the desired velocity is extracted from . The distributions of for both cases are shown in Fig. 5. As one can see, for the Fourierspectral solutions the divergences are reduced to the machine accuracies (see the left panne of Fig. 5). On the other hand, for the homogenous boundary conditions (which Eqns. (3) are solved using the finitedifference method), although the divergences are no longer vanished completely, however, reduced appreciably. Table LABEL:tab1 summarizes the maximum values of the divergences in the above cases.
It is worth mentioning that although the above procedure improves the solenoidality of the flow, it may change the velocity boundary conditions. In fact, the velocities are changed all over the field, especially near the boundaries. This issue will be discussed in detail in the next section, because for the present test case boundary condition changes have not been sensible.
4.2 Moving square cylinder
In this section, the capability of the method in dealing with multiplyconnected domains is examined. In this context, the practically important (and academically challenging) problem of a moving obstacle in a quiescent fluid is considered. Since the objective of the present paper is not implementation of the nonregular boundaries, a square cylinder, coinciding with the Cartesian grid, is considered. For the problem of nonregular immersed boundaries one can see the related references (see e.g. [14, 15], or [21]).
Consider a regular fluid–solid domain , consists of a square cylinder , and the flow domain . Moreover, assume that it is desired to implement the velocity boundary condition to the square cylinder (which will be treated here as an immersed boundary condition). Implementation of such a boundary condition imposes some discontinuities to the solution for , the problem which is known as the corner singularity, or nonconsistency of the boundary and initial conditions and were studied previously (e.g. see [17, 1, 18, 8]). The initial velocity and the distribution of the divergences are shown in Fig. 6, calculated on a – point grid. As one can see, the continuity is lost, especially on the boundaries of the moving obstacle, because of implementation of the sharp immersed velocity boundary conditions.
Now the method is applied to this flow for the homogeneous and periodic boundary conditions. It should be noted that since the velocities are already zero in the vicinity of the regular boundary , the zero velocity margin is not needed indeed. First we consider the homogeneous boundary conditions. In this case, the vorticity is found, and the Poisson’s equations (3) are discretized on using the finitedifference method. The homogeneous Dirichlet boundary conditions are implemented to the discretized equations, and the resulting linear problem is solved using the point successive over relaxation (PSOR) method. The modified velocity field and the divergences are shown in Fig. 7. As one can see in the left panel, the initial immersed nonsolenoidal velocities are now spread across the domain , such that the nonzero velocities are observable near the homogeneous outer boundaries (and therefore, the divergences are reduced near the square).
As a consequence, the divergences in the right panel show an appreciable reduction in the maximum values (with the factor of in comparison to Fig. 6), and furthermore, these maximum values are now occurred on the outer (regular) boundaries.
For the periodic boundary conditions, since the initial velocities are vanishing in the vicinity of the regular boundary , the Fourier spectral solver may be used easily [2, 14, 15]. The vorticity (on ) is obtained in the Fourier space, and the modified velocities are found via solution of Eqns. (3). The modified velocities and the divergences are shown in Fig. 8. As one can see, solenoidality of the velocity field is achieved up to the machine accuracy. With regard to the modified velocities (the left panel of Fig. 8), we should emphasized that since the periodic boundary conditions (not homogenous ones) are implemented, the modified velocities are not vanished on . In fact, nonzero velocities are induced on the by the moving obstacle. However, note that the overall mass balance is satisfied, since the divergences are negligibly small.
On the other hand, it should be emphasized that satisfying the solenoidality changes the immersed velocity boundary conditions. Similar problem was observed in many other immersed boundary methods for the incompressible flow, and motivated emergence of some methods like the multidirect forcing [22, 23], or the implicit forcing method [10, 20] in the last years. To elucidate the issue, the modified velocities on the edge of the square cylinder are compared with the desired boundary conditions in Fig. 9. As one can see, both components of the velocity (i.e. and ) are changed on all the immersed boundaries. In fact this is the cost that we have to pay for working with physical (solenoidal) velocities. However, a body of evidence show that these discrepancies eliminates by developing the solution in the next times [15, 22, 23]; and furthermore, there are some remedies for overcoming this difficulty [15, 23, 20].
4.3 Extension to the threedimensional problems
Although merely the twodimensional problems are considered up to now, the method may be extended easily to the threedimensional flow. However, in these problems all three components of the vorticity vector present, and consequently we are faced with three Poisson’s problems. Therefore, the method seems to be not so efficient (in comparison to e.g. some primitive variables formulations of the NSE). Nevertheless, since the boundary conditions, and the right hand sides of three Poisson’s problems (3) are uncoupled, the problems are completely parallelizable. Therefore, the method has the potential of being competing, at least in principle.
In the present work, our numerical experiment is restricted to the periodic boundary conditions (for which we can use the efficient Fourier–spectral solvers). The problem set up (which is an extension to the last moving square cylinder problem), is shown in Fig. 10. A moving cube with velocity is placed in a regular solution domain , occupied by a quiescent fluid. The extended vorticity vector is obtained and the modified velocities are recalculated on a uniform grid, using the Fourier–spectral method. The results are shown in Fig. 10. As one can see, the cube velocity is induced to the whole of the solution domain, and the modified velocity vector is perfectly solenoidal.
4.4 Conclusions
A method was suggested for construction of (nearly) solenoidal velocity vectors which satisfy (approximately) the desired immersed velocity boundary conditions. The method consists of recalculation of the velocities from an extended vorticity field, obtained from an extended initial nonsolenoidal velocity vector. It was shown that the solenoidality of the modified velocities are depended on the regular boundary conditions, and two different boundary conditions were considered. For the homogeneous boundary conditions, the method improves the solenoidality, and for the periodic boundary conditions (that is, the noboundary problems), the solenoidality satisfies, up to the machine accuracy. The modified solenoidal velocities are nolonger satisfy the desired immersed velocity boundary conditions, however (as it is shown in other references), the discrepancies are eliminated by developing the solution. The method was applied to the two and threedimensional problems, and effects of the method on the velocity fields were discussed.
References
 [1] J.P. Boyd, N. Flyer, Compatibility conditions for timedependent partial differential equations and the rate of convergence of Chebyshev and Fourier spectral methods, Comput. Methods Appl. Mech. Engrg. 175 (1999) 281309.
 [2] J. P. Boyd, A Comparison of Numerical Algorithms for Fourier Extension of the First, Second, and Third Kinds, J. Comput. Phys. 178 (2002) 118160 .
 [3] D. Calhoun, A Cartesian grid method for solving the streamfunction–vorticity equations in irregular geometries, Ph.D. thesis, University of Washington, The USA, 1999.
 [4] C. Canuto, M. Y. Hussaini, A. Quarteroni and T.A Zang, Spectral Methods: Fundamentals in Single Domains. SpringerVerlag, Berlin Heidelberg, 2006.
 [5] S. Chandrasekhar, W. H. Reid, On the expansion of functions which satisfy four boundary conditions, Proc. Nat. Acad. Sci. 43 (1957) 521527.
 [6] W. E, J. Liu, Vorticity Boundary Condition and Related Issues for Finite Difference Schemes, J. Comput. Phys. 124 (1996) 368382.
 [7] C. Foias, O. Manely, R. Rosa, R. Temam, Navier–Stokes Equations and Turbulence, Cambridge University Press, 2001.
 [8] N. Flyer, P. N. Swarztrauber, The Convergence of Spectral and Finite Difference Methods for InitialBoundary Value Problems, SIAM J. Sci. Comput., 23 (5) (2002) 1731–1751.
 [9] P.M. Gresho, R.L. Sani, On pressure boundaryconditions for the incompressible Navier–Stokes equations, Int. J. Numer. Methods Fluids 7 (10) (1987) 11111145.
 [10] D.V. Le, B.C. Khoo, K.M. Lim, An implicitforcing immersed boundary method for simulating viscous flows in irregular domains, Comput. Methods Appl. Mech. Engrg. 197 (2008) 21192130.
 [11] L. Quartapelle, Vorticity Conditioning in the Computation of TwoDimensional Viscous Flows, J. Comput. Phys., 40 (1981) 453477.
 [12] L. Quartapelle, Numerical solution of the incompressible Navierstokes equations. BirkhauerVerlag, Basel, 1993.
 [13] D. Rempfer, On boundary conditions for incompressible NavierStokes problems. Appl. Mech. Rev. 59 (2006) 107125.
 [14] F. Sabetghadam, S. Sharafatmandjoor, F. Norouzi, Fourier spectral embedded boundary solution of the Poisson’s and Laplace equations with Dirichlet boundary conditions, J. Comput. Phys. 228 (1) (2009) 5574.
 [15] F. Sabetghadam, M. Badri, S. Sharafatmandjoor, H. Kor, An Immersed Boundary Fourier Pseudospectral Method for Simulation of Confined Twodimensional Incompressible Flow, Under revision (J. Comput. Phys.), also available fron arXiv (arXiv:1110.5984).
 [16] K. Sugiyama, S. Ii, S. Takeuchi, S. Takagi, Y. Matsumoto, A full Eulerian finite difference approach for solving fluid–structure coupling problems, J. Comptu. Phys. 230 (2011) 596627.
 [17] R. Temam, Behaviour at time t = 0 of the solutions of semilinear evolution equations, J. Differ. Equ. 17 (1982) 7392.
 [18] R. Temam, Suitable initial conditions, J. Comput. Phys. 218 (2006) 443450.
 [19] G. Tryggvason, M. Sussman, M.Y. Hussaini, Immersed boundary methods for fluid interfaces, in: A. Prosperetti, G. Tryggvason (Eds.), Computational Methods for Multiphase Flow, Cambridge University Press, 2007.
 [20] M. Uhlmann, An immersed boundary method with direct forcing for the simulation of particulate flows, J. Comput. Phys., 209 (2005) 448476.
 [21] M. Vanella, P. Rabenold, E. Balaras, A directforcing embeddedboundary method with adaptive mesh refinement for fluidstructure interaction problems, Journal of Computational Physics 229 (2010) 64276449.
 [22] K. Luo, Z. Wang, J. Fan, K. Cen, Fullscale solutions to particleladen flows: Multidirect forcing and immersed boundary method, Physical Review E 76 (2007) 066709(19).
 [23] Z. Wang, J. Fan and K. Cen, Immersed Boundary Method for the Simulation of 2D Viscous Flow Based on VorticityVelocity Formulations, J. Comput. Phys. 228 (2009) 15041520.