- Email: [email protected]
Bivariate splines for fluid flows
Computers & Fluids 33 (2004) 1047–1073 www.elsevier.com/locate/compfluid
Bivariate splines for fluid flows Ming-Jun Lai *, Paul Wenston Department of Mathematics, The University of Georgia, Athens, GA 30602, USA Received 22 January 2003; received in revised form 22 September 2003; accepted 6 October 2003
Abstract We discuss numerical approximations of the 2D steady-state Navier–Stokes equations in stream function formulation using bivariate splines of arbitrary degree d and arbitrary smoothness r with r < d. We derive the discrete Navier–Stokes equations in terms of B-coefficients of bivariate splines over a triangulation, with curved boundary edges, of any given domain. Smoothness conditions and boundary conditions are enforced through Lagrange multipliers. The pressure is computed by solving a Poisson equation with Neumann boundary conditions. We have implemented this approach in MATLAB and our numerical experiments show that our method is effective. Numerical simulations of several fluid flows will be included to demonstrate the effectiveness of the bivariate spline method. Ó 2003 Elsevier Ltd. All rights reserved.
1. Introduction One of the major tasks for applied mathematicians is to develop efficient methods for numerical solutions of partial differential equations. The steady-state Navier–Stokes equations are one of the most important partial differential equations which have various applications. Although there are many computational methods available in the literature for the numerical solution of the Navier– Stokes equations, new and more efficient methods continue to be developed in order to increase the power of computational flow simulations. For example, in a recent paper [3], a collocation Bspline method was developed for the solution of the Navier–Stokes equations. In this paper, we use bivariate spline functions over arbitrary triangulations for numerical solution of 2D Navier– Stokes equations. (See [1] for the trivariate spline approximation of 3D Navier–Stokes equations.) Our approach is like the finite element method using triangles to approximate any given 2D
*
Corresponding author. Tel.: +1-706-542-2065, fax: +1-706-542-2573. E-mail addresses: [email protected] (M.-J. Lai), [email protected] (P. Wenston).
0045-7930/$ - see front matter Ó 2003 Elsevier Ltd. All rights reserved. doi:10.1016/j.compfluid.2003.10.003
1048
M.-J. Lai, P. Wenston / Computers & Fluids 33 (2004) 1047–1073
polygonal domains and using piecewise polynomials over triangulations to approximate the solution of the Navier–Stokes equations. The main different features are (1) no macro-element or locally supported spline functions are constructed; (2) polynomials of high degrees can be easily used to get a better approximation power; (3) smoothness can be imposed in a flexible way across the domain at places where the solution is expected to be smooth. For example, the solution of the steady-state Navier–Stokes equation is H 2 inside the domain and H 1 near the boundary (cf. [18]); (4) the mass and stiffness matrices can be assembled easily and these processes can be done in parallel; (5) the stream function formulation will be used and thus the spline approximation of the solution of Navier–Stokes equations satisfies the divergence-free condition exactly; (6) the matrices that arise are singular which is an important difference from the classical finite element method; (7) our spline method leads to a linear system of special structure. We introduce a special numerical method to solve such particularly structured linear systems. Let us first introduce the stream function formulation. Let X R2 be a simply connected T domain and u ¼ ðu1 ; u2 Þ be the planar velocity of a fluid flow over X. Also, let p be the pressure T T function, f ¼ ðf1 ; f2 Þ be the external body force of the fluid and g ¼ ðg1 ; g2 Þ be the velocity of the fluid flow on the boundary oX. Then the steady-state Navier–Stokes equations are 8 > < mDu þ ðu rÞu þ rp ¼ f; ðx; yÞ 2 X; div u ¼ 0; ðx; yÞ 2 X; ð1:1Þ > : u ¼ g; ðx; yÞ 2 oX; where D denotes the usual Laplace operator and r the gradient operator. After omitting the nonlinear terms, we have the steady-state StokesÕ equations 8 > < mDu þ rp ¼ f; ðx; yÞ 2 X; div u ¼ 0; ðx; yÞ 2 X; ð1:2Þ > : u¼g ðx; yÞ 2 oX: , u2 ¼ ou . Recall the fact that there exists a stream function u such that u ¼ curl u, i.e., u1 ¼ ou oy ox Such u is unique up to a constant (cf. [13, pp. 37–39]). Thus we may simplify the above Stokes and Navier–Stokes equations by cancelling the pressure term. Consider the Stokes equations first. Replacing u by curl u and then differentiating the first equation with respect to y and the second with respect to x, we subtract the first equation from the second one to obtain the following fourth-order equation mD2 u ¼ h with h ¼ ofox2 ofoy1 . Thus, the Stokes equations become a biharmonic equation
M.-J. Lai, P. Wenston / Computers & Fluids 33 (2004) 1047–1073
8 > mD2 u ¼ h > > > > ou > > ¼ g2 < ox ou > > > ¼ g1 > > oy > > : u ¼ b2
1049
in X; on oX; ð1:3Þ on oX; on oX;
where b2 is an anti-derivative of the tangential derivative of u along oX and will be examined in detail later. By a similar calculation, we easily see that the Navier–Stokes equations become the following fourth-order nonlinear equation: 8 o ou o2 u ou o2 u > 2 > > u mD > > oy oy oxoy ox oy 2 > > > > > o ou o2 u ou o2 u > > > ¼ h in X; > < ox oy ox2 ox oxoy ð1:4Þ ou > > ¼ g2 on oX; > > ox > > > > ou > > on oX; ¼ g1 > > > oy > : on oX: u ¼ b2 Let H 2 ðXÞ be the usual Sobolev space and H02 ðXÞ be the subspace of H 2 ðXÞ of functions whose derivatives of order less than or equal to one all vanish on the boundary oX. Define the bilinear form a2 ðu; wÞ and trilinear form qðh; u; wÞ by Z a2 ðu; wÞ ¼ Duðx; yÞDwðx; yÞ dxdy; X Z ouðx; yÞ owðx; yÞ ouðx; yÞ owðx; yÞ Dhðx; yÞ qðh; u; wÞ ¼ dxdy ox oy oy ox X and denote the L2 ðXÞ inner product by Z hh; wi ¼ hðx; yÞwðx; yÞ dxdy: X
We that say u 2 H 2 ðXÞ lowing: 8 ma2 ðu; wÞ ¼ hh; wi > > > > ou > > < ¼ g2 ox ou > > > ¼ g1 > > > : oy u ¼ b2
is a weak solution of the Stokes equations (1.3) if u satisfies the fol8w 2 H02 ðXÞ; on oX; on oX; on oX:
Similarly, a function u 2 H 2 ðXÞ is a weak solution of the Navier–Stokes equations (1.4) if u satisfies
1050
M.-J. Lai, P. Wenston / Computers & Fluids 33 (2004) 1047–1073
8 ma2 ðu; wÞ þ qðu; u; wÞ ¼ hh; wi > > > > > ou > > ¼ g2 < ox ou > > > ¼ g1 > > oy > > : u ¼ b2
8w 2 H02 ðXÞ; on oX; on oX; on oX:
Such weak formulations are referred as the stream function formulation of the Stokes and Navier–Stokes equations, respectively. It is known that the weak solution for StokesÕ equations exists and is unique for any m > 0. For the Navier–Stokes equations, such a weak solution exists for any m > 0, and is unique when m is sufficiently large (see, e.g., [13]). There are two major advantages using the stream function formulation over the traditional velocity–pressure formulation and vorticity–stream function formulation. Indeed, with the stream function formulation, we need to approximate only one stream function. Otherwise, we need to approximate two components of the velocity and one pressure function if the velocity–pressure formulation is used or one vorticity and one stream function if the vorticity–stream function formulation is used. In addition, the stream function formulation eliminates the pressure function which does not have appropriate boundary condition. In the vorticity–stream function formulation, the vorticity function does not have appropriate boundary condition. With bivariate spline functions of higher degrees, we are able to approximate stream functions very well. Thus, in this paper, we will use bivariate splines to approximate the stream function of the Stokes and Navier– Stokes equations. Next we discuss the computation of the pressure functions. By taking divergence, div, of the equations in (1.1) and (1.2), we can easily see that the pressure functions p of the Stokes and Navier–Stokes equations satisfy the following Poisson equations with nonhomogeneous Neumann boundary conditions involving the stream functions 8 in X < Dp ¼ divðfÞ ð1:5Þ op : ¼ mDðn curl uÞ þ n f on oX on for the Stokes equation and 8 < Dp ¼ divðfÞ þ div½ðcurl u rÞcurl ðuÞ op : ¼ mDðn curl uÞ þ n f þ n ½ðcurl u rÞðcurl uÞ on
in X; on oX
ð1:6Þ
for the Navier–Stokes equations. We will use these boundary conditions to compute p after we find the spline approximations of the stream function u. The paper is organized as follows: in Section 2, we first introduce bivariate spline spaces and in particular, the B-form representation of spline functions. The constrained minimization problems that will be used to generate approximate solutions of the stokes and Navier–Stoke equations lead, by way of Lagrange multipliers, to linear systems of the form T F L A k ¼ ; G 0 L c
M.-J. Lai, P. Wenston / Computers & Fluids 33 (2004) 1047–1073
1051
with A singular. We then introduce an iterative method to solve the above system and discuss its convergence. We continue in Section 3 with a discussion of the spline solution of the 2D Stokes equations, an equivalent biharmonic equation. In Section 4, we discuss the spline approximations of the 2D Navier–Stokes equations and prove the convergence of two methods for solving the discrete nonlinear equations. In Section 5, we present numerical results for standard benchmark flows such as the driven cavity flows, backward-facing step flows, and flows around circular obstacles. Also, we also present some other fluid flow experiments.
2. Preliminaries 2.1. Bivariate spline functions Given a bounded polygonal domain X 2 R2 , let M be a triangulation of X. Let d P 1 and r P 1 be two fixed integers. We introduce the spline spaces Sdr ðMÞ ¼ fs 2 C r ðXÞ; sjt 2 Pd 8t 2 Mg; where Pd denotes the space of bivariate polynomials of total degree d. In this paper, the B-form representation of splines on triangulations will be used (cf. [9] or [4]). Let T ¼ hv1 ; v2 ; v3 i be a nondegenerate triangle with vi ¼ ðxi ; yi Þ, i ¼ 1; 2; 3. It is well known that every point v ¼ ðx; yÞ can be written uniquely in the form v ¼ k1 v1 þ k2 v2 þ k3 v3
ð2:1Þ
k1 þ k2 þ k3 ¼ 1;
ð2:2Þ
with where k1 , k2 , and k3 are called the barycentric coordinates of the point v ¼ ðx; yÞ relative to the triangle T . Moreover each ki is a linear polynomial in x; y. Let Bdijk ðvÞ ¼
d! i j k kkk; i!j!k! 1 2 3
i þ j þ k ¼ d:
They are called the Bernstein–Bezier polynomials of degree d. In fact, the set Bd ¼ fBdijk ðx; y; zÞ; i þ j þ k ¼ dg is a basis for the space of polynomials Pd . As a consequence any polynomial p of degree d can be written uniquely in terms of Bdijk Õs, i.e., X p¼ cijk Bdijk : ð2:3Þ iþjþk¼d
The representation (2.3) for polynomials is referred to as the B-form with respect to T . Let iv1 þ jv2 þ kv3 Dd;T ¼ nijk ¼ ; i þ j þ k ¼ d; T 2 M ð2:4Þ d be a set of the domain points of degree d over triangulation M. For each spline function s 2 Sdr ðMÞ, since s restricted to each triangle T 2 M is a polynomial of degree d, we may write
1052
M.-J. Lai, P. Wenston / Computers & Fluids 33 (2004) 1047–1073
sjT ¼
X
cTijk Bdijk ;
T 2 M:
iþjþk¼d
Such a representation is called the B-form representation of the spline function s (cf. [4]). We denote by c :¼ fcTijk ; i þ j þ k ¼ d; T 2 Mg the B-coefficient vector of s. To evaluate a polynomial in B-form, there is the so-called de Casteljau algorithm which we now P ð0Þ describe. For p ¼ iþjþk¼d cijk Bdijk , let us write cijk ¼: cijk ðkÞ with k ¼ ðk1 ; k2 ; k3 Þ being the barycentric coordinates of v ¼ ðx; yÞ with respective to T and define for a positive integer r P 1 ðrÞ
ðr1Þ
ðr1Þ
ðr1Þ
cijk ðkÞ ¼ k1 ciþ1;j;k ðkÞ þ k2 ci;jþ1;k ðkÞ þ k3 ci;j;kþl ðkÞ: We have then X p¼
ðrÞ
cijk ðkÞBdr ijk ;
0 6 r 6 d:
iþjþk¼dr
In particular, for r ¼ d, we have ðdÞ
p ¼ c0;0;0 ðkÞ which is the value of p at v ¼ ðx; yÞ whose barycentric coordinates are k ¼ ðk1 ; k2 ; k3 Þ with k1 þ k2 þ k3 ¼ 1. We next discuss how to take derivatives of polynomials in B-form. We start with formulas for the directional derivatives of p in a direction defined by a vector u. We have X ð1Þ cijk ðaÞBd1 Du p ¼ d ijk iþjþk¼d1
with a ¼ ða1 ; a2 ; a3 Þ the T -coordinates of u; that is u ¼ a1 v1 þ a2 v2 þ a3 v3 with a1 þ a2 þ a3 ¼ 1: Note that if u ¼ v1 v2 is the direction vector, the T -coordinates of u are (1, )1, 0). In general, we have Dmu pðvÞ ¼
d! ðd mÞ!
X
ðmÞ
cijk ðaÞBdm ijkl ðvÞ:
ð2:5Þ
iþjþk¼dm
Note that for arbitrary direction vector u ¼ ðu1 ; u2 Þ 2 R2 Du p ¼ u1
o o p þ u2 p: ox oy
For a triangle T ¼ hv1 ; v2 ; v3 i, if we let u ¼ v2 v1 ¼ ðx2 x1 ; y2 y1 Þ, v ¼ v3 v1 ¼ ðx3 x1 ; y3 y1 Þ, it follows that:
M.-J. Lai, P. Wenston / Computers & Fluids 33 (2004) 1047–1073
ðy3 y1 Þ ðy2 y1 Þ Du p Dv p; 2AT 2AT ðx2 x1 Þ ðx3 x1 Þ Dv p Du p; Dy p ¼ 2AT 2AT
1053
Dx p ¼
ð2:6Þ
where AT is the area of T . There are precise formulas for the integrals and inner products of polynomials in B-form (cf. [5]). Lemma 2.1. Let p be a polynomial of degree d with B-coefficients cijk , i þ j þ k ¼ d on a triangle T . Then Z X AT pðx; yÞ dx dy ¼ cijk ; d þ 2 iþjþk¼d T 2 where AT ¼ area of T . Lemma 2.2. Let q be another polynomial with B-coefficients dijk , i þ j þ k ¼ d, the inner product of p and q over t is given by Z X AT iþr jþs kþt pðx; yÞqðx; yÞ dx dy ¼ cijk drst : ð2:7Þ i j k 2d 2d þ 2 t iþjþk¼d rþsþt¼d d 2 We next discuss the smoothness conditions for a spline function s in Sdr ðMÞ. These are wellknown conditions on the coefficients of s that will assure that s has certain global smoothness properties. Theorem 2.3. Let t ¼ hv1 ; v2 ; v3 i and t0 ¼ hv1 ; v2 ; v4 i be two triangles with common edge hv1 ; v2 i. Then s is of class C r on t [ t0 if and only if X 0 ctiþl;jþm;j Bml;m;j ðv4 Þ; m ¼ 0; . . . ; r; i þ j ¼ d m: ctijm ¼ lþmþj¼m
For a proof, see [9] or [4]. This theorem guarantees the existence of a matrix H such that s is in C ðXÞ if and only if r
H c ¼ 0; where c encodes the B-coefficients of s. 2.2. A matrix iterative method The constrained minimization problems that will be used to generate approximate solutions of the stokes and Navier–Stoke equations will lead, by way of Lagrange multipliers, to linear systems of the form (cf. [10] and [19])
1054
M.-J. Lai, P. Wenston / Computers & Fluids 33 (2004) 1047–1073
LT 0
F A k ¼ ; G L c
ð2:8Þ
with A singular and appropriate matrices L and vectors F and G. Computing a least squares solution of the above system would allow us to compute c when c is unique. However, when the size of a system is very large, any least squares method will be too costly. We present here an iterative method for solving the above system which involves matrices of smaller size. Consider the following sequence of problems: T F A kðkþ1Þ L ¼ ð2:9Þ G kðkÞ I L cðkþ1Þ for k ¼ 0; 1; . . ., with an initial guess kð0Þ , e.g., kð0Þ ¼ 0, and I the identity matrix. Assume that the size of L is m n and the size of A is n n with m < n. Note that (2.9) reads Acðkþ1Þ þ LT kðkþ1Þ ¼ F ; Lcðkþ1Þ kðkþ1Þ ¼ G kðkÞ :
ð2:10Þ
Multiplying on the left the second equation in (2.10) by LT , we get LT Lcðkþ1Þ LT kðkþ1Þ ¼ LT G LT kðkÞ or LT kðkþ1Þ ¼ 1 LT Lcðkþ1Þ 1 LT G þ LT kðkÞ . We substitute this into the first equation in (2.10) to get 1 T 1 A þ L L cðkþ1Þ ¼ F þ LT G LT kðkÞ : ð2:11Þ Using again the first equation in (2.10), i.e., AcðkÞ ¼ F LT kðkÞ to replace F in (2.11), we have 1 T 1 A þ L L cðkþ1Þ ¼ AcðkÞ þ LT G for all k P 1. This leads to the following. Algorithm 2.4. Fix > 0. Given an initial guess kð0Þ 2 ImðLÞ, e.g., kð0Þ ¼ 0, we define cð1Þ by 1 1 T 1 T ð1Þ T ð0Þ F þ L GL k c ¼ Aþ L L and iteratively define 1 1 T 1 T ðkþ1Þ ðkÞ ¼ Aþ L L Ac þ L G c
ð2:12Þ
for k ¼ 1; 2; . . ., where ImðLÞ is the range of L. This algorithm was briefly discussed in [12]. A convergence result was pointed out. That is, kc ckþ1 k 6 Ckc cðkÞ k for some constant C > 0. According to this result for convergence one clearly needs to choose small enough to guarantee that C < 1. Here we present the following convergence result (cf. [2] for a proof) which ensures the convergence for any > 0.
M.-J. Lai, P. Wenston / Computers & Fluids 33 (2004) 1047–1073
1055
Theorem 2.5. Let A ¼ As þ Aa where As ¼ 12ðA þ AT Þ is the symmetric part of A and Aa ¼ 12ðA AT Þ is the anti-symmetric part of A. Furthermore, suppose that As is positive definite with respect to L, that is, xT As x ¼ 0 and Lx ¼ 0 imply that x ¼ 0. Then for any > 0, the matrix 1 A þ LT L is invertible. That is, the Algorithm 2.4 is well defined. Furthermore, there exists positive constants CðÞ > 1 and cðÞ < 1 depending on but independent of k such that kc cðkþ1Þ k 6 CðÞðcðÞÞ
kþ1
for k P 1. Our numerical experiments show that it gives the same result as the least squares method. Since any least squares method is computationally expensive, we recommend the iterative method for large problems. 3. Spline approximations of the Stokes equations In this section, we consider spline approximations of the 2D Stokes equations in the stream function formulation. That is, we consider the following biharmonic equation with nonhomogeneous boundary conditions g1 and g2 : 8 ma2 ðu; wÞ ¼ hh; wi 8w 2 H02 ðXÞ; > > > > ou > > < on oX; ¼ g2 ox ð3:1Þ ou > > on oX; ¼ g1 > > > > : oy u ¼ b2 on oX; where b2 is an anti-derivative of the tangential derivative of u along oX which can be computed based on g1 and g2 . For simplicity, let b1 be the normal derivative of u along oX. Note that b1 can be easily computed using g1 and g2 . It is well known that the weak solution of (3.1) exists and is unique (cf. [13]). We are going to approximate the weak solution of (3.1) by using bivariate spline functions. Let M be a triangulation of X if X is a polygonal domain. Otherwise, we choose vertices v1 ; . . . ; vn on oX which include all the corner points of oX. Instead of the piecewise linear interpolation of the vertices, we use a C 0 piecewise quadratic spline curve soX , which interpolates the vertices, to approximate the boundary oX. Clearly, we may assume that soX is a much better approximation of oX than the piecewise linear interpolant based on an appropriate choice of v1 ; . . . ; vn . Adding interior vertices vnþ1 ; . . . ; vN inside X, we triangulate X using the vertices v1 ; . . . ; vN . For convenience, we let M be the triangulation consisting of all interior triangles and boundary triangles with piecewise quadratic edges. Let d and r be two positive integer with d > r and let S Sd0 ðMÞ be the space of spline functions which are C r inside of X. Note that S H 2 ðXÞ. Recall from Section 2.1 that s 2 Sd0 ðMÞ can be given by
1056
M.-J. Lai, P. Wenston / Computers & Fluids 33 (2004) 1047–1073
X
sðx; yÞ ¼
cti;j;k Bti;j;k ðx; yÞ; 8ðx; yÞ 2 t 2 M:
iþjþk¼d
ðcti;j;k ; i
þ j þ k ¼ d; t 2 MÞ be the B-coefficient vector of s. Recall from Section 2.1 that Let c ¼ there is a matrix H such that if s 2 S, then H c ¼ 0: Although in general we may not be able to find a spline function s 2 S satisfying the boundary conditions exactly, we can find a spline Sb 2 S approximating the boundary conditions. If the the boundary conditions are continuous, the interpolation is the obvious choice. Assuming that v1 ; . . . ; vb are arranged in the counterclockwise direction, Sb interpolates b2 at d þ 1 distinct points on any straight boundary edge hvi ; viþ1 i and at 2d þ 1 distinct points on any quadratic boundary b curved edge vi ^ viþ1 . Also oS interpolates b1 at d distinct points on any straight boundary edge oni hvi ; viþ1 i and at 2d 1 distinct points on any quadratic boundary curved edge vi ^ viþ1 , where ni denotes the outer normal direction of boundary edge hvi ; viþ1 i or vi ^ viþ1 . If the boundary conditions are not continuous, we can still use interpolation except that the interpolation values now come from a conveniently chosen smooth approximation of the boundary conditions. Since the values of any spline and its normal derivative at any of the interpolation points are linear combinations of the B-coefficients of the spline, we may use Bc ¼ b
ð3:2Þ
to express the approximate boundary conditions. Let S0 ¼ S \ H02 ðXÞ. Our spline approximation of the biharmonic equation (3.1) is to find sh;b 2 S satisfying the approximate boundary condition (3.2) such that ma2 ðsh;b ; sÞ ¼ hh; si; 8s 2 S0 :
ð3:3Þ
Using the Lax–Milgram Theorem, problem (3.3) has an unique solution sh;b in S. Based on the theory of elliptic equations (cf. [6]), solving problem (3.3) is equivalent to finding c which minimizes the following energy functional m EðsÞ ¼ a2 ðs; sÞ hh; si 2 subject to s 2 S satisfying the approximation conditions. In terms of B-coefficients, the above energy functional can be approximated by m ~ EðcÞ :¼ cT Ac cT Mh 2 subject to Hc ¼ 0 and Bc ¼ b, where A ¼ diagðAt ; t 2 MÞ is a bending matrix with blocks Z DBti;j;k DBtl;m;n dx dy ; t 2 M; At ¼ iþjþk¼d lþmþn¼d
t
M ¼ diagðMt ; t 2 MÞ is a mass matrix with blocks Z t t Bi;j;k Bl;m;n dx dy ; t2M Mt ¼ t
iþjþk¼d lþmþn¼d
M.-J. Lai, P. Wenston / Computers & Fluids 33 (2004) 1047–1073
1057
and h ¼ ðhtl;m;n ; l þ m þ n ¼ d; t 2 MÞ is a vector with X hl;m;n Btl;m;n dx dy 2 Sd0 ðMÞ sh ¼ lþmþn¼d
the spline interpolation of h over the domain points fntl;m;n ; l þ m þ n ¼ d; t 2 Mg. We note that the mass and bending matrices, M and A, are diagonal block matrices whose blocks can be assembled trivially and in parallel. The matrix H is sparse with each row involving at most ðr þ 1Þðr þ 2Þ=2 þ 1 nonzero elements. Also, the matrix B is very sparse with each row containing at most d þ 1 nonzero elements. Following the Lagrange multiplier method, we let m Lðc; a; bÞ ¼ cT Ac þ aT H c þ bT Bc cT Mh 2 and compute local minimizers. It follows that we need to solve the following: 2 T 32 3 2 3 a Mh H BT mA 4 0 0 H 54 b 5 ¼ 4 0 5: c b 0 0 B We note that the existence and uniqueness of sh;b implies that there exists a unique solution c satisfying the above linear system. Clearly, A is positive definite with respect to ½H; B. Thus, our iterative method introduced in Section 2.2 can be applied to this system and the iterative solutions converge. Next we use the approach presented above to compute the pressure function assuming that we have found the approximation of the stream function u. That is, we need to solve the Poisson problem with Neumann boundary conditions (1.5). Note that we are seeking the pressure in Z 2 2 L0 ðXÞ ¼ p 2 L ðXÞ; p ¼ 0 : X
Again we can use the spline space S as an approximating space. For a spline approximation sp of p with B-coefficient vector cp , its integral over a triangle is simply a weighted sum of its Bcoefficients (cf. Lemma 2.1), hence the spline approximation of the zero mean value condition, R p ¼ 0, can be written X Uc ¼ 0 with U a vector of size N . The pressure is also the solution of a minimization problem, namely that of minimizing the functional Z Z Z 1 2 jDvj ðhÞv ðf nv þ mDðcurl u nÞvÞ P ðvÞ ¼ 2 X X oX over L20 ðXÞ. We now write P in terms of cp . Let Gp interpolate f n þ mðDcurl uÞ n on the boundary. The later is computed by using the spline approximation sh;b of the solution u of the Stokes equations. Let R be the matrix which restricts the B-coefficients c of a spline over X to the B-coefficients of the spline restricted to oX. We have
1058
M.-J. Lai, P. Wenston / Computers & Fluids 33 (2004) 1047–1073
1 P ðcp Þ ¼ ðcp ÞT Kcp hT Mcp GTp Mb Rcp ; 2 where K is the stiffness matrix which is a diagonal block matrix similar to A and Mb is the mass matrix for the bivariate splines of degree d over the edges of oX. Following the Lagrange multiplier method, we need to solve the following linear system: 3 2 32 3 2 k1 U HT K Mh þ RT Mb RGp 5: 40 0 0 H 54 k2 5 ¼ 4 c 0 0 0 U p Again, it is easy to see that K is positive definite with respect to ½H ; U . The iterative method described in Section 2.2 converges.
4. Spline approximations of the Navier–Stokes equations In this section, we consider the bivariate spline approximation of the 2D Navier–Stokes equations in the stream function formulation. Let u 2 H 2 ðXÞ be the weak solution of the Navier– Stokes equations 8 ma2 ðu; wÞ þ qðu; u; wÞ ¼ hh; wi 8w 2 H02 ðXÞ; > > > > ou > > < on oX; ¼ g2 ox ð4:1Þ ou > > on oX; ¼ g > 1 > > > : oy u ¼ b2 on oX: It is known that the weak solution u exists and is unique if m is sufficiently large or h is sufficiently small (cf. [13]). We shall use the approach discussed in the previous section to derive spline approximations of the stream function u. With the same spline space S and the same notation as in the previous section, we have a2 ðu; wÞ ¼ cT Ad; where c and d stand for the B-coefficient vectors of spline functions u; w 2 S. The difference with the previous section is the presence of the nonlinear term, the trilinear form q defined by Z ouðx; yÞ owðx; yÞ ouðx; yÞ owðx; yÞ Dhðx; yÞ dxdy: ð4:2Þ qðh; u; wÞ ¼ ox oy oy ox X Let us first get acquainted with the trilinear form (4.2). Let e encode the B-coefficient vector of the spline function h in S. We can, by tedious computations, show that there is a matrix QðeÞ, which is a linear function of e, such that qðh; u; wÞ ¼ dT QðeÞc: There are several useful properties of QðeÞ given below. Lemma 4.1. QðeÞ is anti-symmetric, i.e., QðeÞT ¼ QðeÞ.
ð4:3Þ
M.-J. Lai, P. Wenston / Computers & Fluids 33 (2004) 1047–1073
1059
Proof. We first note that qðh; u; wÞ ¼ qðh; w; uÞ. In terms of B-coefficient vectors, we have dT QðeÞc ¼ cT QðeÞd: T
T
Thus, it follows from dT QðeÞ c ¼ cT QðeÞd that dT QðeÞ c ¼ dT QðeÞc for all c and d. Hence, QðeÞT ¼ QðeÞ. h Consequently, we have Lemma 4.2. Fix e. dT QðeÞd ¼ 0 for any d satisfying H d ¼ 0. Next we will need the following. Lemma 4.3. There exists a constant C1 such that jdT QðeÞcj 6 C1 kck kdk kek:
Proof. The Markov inequality and a tedious computation yields the above estimate. h The coefficient vector c of the spline approximation of the weak solution satisfies Hc ¼ 0 and Bc ¼ b and mcT Ad þ cT QðcÞd ¼ hT Md
ð4:4Þ
for all coefficient vectors d of splines in S0 . The existence of c can be shown by using the same argument for the existence of the weak solution / satisfying (4.1) (cf. [13,16]). We are mainly interested in computing c. We now observe that if c and k are chosen so that 8 < mAc QðcÞc þ ½ H T BT k ¼ Mh ð4:5Þ Hc ¼ 0 : Bc ¼ b then mcT Ad þ cT QðcÞd þ kT
H d hT Md ¼ 0 B
for all vectors d. Since both H d ¼ 0 and Bd ¼ 0 for all coefficient vectors d of splines in S0 (4.4) follows. We next derive two methods to linearize the nonlinear equations (4.5), using a simple iteration algorithm and NewtonÕs method. Algorithm 4.4 (A simple iteration algorithm). Let ðcð0Þ be the solution of the linear problem (i.e. the associated Stokes equations) and for n ¼ 0; 1; . . ., define ðcðnþ1Þ ; kðnþ1Þ Þ as the solution of
1060
M.-J. Lai, P. Wenston / Computers & Fluids 33 (2004) 1047–1073
mAcðnþ1Þ QðcðnÞ Þcðnþ1Þ þ ½ H T
BT kðnþ1Þ ¼ Mh;
H cðnþ1Þ ¼ 0;
ð4:6Þ
Bcðnþ1Þ ¼ b: We use Algorithm 2.4 to find the new iterates cðnþ1Þ of Algorithm 4.4. Note that A is symmetric and positive definite with respect to H and B and that QðcðnÞ Þ is anti-symmetric by Lemma 4.1. Theorem 2.5 ensures that the matrix iterative method converges for the linear system in (4.6). To apply NewtonÕs method to solve (4.5), we consider the mapping C defined by C : ðc; kÞ 7! ðmAc QðcÞc þ ½ H T
BT k Mh; H c; Bc bÞ
and seek to solve Cðc; kÞ ¼ 0. We write XðnÞ ¼ ðcðnÞ ; kðnÞ Þ and let Xð0Þ be the solution of the linear Stokes problem. Then we define Xðnþ1Þ by the equation C0 ðXðnÞ ÞðXðnþ1Þ XðnÞ Þ ¼ CðXðnÞ Þ: Thanks to the bilinearity of the mapping ðc; dÞ ! QðcÞd, the above equation leads to mAcðnþ1Þ QðcðnÞ Þcðnþ1Þ Qðcðnþ1Þ cðnÞ ÞcðnÞ þ ½ H T
BT kðnþ1Þ ¼ Mh;
H cðnþ1Þ ¼ 0; ðnþ1Þ
Bc
ð4:7Þ
¼ b:
e ðnÞ Þðcðnþ1Þ cðnÞ Þ with Qðc e ðnÞ ÞcðnÞ ¼ It can be shown that Qðcðnþ1Þ cðnÞ ÞcðnÞ can be written Qðc e ðnÞ Þ. We have our second algorithm. QðcðnÞ ÞcðnÞ for some matrix Qðc Algorithm 4.5 (NewtonÕs method). Let ðcð0Þ Þ be the solution of the Stokes equations and for n ¼ 0; 1; 2; . . . define ðcðnþ1Þ ; kðnþ1Þ Þ to be the solution of e ðnÞ Þcðnþ1Þ þ ½ H T mAcðnþ1Þ QðcðnÞ Þcðnþ1Þ Qðc
e ðnÞ ÞcðnÞ ; BT kðnþ1Þ ¼ Mh Qðc
H cðnþ1Þ ¼ 0; ðnþ1Þ
Bc
ð4:8Þ
¼ b:
For the linear system (4.8), Algorithm 2.4 converges by Theorem 2.5 because QðcðnÞ Þ is antie ðnÞ Þ is positive definite with respect to ½H ; B for m sufficiently large. Indeed, symmetric and mA þ Qðc since A is positive definite with respect to H , we have xT Ax P a0 kxk2 for any x satisfying H x ¼ 0. It follows from Lemma 4.3, e ðnÞ Þx P ma0 kxk2 C1 kcðnÞ k kxk2 : mxT Ax þ xT Qðc By Lemma 4.6, the sequence kcðnÞ k is bounded and hence e ðnÞ Þx P ðma0 C1 C2 Þkxk2 : mxT Ax þ xT Qðc e ðnÞ Þ is positive definite with respect to ½H ; B. Therefore, for m sufficiently large, mA þ Qðc
M.-J. Lai, P. Wenston / Computers & Fluids 33 (2004) 1047–1073
1061
Lemma 4.6. Suppose that cðnþ1Þ is the solution of (4.8). Then kcðnþ1Þ k 6 C for a positive constant C independent of n. Proof. Let c be the solution of (4.5). Let eðnþ1Þ ¼ c cðnþ1Þ and sðnþ1Þ ¼ k kðnþ1Þ . Then we have mAeðnþ1Þ þ QðcðnÞ Þcðnþ1Þ þ Qðcðnþ1Þ ÞcðnÞ QðcÞc þ ½H ; BT sðnþ1Þ ¼ QðcðnÞ ÞcðnÞ ; which gives mAeðnþ1Þ QðcðnÞ Þeðnþ1Þ QðeðnÞ Þc þ Qðcðnþ1Þ ÞcðnÞ þ ½H ; BT sðnþ1Þ ¼ QðcðnÞ ÞcðnÞ or mAeðnþ1Þ QðcðnÞ Þeðnþ1Þ QðeðnÞ Þc Qðeðnþ1Þ ÞcðnÞ þ QðeðnÞ ÞcðnÞ þ ½H ; BT sðnþ1Þ ¼ 0: We multiply this equation on the left by ðeðnþ1Þ ÞT and get mðeðnþ1Þ ÞT Aeðnþ1Þ ¼ ðeðnþ1Þ ÞT QðeðnÞ ÞeðnÞ þ ðeðnþ1Þ ÞT Qðeðnþ1Þ ÞcðnÞ : It follows that
ma0 keðnþ1Þ k 6 C1 keðnÞ k2 þ keðnþ1Þ kkcðnÞ k 6 C1 keðnÞ k2 þ keðnþ1Þ kðkeðnÞ k þ kckÞ :
1 kck We now use an induction to prove the desired result. Assume that keðnÞ k 6 K with K 6 12ma0 C . C1 This is always possible if m is large enough (cf. [14]). Then we have
ðma0 C1 ðK þ kckÞÞkeðnþ1Þ k 6 C1 keðnÞ k2 or keðnþ1Þ k 6
C1 K 2 6 K: ma0 C1 ðK þ kckÞ
This completes the proof with C2 ¼ K þ kck. h 5. Numerical simulation of various fluid flows In this section we first present our numerical experiments with three benchmark flow examples: the driven cavity flow, the backward facing step flow, and the flow around a circular object. We show the streamlines, cavity and pressure of the driven cavity flows, the streamlines for other two kinds of flows. Then we present some other examples of fluid flows. All the figures are shown in the Cartesian coordinates with the horizontal direction being the x-axis and the vertical direction being the y-axis. 5.1. Computational experiments on driven cavity flows Example 5.1. Let us consider the standard cavity flow over unit square domain ½0; 1 ½0; 1. The boundary conditions are u ¼ ðu1 ; u2 Þ ¼ ð0; 0Þ for all four line boundary pieces except for
1062
M.-J. Lai, P. Wenston / Computers & Fluids 33 (2004) 1047–1073
u1 ðx; 1Þ ¼ 1 when 0 6 x 6 1. With Reynolds numbers are 10, 100, 1000, 5000, 10000 and 20,000, the stream lines of the cavity flow are shown in Figs. 2–5. We use several different triangulations for the computation. Two typical triangulations are shown in Fig. 1. We use the splines of degree 8 or higher and smoothness 2. Our main reason to use splines of higher order is to achieve the best approximation property (cf. [15]). All the computation is done using a PC with MS Window 2000. Certainly, our method can compute the stream functions for much higher Reynolds numbers. The stream lines, vorticity, and pressure are similar based on different triangulations, especially when Reynolds numbers are small. We use the triangulation (the right one in Fig. 1) to capture the eddies at the corners.
1
1
0.9
0.9
0.8
0.8
0.7
0.7
0.6
0.6
0.5
0.5
0.4
0.4
0.3
0.3
0.2
0.2
0.1
0.1
0
0 0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
Fig. 1. Two typical triangulations of the square domain.
Fig. 2. The stream lines of driven cavity flow (Reynolds number ¼ 10 and 100).
0.9
1
M.-J. Lai, P. Wenston / Computers & Fluids 33 (2004) 1047–1073
1063
Fig. 3. The stream lines of driven cavity flow (Reynolds number ¼ 1000 and 5000).
Fig. 4. The stream lines of driven cavity flow (Reynolds number ¼ 10,000 and 20,000).
It can be seen that the stream lines in Figs. 2 and 3 are similar to those computed by using the finite difference method given in [11]. The stream lines in Fig. 4 are very similar to the one in [7] which is computed based on time dependent Navier–Stokes equations solved by a finite difference method. Also, the vorticity of the flows wðx; yÞ ¼
o2 o2 /ðx; yÞ þ /ðx; yÞ; ox2 oy 2
1064
M.-J. Lai, P. Wenston / Computers & Fluids 33 (2004) 1047–1073
are shown as in Figs. 5–7. They are very similar to those in [11]. The vorticity in Fig. 7 with Reynolds number ¼ 10,000 is similar to the one in [8] which is computed based on time dependent Navier–Stokes equations and a fourth order finite difference method. Next we compute the pressure of the flows which is normalized such that the integral of the pressure over the domain is zero. As we pointed out before, the pressure is computed by solving the Neumann problem of the standard Poisson equation after the stream function is obtained. Again we set up as a minimization problem which minimizes the stiff matrix subject to the smoothness conditions, boundary conditions, and a normalization condition (cf. [2] for multivariate spline method for numerical solution of partial differential equations). The pressures for flows with different Reynolds numbers are given in Figs. 8–10.
Fig. 5. The contour of the vorticity of the driven cavity flow (Reynolds number ¼ 10 and 100).
Fig. 6. The contour of the vorticity of the driven cavity flow (Reynolds number ¼ 1000, 5000).
M.-J. Lai, P. Wenston / Computers & Fluids 33 (2004) 1047–1073
1065
Fig. 7. The contour of the vorticity of the driven cavity flow (Reynolds number ¼ 10,000 and 20,000).
Fig. 8. The contour of the pressure of the driven cavity flow (Reynolds number ¼ 10 and 100).
5.2. Computational experiments on backward facing step flows Next we consider the well-known backward facing step flows. The inflow condition at the left inlet is the parabolic flow with maximum value 1. The outflow condition at the right outlet is Outflow condition. Both velocity components remain the same in the direction normal to the boundary. That is o2 u¼0 on2
and
o o u¼0 on os
1066
M.-J. Lai, P. Wenston / Computers & Fluids 33 (2004) 1047–1073
Fig. 9. The contour of the pressure of the driven cavity flow (Reynolds number ¼ 1000 and 5000).
Fig. 10. The contour of the pressure of the driven cavity flow (Reynolds number ¼ 10,000 and 20,000).
at the outlet, where s; n denote the tangential and normal direction of the outlet. Figs. 11–14 are some numerical simulations of the backward-facing step flows. For detailed explanation of these flows, see, e.g. [11] and the reference therein. The streamlines of the flows from our numerical experiments match with those in [11]. 5.3. Computational experiments of flows around circular obstacle In addition to the no-slip and outflow boundary conditions in the previous sections, we need the following freeslip conditions for simulating the flows around some obstacles.
M.-J. Lai, P. Wenston / Computers & Fluids 33 (2004) 1047–1073
Fig. 11. Flows over a backward facing step (Reynolds number ¼ 1).
Fig. 12. Flows over a backward facing step (Reynolds number ¼ 100).
Fig. 13. Flows over a backward facing step (Reynolds numbers ¼ 250).
Fig. 14. Flows over a backward facing step (Reynolds number ¼ 500).
Fig. 15. Flow around a disk with Reynolds number ¼ 10.
Fig. 16. Flow around a disk with Reynolds number ¼ 50.
Fig. 17. Flow around a disk with Reynolds number ¼ 100.
1067
1068
M.-J. Lai, P. Wenston / Computers & Fluids 33 (2004) 1047–1073
Fig. 18. Triangulation of the triangular domain and the stream lines of the triangular driven cavity flow (Reynolds number ¼ 100).
Fig. 19. Contour of the vorticity and pressure of the driven cavity flow (Reynolds number ¼ 100).
Fig. 20. The stream lines of the triangular driven cavity flow (Reynolds number ¼ 1000).
M.-J. Lai, P. Wenston / Computers & Fluids 33 (2004) 1047–1073
1069
Fig. 21. Contour of the vorticity and pressure of the triangular driven cavity flow (Reynolds number ¼ 1000).
Fig. 22. The stream lines of the triangular driven cavity flow (Reynolds number ¼ 5000).
Fig. 23. Contour of the vorticity and pressure of the triangular driven cavity flow (Reynolds number ¼ 5000).
1070
M.-J. Lai, P. Wenston / Computers & Fluids 33 (2004) 1047–1073
Fig. 24. The flows around a boxy car model and a streamlined car model.
Fig. 25. The pressures around a boxy car model and a streamlined car model.
M.-J. Lai, P. Wenston / Computers & Fluids 33 (2004) 1047–1073
1071
Fig. 26. Flows through sudden contraction and enlargement.
Freeslip conditions. No fluid penetrates the boundary, i.e. o u¼0 on
and
o o u¼0 on os
along a part of the boundary oX. Here are some numerical simulations of flows around a circular disk for various Reynolds numbers (Figs. 11–17). The in-flows are at a constant velocity at the left vertical boundary, i.e., u ¼ 1 and v ¼ 0. These flows can be compared with the ones in [11,20]. 5.4. Computational experiments on other interested flows We first show driven cavity flows over a triangular domain. Example 5.2. Let us consider the cavity flow over triangular domain. We uniformly refine the triangle as shown in Fig. 18 and use bivariate splines of degree 8 and smoothness 2. The boundary conditions are u ¼ ðu1 ; u2 Þ ¼ ð0; 0Þ for all three line boundary pieces except for u1 ðx; 1Þ ¼ 1 when 0 6 x 6 1. With Reynolds numbers are 100, 1000, and 5000, the stream lines of the cavity flow are shown in Figs. 18, 20 and 22, the contours of vorticity and pressure are shown in Figs. 19, 21 and 23. Next we present some additional examples which need inflow and outflow conditions. Example 5.3. In this example, we compare with the flows around a boxy car model and streamlined car model as in Fig. 24. The flows are similar to the real simulations as compiled in [17]. The pressures of two flows are computed as shown in Fig. 25. The highest pressure in the front of the boxy car is about 140 while that in the front of the smooth car is about 80. These clearly show that the boxy car has more wind resistance than the smooth car. Example 5.4. In this example, we consider the flow through sudden contraction and enlargement. We compare with the flows with two levels of contraction and enlargement. The flows are constant at the left vertical boundary and are the same for the two flows in Fig. 26, but result two different flow patterns and pressure levels as shown in Fig. 27.
1072
M.-J. Lai, P. Wenston / Computers & Fluids 33 (2004) 1047–1073
Fig. 27. Pressures of the flows through sudden contraction and enlargement.
References [1] Awanou G, Lai MJ. Trivariate spline approximations of the 3D Navier–Stokes equations, submitted. [2] Awanou G, Lai MJ, Wenston P. The multivariate spline method for numerical solutions of partial differential equations and scattered data fitting, submitted. [3] Botella O. On a collocation B-spline method for the solution of the Navier–Stokes equations. Comput Fluids 2002;13:397–420. [4] de Boor C. B-form basics. In: Farin G, editor. Geometric modelling. Philadelphia: SIAM Publication; 1987. p. 131– 48. [5] Chui CK, Lai MJ. Multivariate vertex splines and finite elements. J Approx Theory 1990;60:245–343. [6] Evans L. Partial Differential Equations. Providence: American Mathematical Society; 1998.
M.-J. Lai, P. Wenston / Computers & Fluids 33 (2004) 1047–1073
1073
[7] Weinan E, Liu J-G. Vorticity boundary condition and related issues for finite difference schemes. J Comput Phys 1996;124:368–82. [8] Weinan E, Liu J-G. Finite difference schemes for incompressible flows in the velocity impulse density formulation. J Comput Phys 1997;130:67–76. [9] Farin G. Triangular Bernstein–Bezier patches. Comput Aided Geom Design 1986;3:83–127. [10] Fortin M, Glowinski R. Augmented Lagrangian methods. Amsterdam: Elsevier; 1983. [11] Griebel M, Dornseifer T, Neunhoeffer T. Numerical simulation in fluid dynamics. Philadelphia: SIAM Publications; 1998. [12] Gunzburger M. Finite element methods for viscous incompressible flows. New York: Academic Press; 1989. [13] Giraultt V, Raviart PA. Finite element method for Navier–Stokes equations. Berlin: Springer-Verlag; 1986. [14] Karakashian O. On a Galerkin–Lagrange multiplier method for the stationary Navier–Stokes equations. SIAM J Numer Anal 1982;19:909–23. [15] Lai MJ, Schumaker LL. Approximation power of bivariate splines. Adv Comput Math 1998;9:251–79. [16] Lai MJ, Wenston P. Bivariate spline method for numerical solution of Navier–Stokes equations over polygons in stream function formulation. Numer Methods PDE 2000;16:147–83. [17] Nakayama Y. Visualized flow. Oxford: Pergamon Press; 1978. [18] Serrin J. On the interior regularity of weak solution of Navier–Stokes equations. Arch Rat Mech Anal 1962;9:187– 95. [19] Temam R. Navier–Stokes equations. In: Theory and numerical analysis. Amsterdam: North-Holland Publishing Co.; 1984. [20] Van Dyke M. An album of fluid motion. Stanford, CA: Parabolic Press; 1982.
Bivariate splines for fluid flows Ming-Jun Lai *, Paul Wenston Department of Mathematics, The University of Georgia, Athens, GA 30602, USA Received 22 January 2003; received in revised form 22 September 2003; accepted 6 October 2003
Abstract We discuss numerical approximations of the 2D steady-state Navier–Stokes equations in stream function formulation using bivariate splines of arbitrary degree d and arbitrary smoothness r with r < d. We derive the discrete Navier–Stokes equations in terms of B-coefficients of bivariate splines over a triangulation, with curved boundary edges, of any given domain. Smoothness conditions and boundary conditions are enforced through Lagrange multipliers. The pressure is computed by solving a Poisson equation with Neumann boundary conditions. We have implemented this approach in MATLAB and our numerical experiments show that our method is effective. Numerical simulations of several fluid flows will be included to demonstrate the effectiveness of the bivariate spline method. Ó 2003 Elsevier Ltd. All rights reserved.
1. Introduction One of the major tasks for applied mathematicians is to develop efficient methods for numerical solutions of partial differential equations. The steady-state Navier–Stokes equations are one of the most important partial differential equations which have various applications. Although there are many computational methods available in the literature for the numerical solution of the Navier– Stokes equations, new and more efficient methods continue to be developed in order to increase the power of computational flow simulations. For example, in a recent paper [3], a collocation Bspline method was developed for the solution of the Navier–Stokes equations. In this paper, we use bivariate spline functions over arbitrary triangulations for numerical solution of 2D Navier– Stokes equations. (See [1] for the trivariate spline approximation of 3D Navier–Stokes equations.) Our approach is like the finite element method using triangles to approximate any given 2D
*
Corresponding author. Tel.: +1-706-542-2065, fax: +1-706-542-2573. E-mail addresses: [email protected] (M.-J. Lai), [email protected] (P. Wenston).
0045-7930/$ - see front matter Ó 2003 Elsevier Ltd. All rights reserved. doi:10.1016/j.compfluid.2003.10.003
1048
M.-J. Lai, P. Wenston / Computers & Fluids 33 (2004) 1047–1073
polygonal domains and using piecewise polynomials over triangulations to approximate the solution of the Navier–Stokes equations. The main different features are (1) no macro-element or locally supported spline functions are constructed; (2) polynomials of high degrees can be easily used to get a better approximation power; (3) smoothness can be imposed in a flexible way across the domain at places where the solution is expected to be smooth. For example, the solution of the steady-state Navier–Stokes equation is H 2 inside the domain and H 1 near the boundary (cf. [18]); (4) the mass and stiffness matrices can be assembled easily and these processes can be done in parallel; (5) the stream function formulation will be used and thus the spline approximation of the solution of Navier–Stokes equations satisfies the divergence-free condition exactly; (6) the matrices that arise are singular which is an important difference from the classical finite element method; (7) our spline method leads to a linear system of special structure. We introduce a special numerical method to solve such particularly structured linear systems. Let us first introduce the stream function formulation. Let X R2 be a simply connected T domain and u ¼ ðu1 ; u2 Þ be the planar velocity of a fluid flow over X. Also, let p be the pressure T T function, f ¼ ðf1 ; f2 Þ be the external body force of the fluid and g ¼ ðg1 ; g2 Þ be the velocity of the fluid flow on the boundary oX. Then the steady-state Navier–Stokes equations are 8 > < mDu þ ðu rÞu þ rp ¼ f; ðx; yÞ 2 X; div u ¼ 0; ðx; yÞ 2 X; ð1:1Þ > : u ¼ g; ðx; yÞ 2 oX; where D denotes the usual Laplace operator and r the gradient operator. After omitting the nonlinear terms, we have the steady-state StokesÕ equations 8 > < mDu þ rp ¼ f; ðx; yÞ 2 X; div u ¼ 0; ðx; yÞ 2 X; ð1:2Þ > : u¼g ðx; yÞ 2 oX: , u2 ¼ ou . Recall the fact that there exists a stream function u such that u ¼ curl u, i.e., u1 ¼ ou oy ox Such u is unique up to a constant (cf. [13, pp. 37–39]). Thus we may simplify the above Stokes and Navier–Stokes equations by cancelling the pressure term. Consider the Stokes equations first. Replacing u by curl u and then differentiating the first equation with respect to y and the second with respect to x, we subtract the first equation from the second one to obtain the following fourth-order equation mD2 u ¼ h with h ¼ ofox2 ofoy1 . Thus, the Stokes equations become a biharmonic equation
M.-J. Lai, P. Wenston / Computers & Fluids 33 (2004) 1047–1073
8 > mD2 u ¼ h > > > > ou > > ¼ g2 < ox ou > > > ¼ g1 > > oy > > : u ¼ b2
1049
in X; on oX; ð1:3Þ on oX; on oX;
where b2 is an anti-derivative of the tangential derivative of u along oX and will be examined in detail later. By a similar calculation, we easily see that the Navier–Stokes equations become the following fourth-order nonlinear equation: 8 o ou o2 u ou o2 u > 2 > > u mD > > oy oy oxoy ox oy 2 > > > > > o ou o2 u ou o2 u > > > ¼ h in X; > < ox oy ox2 ox oxoy ð1:4Þ ou > > ¼ g2 on oX; > > ox > > > > ou > > on oX; ¼ g1 > > > oy > : on oX: u ¼ b2 Let H 2 ðXÞ be the usual Sobolev space and H02 ðXÞ be the subspace of H 2 ðXÞ of functions whose derivatives of order less than or equal to one all vanish on the boundary oX. Define the bilinear form a2 ðu; wÞ and trilinear form qðh; u; wÞ by Z a2 ðu; wÞ ¼ Duðx; yÞDwðx; yÞ dxdy; X Z ouðx; yÞ owðx; yÞ ouðx; yÞ owðx; yÞ Dhðx; yÞ qðh; u; wÞ ¼ dxdy ox oy oy ox X and denote the L2 ðXÞ inner product by Z hh; wi ¼ hðx; yÞwðx; yÞ dxdy: X
We that say u 2 H 2 ðXÞ lowing: 8 ma2 ðu; wÞ ¼ hh; wi > > > > ou > > < ¼ g2 ox ou > > > ¼ g1 > > > : oy u ¼ b2
is a weak solution of the Stokes equations (1.3) if u satisfies the fol8w 2 H02 ðXÞ; on oX; on oX; on oX:
Similarly, a function u 2 H 2 ðXÞ is a weak solution of the Navier–Stokes equations (1.4) if u satisfies
1050
M.-J. Lai, P. Wenston / Computers & Fluids 33 (2004) 1047–1073
8 ma2 ðu; wÞ þ qðu; u; wÞ ¼ hh; wi > > > > > ou > > ¼ g2 < ox ou > > > ¼ g1 > > oy > > : u ¼ b2
8w 2 H02 ðXÞ; on oX; on oX; on oX:
Such weak formulations are referred as the stream function formulation of the Stokes and Navier–Stokes equations, respectively. It is known that the weak solution for StokesÕ equations exists and is unique for any m > 0. For the Navier–Stokes equations, such a weak solution exists for any m > 0, and is unique when m is sufficiently large (see, e.g., [13]). There are two major advantages using the stream function formulation over the traditional velocity–pressure formulation and vorticity–stream function formulation. Indeed, with the stream function formulation, we need to approximate only one stream function. Otherwise, we need to approximate two components of the velocity and one pressure function if the velocity–pressure formulation is used or one vorticity and one stream function if the vorticity–stream function formulation is used. In addition, the stream function formulation eliminates the pressure function which does not have appropriate boundary condition. In the vorticity–stream function formulation, the vorticity function does not have appropriate boundary condition. With bivariate spline functions of higher degrees, we are able to approximate stream functions very well. Thus, in this paper, we will use bivariate splines to approximate the stream function of the Stokes and Navier– Stokes equations. Next we discuss the computation of the pressure functions. By taking divergence, div, of the equations in (1.1) and (1.2), we can easily see that the pressure functions p of the Stokes and Navier–Stokes equations satisfy the following Poisson equations with nonhomogeneous Neumann boundary conditions involving the stream functions 8 in X < Dp ¼ divðfÞ ð1:5Þ op : ¼ mDðn curl uÞ þ n f on oX on for the Stokes equation and 8 < Dp ¼ divðfÞ þ div½ðcurl u rÞcurl ðuÞ op : ¼ mDðn curl uÞ þ n f þ n ½ðcurl u rÞðcurl uÞ on
in X; on oX
ð1:6Þ
for the Navier–Stokes equations. We will use these boundary conditions to compute p after we find the spline approximations of the stream function u. The paper is organized as follows: in Section 2, we first introduce bivariate spline spaces and in particular, the B-form representation of spline functions. The constrained minimization problems that will be used to generate approximate solutions of the stokes and Navier–Stoke equations lead, by way of Lagrange multipliers, to linear systems of the form T F L A k ¼ ; G 0 L c
M.-J. Lai, P. Wenston / Computers & Fluids 33 (2004) 1047–1073
1051
with A singular. We then introduce an iterative method to solve the above system and discuss its convergence. We continue in Section 3 with a discussion of the spline solution of the 2D Stokes equations, an equivalent biharmonic equation. In Section 4, we discuss the spline approximations of the 2D Navier–Stokes equations and prove the convergence of two methods for solving the discrete nonlinear equations. In Section 5, we present numerical results for standard benchmark flows such as the driven cavity flows, backward-facing step flows, and flows around circular obstacles. Also, we also present some other fluid flow experiments.
2. Preliminaries 2.1. Bivariate spline functions Given a bounded polygonal domain X 2 R2 , let M be a triangulation of X. Let d P 1 and r P 1 be two fixed integers. We introduce the spline spaces Sdr ðMÞ ¼ fs 2 C r ðXÞ; sjt 2 Pd 8t 2 Mg; where Pd denotes the space of bivariate polynomials of total degree d. In this paper, the B-form representation of splines on triangulations will be used (cf. [9] or [4]). Let T ¼ hv1 ; v2 ; v3 i be a nondegenerate triangle with vi ¼ ðxi ; yi Þ, i ¼ 1; 2; 3. It is well known that every point v ¼ ðx; yÞ can be written uniquely in the form v ¼ k1 v1 þ k2 v2 þ k3 v3
ð2:1Þ
k1 þ k2 þ k3 ¼ 1;
ð2:2Þ
with where k1 , k2 , and k3 are called the barycentric coordinates of the point v ¼ ðx; yÞ relative to the triangle T . Moreover each ki is a linear polynomial in x; y. Let Bdijk ðvÞ ¼
d! i j k kkk; i!j!k! 1 2 3
i þ j þ k ¼ d:
They are called the Bernstein–Bezier polynomials of degree d. In fact, the set Bd ¼ fBdijk ðx; y; zÞ; i þ j þ k ¼ dg is a basis for the space of polynomials Pd . As a consequence any polynomial p of degree d can be written uniquely in terms of Bdijk Õs, i.e., X p¼ cijk Bdijk : ð2:3Þ iþjþk¼d
The representation (2.3) for polynomials is referred to as the B-form with respect to T . Let iv1 þ jv2 þ kv3 Dd;T ¼ nijk ¼ ; i þ j þ k ¼ d; T 2 M ð2:4Þ d be a set of the domain points of degree d over triangulation M. For each spline function s 2 Sdr ðMÞ, since s restricted to each triangle T 2 M is a polynomial of degree d, we may write
1052
M.-J. Lai, P. Wenston / Computers & Fluids 33 (2004) 1047–1073
sjT ¼
X
cTijk Bdijk ;
T 2 M:
iþjþk¼d
Such a representation is called the B-form representation of the spline function s (cf. [4]). We denote by c :¼ fcTijk ; i þ j þ k ¼ d; T 2 Mg the B-coefficient vector of s. To evaluate a polynomial in B-form, there is the so-called de Casteljau algorithm which we now P ð0Þ describe. For p ¼ iþjþk¼d cijk Bdijk , let us write cijk ¼: cijk ðkÞ with k ¼ ðk1 ; k2 ; k3 Þ being the barycentric coordinates of v ¼ ðx; yÞ with respective to T and define for a positive integer r P 1 ðrÞ
ðr1Þ
ðr1Þ
ðr1Þ
cijk ðkÞ ¼ k1 ciþ1;j;k ðkÞ þ k2 ci;jþ1;k ðkÞ þ k3 ci;j;kþl ðkÞ: We have then X p¼
ðrÞ
cijk ðkÞBdr ijk ;
0 6 r 6 d:
iþjþk¼dr
In particular, for r ¼ d, we have ðdÞ
p ¼ c0;0;0 ðkÞ which is the value of p at v ¼ ðx; yÞ whose barycentric coordinates are k ¼ ðk1 ; k2 ; k3 Þ with k1 þ k2 þ k3 ¼ 1. We next discuss how to take derivatives of polynomials in B-form. We start with formulas for the directional derivatives of p in a direction defined by a vector u. We have X ð1Þ cijk ðaÞBd1 Du p ¼ d ijk iþjþk¼d1
with a ¼ ða1 ; a2 ; a3 Þ the T -coordinates of u; that is u ¼ a1 v1 þ a2 v2 þ a3 v3 with a1 þ a2 þ a3 ¼ 1: Note that if u ¼ v1 v2 is the direction vector, the T -coordinates of u are (1, )1, 0). In general, we have Dmu pðvÞ ¼
d! ðd mÞ!
X
ðmÞ
cijk ðaÞBdm ijkl ðvÞ:
ð2:5Þ
iþjþk¼dm
Note that for arbitrary direction vector u ¼ ðu1 ; u2 Þ 2 R2 Du p ¼ u1
o o p þ u2 p: ox oy
For a triangle T ¼ hv1 ; v2 ; v3 i, if we let u ¼ v2 v1 ¼ ðx2 x1 ; y2 y1 Þ, v ¼ v3 v1 ¼ ðx3 x1 ; y3 y1 Þ, it follows that:
M.-J. Lai, P. Wenston / Computers & Fluids 33 (2004) 1047–1073
ðy3 y1 Þ ðy2 y1 Þ Du p Dv p; 2AT 2AT ðx2 x1 Þ ðx3 x1 Þ Dv p Du p; Dy p ¼ 2AT 2AT
1053
Dx p ¼
ð2:6Þ
where AT is the area of T . There are precise formulas for the integrals and inner products of polynomials in B-form (cf. [5]). Lemma 2.1. Let p be a polynomial of degree d with B-coefficients cijk , i þ j þ k ¼ d on a triangle T . Then Z X AT pðx; yÞ dx dy ¼ cijk ; d þ 2 iþjþk¼d T 2 where AT ¼ area of T . Lemma 2.2. Let q be another polynomial with B-coefficients dijk , i þ j þ k ¼ d, the inner product of p and q over t is given by Z X AT iþr jþs kþt pðx; yÞqðx; yÞ dx dy ¼ cijk drst : ð2:7Þ i j k 2d 2d þ 2 t iþjþk¼d rþsþt¼d d 2 We next discuss the smoothness conditions for a spline function s in Sdr ðMÞ. These are wellknown conditions on the coefficients of s that will assure that s has certain global smoothness properties. Theorem 2.3. Let t ¼ hv1 ; v2 ; v3 i and t0 ¼ hv1 ; v2 ; v4 i be two triangles with common edge hv1 ; v2 i. Then s is of class C r on t [ t0 if and only if X 0 ctiþl;jþm;j Bml;m;j ðv4 Þ; m ¼ 0; . . . ; r; i þ j ¼ d m: ctijm ¼ lþmþj¼m
For a proof, see [9] or [4]. This theorem guarantees the existence of a matrix H such that s is in C ðXÞ if and only if r
H c ¼ 0; where c encodes the B-coefficients of s. 2.2. A matrix iterative method The constrained minimization problems that will be used to generate approximate solutions of the stokes and Navier–Stoke equations will lead, by way of Lagrange multipliers, to linear systems of the form (cf. [10] and [19])
1054
M.-J. Lai, P. Wenston / Computers & Fluids 33 (2004) 1047–1073
LT 0
F A k ¼ ; G L c
ð2:8Þ
with A singular and appropriate matrices L and vectors F and G. Computing a least squares solution of the above system would allow us to compute c when c is unique. However, when the size of a system is very large, any least squares method will be too costly. We present here an iterative method for solving the above system which involves matrices of smaller size. Consider the following sequence of problems: T F A kðkþ1Þ L ¼ ð2:9Þ G kðkÞ I L cðkþ1Þ for k ¼ 0; 1; . . ., with an initial guess kð0Þ , e.g., kð0Þ ¼ 0, and I the identity matrix. Assume that the size of L is m n and the size of A is n n with m < n. Note that (2.9) reads Acðkþ1Þ þ LT kðkþ1Þ ¼ F ; Lcðkþ1Þ kðkþ1Þ ¼ G kðkÞ :
ð2:10Þ
Multiplying on the left the second equation in (2.10) by LT , we get LT Lcðkþ1Þ LT kðkþ1Þ ¼ LT G LT kðkÞ or LT kðkþ1Þ ¼ 1 LT Lcðkþ1Þ 1 LT G þ LT kðkÞ . We substitute this into the first equation in (2.10) to get 1 T 1 A þ L L cðkþ1Þ ¼ F þ LT G LT kðkÞ : ð2:11Þ Using again the first equation in (2.10), i.e., AcðkÞ ¼ F LT kðkÞ to replace F in (2.11), we have 1 T 1 A þ L L cðkþ1Þ ¼ AcðkÞ þ LT G for all k P 1. This leads to the following. Algorithm 2.4. Fix > 0. Given an initial guess kð0Þ 2 ImðLÞ, e.g., kð0Þ ¼ 0, we define cð1Þ by 1 1 T 1 T ð1Þ T ð0Þ F þ L GL k c ¼ Aþ L L and iteratively define 1 1 T 1 T ðkþ1Þ ðkÞ ¼ Aþ L L Ac þ L G c
ð2:12Þ
for k ¼ 1; 2; . . ., where ImðLÞ is the range of L. This algorithm was briefly discussed in [12]. A convergence result was pointed out. That is, kc ckþ1 k 6 Ckc cðkÞ k for some constant C > 0. According to this result for convergence one clearly needs to choose small enough to guarantee that C < 1. Here we present the following convergence result (cf. [2] for a proof) which ensures the convergence for any > 0.
M.-J. Lai, P. Wenston / Computers & Fluids 33 (2004) 1047–1073
1055
Theorem 2.5. Let A ¼ As þ Aa where As ¼ 12ðA þ AT Þ is the symmetric part of A and Aa ¼ 12ðA AT Þ is the anti-symmetric part of A. Furthermore, suppose that As is positive definite with respect to L, that is, xT As x ¼ 0 and Lx ¼ 0 imply that x ¼ 0. Then for any > 0, the matrix 1 A þ LT L is invertible. That is, the Algorithm 2.4 is well defined. Furthermore, there exists positive constants CðÞ > 1 and cðÞ < 1 depending on but independent of k such that kc cðkþ1Þ k 6 CðÞðcðÞÞ
kþ1
for k P 1. Our numerical experiments show that it gives the same result as the least squares method. Since any least squares method is computationally expensive, we recommend the iterative method for large problems. 3. Spline approximations of the Stokes equations In this section, we consider spline approximations of the 2D Stokes equations in the stream function formulation. That is, we consider the following biharmonic equation with nonhomogeneous boundary conditions g1 and g2 : 8 ma2 ðu; wÞ ¼ hh; wi 8w 2 H02 ðXÞ; > > > > ou > > < on oX; ¼ g2 ox ð3:1Þ ou > > on oX; ¼ g1 > > > > : oy u ¼ b2 on oX; where b2 is an anti-derivative of the tangential derivative of u along oX which can be computed based on g1 and g2 . For simplicity, let b1 be the normal derivative of u along oX. Note that b1 can be easily computed using g1 and g2 . It is well known that the weak solution of (3.1) exists and is unique (cf. [13]). We are going to approximate the weak solution of (3.1) by using bivariate spline functions. Let M be a triangulation of X if X is a polygonal domain. Otherwise, we choose vertices v1 ; . . . ; vn on oX which include all the corner points of oX. Instead of the piecewise linear interpolation of the vertices, we use a C 0 piecewise quadratic spline curve soX , which interpolates the vertices, to approximate the boundary oX. Clearly, we may assume that soX is a much better approximation of oX than the piecewise linear interpolant based on an appropriate choice of v1 ; . . . ; vn . Adding interior vertices vnþ1 ; . . . ; vN inside X, we triangulate X using the vertices v1 ; . . . ; vN . For convenience, we let M be the triangulation consisting of all interior triangles and boundary triangles with piecewise quadratic edges. Let d and r be two positive integer with d > r and let S Sd0 ðMÞ be the space of spline functions which are C r inside of X. Note that S H 2 ðXÞ. Recall from Section 2.1 that s 2 Sd0 ðMÞ can be given by
1056
M.-J. Lai, P. Wenston / Computers & Fluids 33 (2004) 1047–1073
X
sðx; yÞ ¼
cti;j;k Bti;j;k ðx; yÞ; 8ðx; yÞ 2 t 2 M:
iþjþk¼d
ðcti;j;k ; i
þ j þ k ¼ d; t 2 MÞ be the B-coefficient vector of s. Recall from Section 2.1 that Let c ¼ there is a matrix H such that if s 2 S, then H c ¼ 0: Although in general we may not be able to find a spline function s 2 S satisfying the boundary conditions exactly, we can find a spline Sb 2 S approximating the boundary conditions. If the the boundary conditions are continuous, the interpolation is the obvious choice. Assuming that v1 ; . . . ; vb are arranged in the counterclockwise direction, Sb interpolates b2 at d þ 1 distinct points on any straight boundary edge hvi ; viþ1 i and at 2d þ 1 distinct points on any quadratic boundary b curved edge vi ^ viþ1 . Also oS interpolates b1 at d distinct points on any straight boundary edge oni hvi ; viþ1 i and at 2d 1 distinct points on any quadratic boundary curved edge vi ^ viþ1 , where ni denotes the outer normal direction of boundary edge hvi ; viþ1 i or vi ^ viþ1 . If the boundary conditions are not continuous, we can still use interpolation except that the interpolation values now come from a conveniently chosen smooth approximation of the boundary conditions. Since the values of any spline and its normal derivative at any of the interpolation points are linear combinations of the B-coefficients of the spline, we may use Bc ¼ b
ð3:2Þ
to express the approximate boundary conditions. Let S0 ¼ S \ H02 ðXÞ. Our spline approximation of the biharmonic equation (3.1) is to find sh;b 2 S satisfying the approximate boundary condition (3.2) such that ma2 ðsh;b ; sÞ ¼ hh; si; 8s 2 S0 :
ð3:3Þ
Using the Lax–Milgram Theorem, problem (3.3) has an unique solution sh;b in S. Based on the theory of elliptic equations (cf. [6]), solving problem (3.3) is equivalent to finding c which minimizes the following energy functional m EðsÞ ¼ a2 ðs; sÞ hh; si 2 subject to s 2 S satisfying the approximation conditions. In terms of B-coefficients, the above energy functional can be approximated by m ~ EðcÞ :¼ cT Ac cT Mh 2 subject to Hc ¼ 0 and Bc ¼ b, where A ¼ diagðAt ; t 2 MÞ is a bending matrix with blocks Z DBti;j;k DBtl;m;n dx dy ; t 2 M; At ¼ iþjþk¼d lþmþn¼d
t
M ¼ diagðMt ; t 2 MÞ is a mass matrix with blocks Z t t Bi;j;k Bl;m;n dx dy ; t2M Mt ¼ t
iþjþk¼d lþmþn¼d
M.-J. Lai, P. Wenston / Computers & Fluids 33 (2004) 1047–1073
1057
and h ¼ ðhtl;m;n ; l þ m þ n ¼ d; t 2 MÞ is a vector with X hl;m;n Btl;m;n dx dy 2 Sd0 ðMÞ sh ¼ lþmþn¼d
the spline interpolation of h over the domain points fntl;m;n ; l þ m þ n ¼ d; t 2 Mg. We note that the mass and bending matrices, M and A, are diagonal block matrices whose blocks can be assembled trivially and in parallel. The matrix H is sparse with each row involving at most ðr þ 1Þðr þ 2Þ=2 þ 1 nonzero elements. Also, the matrix B is very sparse with each row containing at most d þ 1 nonzero elements. Following the Lagrange multiplier method, we let m Lðc; a; bÞ ¼ cT Ac þ aT H c þ bT Bc cT Mh 2 and compute local minimizers. It follows that we need to solve the following: 2 T 32 3 2 3 a Mh H BT mA 4 0 0 H 54 b 5 ¼ 4 0 5: c b 0 0 B We note that the existence and uniqueness of sh;b implies that there exists a unique solution c satisfying the above linear system. Clearly, A is positive definite with respect to ½H; B. Thus, our iterative method introduced in Section 2.2 can be applied to this system and the iterative solutions converge. Next we use the approach presented above to compute the pressure function assuming that we have found the approximation of the stream function u. That is, we need to solve the Poisson problem with Neumann boundary conditions (1.5). Note that we are seeking the pressure in Z 2 2 L0 ðXÞ ¼ p 2 L ðXÞ; p ¼ 0 : X
Again we can use the spline space S as an approximating space. For a spline approximation sp of p with B-coefficient vector cp , its integral over a triangle is simply a weighted sum of its Bcoefficients (cf. Lemma 2.1), hence the spline approximation of the zero mean value condition, R p ¼ 0, can be written X Uc ¼ 0 with U a vector of size N . The pressure is also the solution of a minimization problem, namely that of minimizing the functional Z Z Z 1 2 jDvj ðhÞv ðf nv þ mDðcurl u nÞvÞ P ðvÞ ¼ 2 X X oX over L20 ðXÞ. We now write P in terms of cp . Let Gp interpolate f n þ mðDcurl uÞ n on the boundary. The later is computed by using the spline approximation sh;b of the solution u of the Stokes equations. Let R be the matrix which restricts the B-coefficients c of a spline over X to the B-coefficients of the spline restricted to oX. We have
1058
M.-J. Lai, P. Wenston / Computers & Fluids 33 (2004) 1047–1073
1 P ðcp Þ ¼ ðcp ÞT Kcp hT Mcp GTp Mb Rcp ; 2 where K is the stiffness matrix which is a diagonal block matrix similar to A and Mb is the mass matrix for the bivariate splines of degree d over the edges of oX. Following the Lagrange multiplier method, we need to solve the following linear system: 3 2 32 3 2 k1 U HT K Mh þ RT Mb RGp 5: 40 0 0 H 54 k2 5 ¼ 4 c 0 0 0 U p Again, it is easy to see that K is positive definite with respect to ½H ; U . The iterative method described in Section 2.2 converges.
4. Spline approximations of the Navier–Stokes equations In this section, we consider the bivariate spline approximation of the 2D Navier–Stokes equations in the stream function formulation. Let u 2 H 2 ðXÞ be the weak solution of the Navier– Stokes equations 8 ma2 ðu; wÞ þ qðu; u; wÞ ¼ hh; wi 8w 2 H02 ðXÞ; > > > > ou > > < on oX; ¼ g2 ox ð4:1Þ ou > > on oX; ¼ g > 1 > > > : oy u ¼ b2 on oX: It is known that the weak solution u exists and is unique if m is sufficiently large or h is sufficiently small (cf. [13]). We shall use the approach discussed in the previous section to derive spline approximations of the stream function u. With the same spline space S and the same notation as in the previous section, we have a2 ðu; wÞ ¼ cT Ad; where c and d stand for the B-coefficient vectors of spline functions u; w 2 S. The difference with the previous section is the presence of the nonlinear term, the trilinear form q defined by Z ouðx; yÞ owðx; yÞ ouðx; yÞ owðx; yÞ Dhðx; yÞ dxdy: ð4:2Þ qðh; u; wÞ ¼ ox oy oy ox X Let us first get acquainted with the trilinear form (4.2). Let e encode the B-coefficient vector of the spline function h in S. We can, by tedious computations, show that there is a matrix QðeÞ, which is a linear function of e, such that qðh; u; wÞ ¼ dT QðeÞc: There are several useful properties of QðeÞ given below. Lemma 4.1. QðeÞ is anti-symmetric, i.e., QðeÞT ¼ QðeÞ.
ð4:3Þ
M.-J. Lai, P. Wenston / Computers & Fluids 33 (2004) 1047–1073
1059
Proof. We first note that qðh; u; wÞ ¼ qðh; w; uÞ. In terms of B-coefficient vectors, we have dT QðeÞc ¼ cT QðeÞd: T
T
Thus, it follows from dT QðeÞ c ¼ cT QðeÞd that dT QðeÞ c ¼ dT QðeÞc for all c and d. Hence, QðeÞT ¼ QðeÞ. h Consequently, we have Lemma 4.2. Fix e. dT QðeÞd ¼ 0 for any d satisfying H d ¼ 0. Next we will need the following. Lemma 4.3. There exists a constant C1 such that jdT QðeÞcj 6 C1 kck kdk kek:
Proof. The Markov inequality and a tedious computation yields the above estimate. h The coefficient vector c of the spline approximation of the weak solution satisfies Hc ¼ 0 and Bc ¼ b and mcT Ad þ cT QðcÞd ¼ hT Md
ð4:4Þ
for all coefficient vectors d of splines in S0 . The existence of c can be shown by using the same argument for the existence of the weak solution / satisfying (4.1) (cf. [13,16]). We are mainly interested in computing c. We now observe that if c and k are chosen so that 8 < mAc QðcÞc þ ½ H T BT k ¼ Mh ð4:5Þ Hc ¼ 0 : Bc ¼ b then mcT Ad þ cT QðcÞd þ kT
H d hT Md ¼ 0 B
for all vectors d. Since both H d ¼ 0 and Bd ¼ 0 for all coefficient vectors d of splines in S0 (4.4) follows. We next derive two methods to linearize the nonlinear equations (4.5), using a simple iteration algorithm and NewtonÕs method. Algorithm 4.4 (A simple iteration algorithm). Let ðcð0Þ be the solution of the linear problem (i.e. the associated Stokes equations) and for n ¼ 0; 1; . . ., define ðcðnþ1Þ ; kðnþ1Þ Þ as the solution of
1060
M.-J. Lai, P. Wenston / Computers & Fluids 33 (2004) 1047–1073
mAcðnþ1Þ QðcðnÞ Þcðnþ1Þ þ ½ H T
BT kðnþ1Þ ¼ Mh;
H cðnþ1Þ ¼ 0;
ð4:6Þ
Bcðnþ1Þ ¼ b: We use Algorithm 2.4 to find the new iterates cðnþ1Þ of Algorithm 4.4. Note that A is symmetric and positive definite with respect to H and B and that QðcðnÞ Þ is anti-symmetric by Lemma 4.1. Theorem 2.5 ensures that the matrix iterative method converges for the linear system in (4.6). To apply NewtonÕs method to solve (4.5), we consider the mapping C defined by C : ðc; kÞ 7! ðmAc QðcÞc þ ½ H T
BT k Mh; H c; Bc bÞ
and seek to solve Cðc; kÞ ¼ 0. We write XðnÞ ¼ ðcðnÞ ; kðnÞ Þ and let Xð0Þ be the solution of the linear Stokes problem. Then we define Xðnþ1Þ by the equation C0 ðXðnÞ ÞðXðnþ1Þ XðnÞ Þ ¼ CðXðnÞ Þ: Thanks to the bilinearity of the mapping ðc; dÞ ! QðcÞd, the above equation leads to mAcðnþ1Þ QðcðnÞ Þcðnþ1Þ Qðcðnþ1Þ cðnÞ ÞcðnÞ þ ½ H T
BT kðnþ1Þ ¼ Mh;
H cðnþ1Þ ¼ 0; ðnþ1Þ
Bc
ð4:7Þ
¼ b:
e ðnÞ Þðcðnþ1Þ cðnÞ Þ with Qðc e ðnÞ ÞcðnÞ ¼ It can be shown that Qðcðnþ1Þ cðnÞ ÞcðnÞ can be written Qðc e ðnÞ Þ. We have our second algorithm. QðcðnÞ ÞcðnÞ for some matrix Qðc Algorithm 4.5 (NewtonÕs method). Let ðcð0Þ Þ be the solution of the Stokes equations and for n ¼ 0; 1; 2; . . . define ðcðnþ1Þ ; kðnþ1Þ Þ to be the solution of e ðnÞ Þcðnþ1Þ þ ½ H T mAcðnþ1Þ QðcðnÞ Þcðnþ1Þ Qðc
e ðnÞ ÞcðnÞ ; BT kðnþ1Þ ¼ Mh Qðc
H cðnþ1Þ ¼ 0; ðnþ1Þ
Bc
ð4:8Þ
¼ b:
For the linear system (4.8), Algorithm 2.4 converges by Theorem 2.5 because QðcðnÞ Þ is antie ðnÞ Þ is positive definite with respect to ½H ; B for m sufficiently large. Indeed, symmetric and mA þ Qðc since A is positive definite with respect to H , we have xT Ax P a0 kxk2 for any x satisfying H x ¼ 0. It follows from Lemma 4.3, e ðnÞ Þx P ma0 kxk2 C1 kcðnÞ k kxk2 : mxT Ax þ xT Qðc By Lemma 4.6, the sequence kcðnÞ k is bounded and hence e ðnÞ Þx P ðma0 C1 C2 Þkxk2 : mxT Ax þ xT Qðc e ðnÞ Þ is positive definite with respect to ½H ; B. Therefore, for m sufficiently large, mA þ Qðc
M.-J. Lai, P. Wenston / Computers & Fluids 33 (2004) 1047–1073
1061
Lemma 4.6. Suppose that cðnþ1Þ is the solution of (4.8). Then kcðnþ1Þ k 6 C for a positive constant C independent of n. Proof. Let c be the solution of (4.5). Let eðnþ1Þ ¼ c cðnþ1Þ and sðnþ1Þ ¼ k kðnþ1Þ . Then we have mAeðnþ1Þ þ QðcðnÞ Þcðnþ1Þ þ Qðcðnþ1Þ ÞcðnÞ QðcÞc þ ½H ; BT sðnþ1Þ ¼ QðcðnÞ ÞcðnÞ ; which gives mAeðnþ1Þ QðcðnÞ Þeðnþ1Þ QðeðnÞ Þc þ Qðcðnþ1Þ ÞcðnÞ þ ½H ; BT sðnþ1Þ ¼ QðcðnÞ ÞcðnÞ or mAeðnþ1Þ QðcðnÞ Þeðnþ1Þ QðeðnÞ Þc Qðeðnþ1Þ ÞcðnÞ þ QðeðnÞ ÞcðnÞ þ ½H ; BT sðnþ1Þ ¼ 0: We multiply this equation on the left by ðeðnþ1Þ ÞT and get mðeðnþ1Þ ÞT Aeðnþ1Þ ¼ ðeðnþ1Þ ÞT QðeðnÞ ÞeðnÞ þ ðeðnþ1Þ ÞT Qðeðnþ1Þ ÞcðnÞ : It follows that
ma0 keðnþ1Þ k 6 C1 keðnÞ k2 þ keðnþ1Þ kkcðnÞ k 6 C1 keðnÞ k2 þ keðnþ1Þ kðkeðnÞ k þ kckÞ :
1 kck We now use an induction to prove the desired result. Assume that keðnÞ k 6 K with K 6 12ma0 C . C1 This is always possible if m is large enough (cf. [14]). Then we have
ðma0 C1 ðK þ kckÞÞkeðnþ1Þ k 6 C1 keðnÞ k2 or keðnþ1Þ k 6
C1 K 2 6 K: ma0 C1 ðK þ kckÞ
This completes the proof with C2 ¼ K þ kck. h 5. Numerical simulation of various fluid flows In this section we first present our numerical experiments with three benchmark flow examples: the driven cavity flow, the backward facing step flow, and the flow around a circular object. We show the streamlines, cavity and pressure of the driven cavity flows, the streamlines for other two kinds of flows. Then we present some other examples of fluid flows. All the figures are shown in the Cartesian coordinates with the horizontal direction being the x-axis and the vertical direction being the y-axis. 5.1. Computational experiments on driven cavity flows Example 5.1. Let us consider the standard cavity flow over unit square domain ½0; 1 ½0; 1. The boundary conditions are u ¼ ðu1 ; u2 Þ ¼ ð0; 0Þ for all four line boundary pieces except for
1062
M.-J. Lai, P. Wenston / Computers & Fluids 33 (2004) 1047–1073
u1 ðx; 1Þ ¼ 1 when 0 6 x 6 1. With Reynolds numbers are 10, 100, 1000, 5000, 10000 and 20,000, the stream lines of the cavity flow are shown in Figs. 2–5. We use several different triangulations for the computation. Two typical triangulations are shown in Fig. 1. We use the splines of degree 8 or higher and smoothness 2. Our main reason to use splines of higher order is to achieve the best approximation property (cf. [15]). All the computation is done using a PC with MS Window 2000. Certainly, our method can compute the stream functions for much higher Reynolds numbers. The stream lines, vorticity, and pressure are similar based on different triangulations, especially when Reynolds numbers are small. We use the triangulation (the right one in Fig. 1) to capture the eddies at the corners.
1
1
0.9
0.9
0.8
0.8
0.7
0.7
0.6
0.6
0.5
0.5
0.4
0.4
0.3
0.3
0.2
0.2
0.1
0.1
0
0 0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
Fig. 1. Two typical triangulations of the square domain.
Fig. 2. The stream lines of driven cavity flow (Reynolds number ¼ 10 and 100).
0.9
1
M.-J. Lai, P. Wenston / Computers & Fluids 33 (2004) 1047–1073
1063
Fig. 3. The stream lines of driven cavity flow (Reynolds number ¼ 1000 and 5000).
Fig. 4. The stream lines of driven cavity flow (Reynolds number ¼ 10,000 and 20,000).
It can be seen that the stream lines in Figs. 2 and 3 are similar to those computed by using the finite difference method given in [11]. The stream lines in Fig. 4 are very similar to the one in [7] which is computed based on time dependent Navier–Stokes equations solved by a finite difference method. Also, the vorticity of the flows wðx; yÞ ¼
o2 o2 /ðx; yÞ þ /ðx; yÞ; ox2 oy 2
1064
M.-J. Lai, P. Wenston / Computers & Fluids 33 (2004) 1047–1073
are shown as in Figs. 5–7. They are very similar to those in [11]. The vorticity in Fig. 7 with Reynolds number ¼ 10,000 is similar to the one in [8] which is computed based on time dependent Navier–Stokes equations and a fourth order finite difference method. Next we compute the pressure of the flows which is normalized such that the integral of the pressure over the domain is zero. As we pointed out before, the pressure is computed by solving the Neumann problem of the standard Poisson equation after the stream function is obtained. Again we set up as a minimization problem which minimizes the stiff matrix subject to the smoothness conditions, boundary conditions, and a normalization condition (cf. [2] for multivariate spline method for numerical solution of partial differential equations). The pressures for flows with different Reynolds numbers are given in Figs. 8–10.
Fig. 5. The contour of the vorticity of the driven cavity flow (Reynolds number ¼ 10 and 100).
Fig. 6. The contour of the vorticity of the driven cavity flow (Reynolds number ¼ 1000, 5000).
M.-J. Lai, P. Wenston / Computers & Fluids 33 (2004) 1047–1073
1065
Fig. 7. The contour of the vorticity of the driven cavity flow (Reynolds number ¼ 10,000 and 20,000).
Fig. 8. The contour of the pressure of the driven cavity flow (Reynolds number ¼ 10 and 100).
5.2. Computational experiments on backward facing step flows Next we consider the well-known backward facing step flows. The inflow condition at the left inlet is the parabolic flow with maximum value 1. The outflow condition at the right outlet is Outflow condition. Both velocity components remain the same in the direction normal to the boundary. That is o2 u¼0 on2
and
o o u¼0 on os
1066
M.-J. Lai, P. Wenston / Computers & Fluids 33 (2004) 1047–1073
Fig. 9. The contour of the pressure of the driven cavity flow (Reynolds number ¼ 1000 and 5000).
Fig. 10. The contour of the pressure of the driven cavity flow (Reynolds number ¼ 10,000 and 20,000).
at the outlet, where s; n denote the tangential and normal direction of the outlet. Figs. 11–14 are some numerical simulations of the backward-facing step flows. For detailed explanation of these flows, see, e.g. [11] and the reference therein. The streamlines of the flows from our numerical experiments match with those in [11]. 5.3. Computational experiments of flows around circular obstacle In addition to the no-slip and outflow boundary conditions in the previous sections, we need the following freeslip conditions for simulating the flows around some obstacles.
M.-J. Lai, P. Wenston / Computers & Fluids 33 (2004) 1047–1073
Fig. 11. Flows over a backward facing step (Reynolds number ¼ 1).
Fig. 12. Flows over a backward facing step (Reynolds number ¼ 100).
Fig. 13. Flows over a backward facing step (Reynolds numbers ¼ 250).
Fig. 14. Flows over a backward facing step (Reynolds number ¼ 500).
Fig. 15. Flow around a disk with Reynolds number ¼ 10.
Fig. 16. Flow around a disk with Reynolds number ¼ 50.
Fig. 17. Flow around a disk with Reynolds number ¼ 100.
1067
1068
M.-J. Lai, P. Wenston / Computers & Fluids 33 (2004) 1047–1073
Fig. 18. Triangulation of the triangular domain and the stream lines of the triangular driven cavity flow (Reynolds number ¼ 100).
Fig. 19. Contour of the vorticity and pressure of the driven cavity flow (Reynolds number ¼ 100).
Fig. 20. The stream lines of the triangular driven cavity flow (Reynolds number ¼ 1000).
M.-J. Lai, P. Wenston / Computers & Fluids 33 (2004) 1047–1073
1069
Fig. 21. Contour of the vorticity and pressure of the triangular driven cavity flow (Reynolds number ¼ 1000).
Fig. 22. The stream lines of the triangular driven cavity flow (Reynolds number ¼ 5000).
Fig. 23. Contour of the vorticity and pressure of the triangular driven cavity flow (Reynolds number ¼ 5000).
1070
M.-J. Lai, P. Wenston / Computers & Fluids 33 (2004) 1047–1073
Fig. 24. The flows around a boxy car model and a streamlined car model.
Fig. 25. The pressures around a boxy car model and a streamlined car model.
M.-J. Lai, P. Wenston / Computers & Fluids 33 (2004) 1047–1073
1071
Fig. 26. Flows through sudden contraction and enlargement.
Freeslip conditions. No fluid penetrates the boundary, i.e. o u¼0 on
and
o o u¼0 on os
along a part of the boundary oX. Here are some numerical simulations of flows around a circular disk for various Reynolds numbers (Figs. 11–17). The in-flows are at a constant velocity at the left vertical boundary, i.e., u ¼ 1 and v ¼ 0. These flows can be compared with the ones in [11,20]. 5.4. Computational experiments on other interested flows We first show driven cavity flows over a triangular domain. Example 5.2. Let us consider the cavity flow over triangular domain. We uniformly refine the triangle as shown in Fig. 18 and use bivariate splines of degree 8 and smoothness 2. The boundary conditions are u ¼ ðu1 ; u2 Þ ¼ ð0; 0Þ for all three line boundary pieces except for u1 ðx; 1Þ ¼ 1 when 0 6 x 6 1. With Reynolds numbers are 100, 1000, and 5000, the stream lines of the cavity flow are shown in Figs. 18, 20 and 22, the contours of vorticity and pressure are shown in Figs. 19, 21 and 23. Next we present some additional examples which need inflow and outflow conditions. Example 5.3. In this example, we compare with the flows around a boxy car model and streamlined car model as in Fig. 24. The flows are similar to the real simulations as compiled in [17]. The pressures of two flows are computed as shown in Fig. 25. The highest pressure in the front of the boxy car is about 140 while that in the front of the smooth car is about 80. These clearly show that the boxy car has more wind resistance than the smooth car. Example 5.4. In this example, we consider the flow through sudden contraction and enlargement. We compare with the flows with two levels of contraction and enlargement. The flows are constant at the left vertical boundary and are the same for the two flows in Fig. 26, but result two different flow patterns and pressure levels as shown in Fig. 27.
1072
M.-J. Lai, P. Wenston / Computers & Fluids 33 (2004) 1047–1073
Fig. 27. Pressures of the flows through sudden contraction and enlargement.
References [1] Awanou G, Lai MJ. Trivariate spline approximations of the 3D Navier–Stokes equations, submitted. [2] Awanou G, Lai MJ, Wenston P. The multivariate spline method for numerical solutions of partial differential equations and scattered data fitting, submitted. [3] Botella O. On a collocation B-spline method for the solution of the Navier–Stokes equations. Comput Fluids 2002;13:397–420. [4] de Boor C. B-form basics. In: Farin G, editor. Geometric modelling. Philadelphia: SIAM Publication; 1987. p. 131– 48. [5] Chui CK, Lai MJ. Multivariate vertex splines and finite elements. J Approx Theory 1990;60:245–343. [6] Evans L. Partial Differential Equations. Providence: American Mathematical Society; 1998.
M.-J. Lai, P. Wenston / Computers & Fluids 33 (2004) 1047–1073
1073
[7] Weinan E, Liu J-G. Vorticity boundary condition and related issues for finite difference schemes. J Comput Phys 1996;124:368–82. [8] Weinan E, Liu J-G. Finite difference schemes for incompressible flows in the velocity impulse density formulation. J Comput Phys 1997;130:67–76. [9] Farin G. Triangular Bernstein–Bezier patches. Comput Aided Geom Design 1986;3:83–127. [10] Fortin M, Glowinski R. Augmented Lagrangian methods. Amsterdam: Elsevier; 1983. [11] Griebel M, Dornseifer T, Neunhoeffer T. Numerical simulation in fluid dynamics. Philadelphia: SIAM Publications; 1998. [12] Gunzburger M. Finite element methods for viscous incompressible flows. New York: Academic Press; 1989. [13] Giraultt V, Raviart PA. Finite element method for Navier–Stokes equations. Berlin: Springer-Verlag; 1986. [14] Karakashian O. On a Galerkin–Lagrange multiplier method for the stationary Navier–Stokes equations. SIAM J Numer Anal 1982;19:909–23. [15] Lai MJ, Schumaker LL. Approximation power of bivariate splines. Adv Comput Math 1998;9:251–79. [16] Lai MJ, Wenston P. Bivariate spline method for numerical solution of Navier–Stokes equations over polygons in stream function formulation. Numer Methods PDE 2000;16:147–83. [17] Nakayama Y. Visualized flow. Oxford: Pergamon Press; 1978. [18] Serrin J. On the interior regularity of weak solution of Navier–Stokes equations. Arch Rat Mech Anal 1962;9:187– 95. [19] Temam R. Navier–Stokes equations. In: Theory and numerical analysis. Amsterdam: North-Holland Publishing Co.; 1984. [20] Van Dyke M. An album of fluid motion. Stanford, CA: Parabolic Press; 1982.