Construction of Solenoidal Immersed Velocity Vectors Using the Kinematic Velocity–Vorticity Relation

Construction of Solenoidal Immersed Velocity Vectors Using the Kinematic Velocity–Vorticity Relation

Fereidoun Sabetghadam Shervin Sharafatmandjoor, Mehdi Badri Mechanical and Aerospace Eng. Dept., Science and Research Branch, Islamic Azad University (IAU), Tehran, Iran

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 non-solenoidal velocity field is extended to a regular domain via a zero-velocity margin, where an extended vorticity is found. Re-calculation of the velocities (subjected to appropriate boundary conditions), yields the desired solenoidal velocity vector. The method is applied to the two- and three-dimensional 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.

solenoidal velocity vectors; immersed boundary conditions; velocity–vorticity relation; harmonic functions; periodic and homogenous boundary conditions

1 Introduction

The Navier-Stokes Equations (NSE) are derived for a continuum. Therefore, many difficulties arise in dealing with any kind of discontinuity, including finite-size solution domains (simply-connected or multiply-connected). Particularly in an incompressible flow, the problem usually appears as non-solenoidal 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 Valz-Gris [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 well-known 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 zero-velocity margin in the vicinity of the outer regular boundaries.
Embedding non-regular 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 no-boundary problems), and also improvements in the solenoidality is possible for the homogeneous Dirichlet boundary conditions. The method can be used in recently popularized one-fluid [19], or one-continuum [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 three-dimensional fixed and moving boundary problems are examined, with homogeneous as well as periodic boundary conditions.

2 The suggested procedure

Given an arbitrary (may be non-solenoidal) velocity vector , defined on the flow domain , with arbitrary boundary , where ; one can find the corresponding vorticity vector on the . Now, it is well-known that solution of the Poisson equations


on implies


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

Figure 1: The flow domain , together with fixed or moving obstacle(s) , and other given-velocity regions , are embedded in the regular solution domain (with regular boundary ), via a zero-velocity margin .

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 zero-velocity 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 zero-velocity margin .
With the above definitions in mind, the following procedure is suggested:

Figure 2: Flowchart of the proposed algorithm. The initial velocities on the flow domain is extended to a bigger regular domain , where the extended vorticity is obtained and the Poisson’s equations are solved on. Then the modified velocities are extracted from .
  • Given the arbitrary (may be non-solenoidal) velocity vector defined on , the extended velocity , defined on the regular domain , is found via addition of the zero-velocity 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


    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 finite-difference or finite-volume solvers.

3.1 Periodic Boundary Conditions

In the case of periodic boundary conditions (that is, no-boundary problems), equation (3) can be transformed into the Fourier space, which leads to


in which ; and stands for the wavenumber vector, and . Obviously, is perpendicular to the wavenumbers vector , that is


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 no-boundary 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


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, finite-difference solution of Eqns. (3) with homogeneous boundary conditions, modifies the solenoidality of an initially non-solenoidal 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 Fourier-spectral, and the finite-difference solutions.

4.1 The modified Chandrashkar–Reid flow

In the study of the orthogonal basis functions that satisfy the no-slip conditions, Chandrashkar and Reid proposed a velocity field that satisfies the no-slip 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 no-slip conditions in both directions. Tho this end, we first introduce a one-dimensional characteristic function

Figure 3: Left: The modified Chandrashkar–Reid flow, defined on . The no-slip condition is satisfied formally on the boundaries . Right: The solution domain is extended to , via addition of the zero velocity margin .

where , and are chosen such that , within the machine accuracy. Using this characteristic function, we construct the desired streamfunction


in which the normalization factor is chosen such that the unit maximum velocity occurs in the field, that is


Now the velocities are directly obtainable from this streamfunction as

Figure 4: Distribution of the divergence for the initial modified Chandrashkar–Reid flow. Although the velocities are obtained from a streamfunction, the divergence is not zero and its maximum values are occurred at the boundaries.

The left panel of Fig. 3 shows the above velocity filed on a -point grid, which formally satisfies the no-slip 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


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 over-determined [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 finite-difference 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 Fourier-spectral 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 finite-difference 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.

Figure 5: The divergences for the extended modified Chandrashkar–Reid flow. Left: Fourier spectral solution for the periodic boundary conditions on . The divergence is reduced to the machine accuracy. Right: Finite-difference solution for the homogenous boundary conditions on (the dashed lines show the divergence of the extended velocity , and the solid lines show the desired final velocity ). The divergences are no longer within the machine accuracy, however, appreciable reduction is observable (c.f. Fig. 4).

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 multiply-connected 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 non-regular boundaries, a square cylinder, coinciding with the Cartesian grid, is considered. For the problem of non-regular immersed boundaries one can see the related references (see e.g. [14, 15], or [21]).

Figure 6: Left: Moving square cylinder in a quiescent fluid. The immersed velocity boundary conditions are implemented sharply on the square cylinder , placed in the flow domain . Right: The resulting divergences. As one can see, the continuity is lost because of implementation of the boundary conditions.

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 non-consistency 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 finite-difference 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 non-solenoidal velocities are now spread across the domain , such that the non-zero velocities are observable near the homogeneous outer boundaries (and therefore, the divergences are reduced near the square).

Figure 7: The modified velocities for the moving square cylinder with homogenous boundary conditions. Left: The induced velocities are observable almost on the whole of the regular domain . Right: The maximum values of the divergences are reduced appreciably (c.f. Fig. 6), and occur on the outer boundaries.

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.

Figure 8: The modified velocities for the moving square cylinder with periodic boundary conditions. Left: The periodic boundary conditions induced non-zero velocities on . Right: However, the velocity vector is solenoidal within the machine accuracies.

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, non-zero velocities are induced on the by the moving obstacle. However, note that the overall mass balance is satisfied, since the divergences are negligibly small.

Figure 9: Errors in immersed boundary conditions of the moving square cylinder due to imposing the solenoidality. The positions are presented with respect to the lower left corner of the square . Left: Errors in and on the left edge of the square cylinder. Right: Errors on the upper edge of the square cylinder.

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 multi-direct 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 three-dimensional problems

Although merely the two-dimensional problems are considered up to now, the method may be extended easily to the three-dimensional 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 un-coupled, the problems are completely parallelizable. Therefore, the method has the potential of being competing, at least in principle.

Figure 10: The modified flow field for a moving cube (with velocity ) in the quiescent fluid with periodic boundary conditions. Left: Streamlines and the contour plots of on half of the solution domain. Right: Iso-surfaces of level of the divergences of the flow. As one can see, the flow is divergence-free within the machine accuracy.

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 re-calculated 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 re-calculation of the velocities from an extended vorticity field, obtained from an extended initial non-solenoidal 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 no-boundary problems), the solenoidality satisfies, up to the machine accuracy. The modified solenoidal velocities are no-longer 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 three-dimensional problems, and effects of the method on the velocity fields were discussed.


  • [1] J.P. Boyd, N. Flyer, Compatibility conditions for time-dependent partial differential equations and the rate of convergence of Chebyshev and Fourier spectral methods, Comput. Methods Appl. Mech. Engrg. 175 (1999) 281-309.
  • [2] J. P. Boyd, A Comparison of Numerical Algorithms for Fourier Extension of the First, Second, and Third Kinds, J. Comput. Phys. 178 (2002) 118-160 .
  • [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. Springer-Verlag, 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) 521-527.
  • [6] W. E, J. Liu, Vorticity Boundary Condition and Related Issues for Finite Difference Schemes, J. Comput. Phys. 124 (1996) 368-382.
  • [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 Initial-Boundary Value Problems, SIAM J. Sci. Comput., 23 (5) (2002) 1731–1751.
  • [9] P.M. Gresho, R.L. Sani, On pressure boundary-conditions for the incompressible Navier–Stokes equations, Int. J. Numer. Methods Fluids 7 (10) (1987) 1111-1145.
  • [10] D.V. Le, B.C. Khoo, K.M. Lim, An implicit-forcing immersed boundary method for simulating viscous flows in irregular domains, Comput. Methods Appl. Mech. Engrg. 197 (2008) 2119-2130.
  • [11] L. Quartapelle, Vorticity Conditioning in the Computation of Two-Dimensional Viscous Flows, J. Comput. Phys., 40 (1981) 453-477.
  • [12] L. Quartapelle, Numerical solution of the incompressible Navier-stokes equations. Birkhauer-Verlag, Basel, 1993.
  • [13] D. Rempfer, On boundary conditions for incompressible Navier-Stokes problems. Appl. Mech. Rev. 59 (2006) 107-125.
  • [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) 55-74.
  • [15] F. Sabetghadam, M. Badri, S. Sharafatmandjoor, H. Kor, An Immersed Boundary Fourier Pseudo-spectral Method for Simulation of Confined Two-dimensional 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) 596-627.
  • [17] R. Temam, Behaviour at time t = 0 of the solutions of semi-linear evolution equations, J. Differ. Equ. 17 (1982) 73-92.
  • [18] R. Temam, Suitable initial conditions, J. Comput. Phys. 218 (2006) 443-450.
  • [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) 448-476.
  • [21] M. Vanella, P. Rabenold, E. Balaras, A direct-forcing embedded-boundary method with adaptive mesh refinement for fluid-structure interaction problems, Journal of Computational Physics 229 (2010) 6427-6449.
  • [22] K. Luo, Z. Wang, J. Fan, K. Cen, Full-scale solutions to particle-laden flows: Multidirect forcing and immersed boundary method, Physical Review E 76 (2007) 066709(1-9).
  • [23] Z. Wang, J. Fan and K. Cen, Immersed Boundary Method for the Simulation of 2D Viscous Flow Based on Vorticity-Velocity Formulations, J. Comput. Phys. 228 (2009) 1504-1520.
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description