# Bivariate splines for fluid flows

## Bivariate splines for fluid flows

Computers & Fluids 33 (2004) 1047–1073 www.elsevier.com/locate/compﬂuid Bivariate splines for ﬂuid ﬂows Ming-Jun Lai *, Paul Wenston Department of Ma...

Computers & Fluids 33 (2004) 1047–1073 www.elsevier.com/locate/compﬂuid

Bivariate splines for ﬂuid ﬂows 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-coeﬃcients 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 eﬀective. Numerical simulations of several ﬂuid ﬂows will be included to demonstrate the eﬀectiveness of the bivariate spline method. Ó 2003 Elsevier Ltd. All rights reserved.

1. Introduction One of the major tasks for applied mathematicians is to develop eﬃcient methods for numerical solutions of partial diﬀerential equations. The steady-state Navier–Stokes equations are one of the most important partial diﬀerential 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 eﬃcient methods continue to be developed in order to increase the power of computational ﬂow 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 ﬁnite 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).

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 diﬀerent 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 ﬂexible 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 stiﬀness 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 satisﬁes the divergence-free condition exactly; (6) the matrices that arise are singular which is an important diﬀerence from the classical ﬁnite 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 ﬁrst introduce the stream function formulation. Let X  R2 be a simply connected T domain and u ¼ ðu1 ; u2 Þ be the planar velocity of a ﬂuid ﬂow over X. Also, let p be the pressure T T function, f ¼ ðf1 ; f2 Þ be the external body force of the ﬂuid and g ¼ ðg1 ; g2 Þ be the velocity of the ﬂuid ﬂow 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 ﬁrst. Replacing u by curl u and then diﬀerentiating the ﬁrst equation with respect to y and the second with respect to x, we subtract the ﬁrst 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. Deﬁne 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 satisﬁes 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 satisﬁes

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 suﬃciently 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 ﬁnd the spline approximations of the stream function u. The paper is organized as follows: in Section 2, we ﬁrst 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 ﬂows such as the driven cavity ﬂows, backward-facing step ﬂows, and ﬂows around circular obstacles. Also, we also present some other ﬂuid ﬂow 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 ﬁxed 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–Bezier 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-coeﬃcient 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 deﬁne for a positive integer r P 1 ðrÞ

ðr1Þ

ðr1Þ

ðr1Þ

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ÞBdr ijk ;

0 6 r 6 d:

iþjþk¼dr

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 deﬁned by a vector u. We have X ð1Þ cijk ðaÞBd1 Du p ¼ d ijk iþjþk¼d1

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ÞBdm ijkl ðvÞ:

ð2:5Þ

iþjþk¼dm

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 coeﬃcients 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-coeﬃcients 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 ﬁrst 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 ﬁrst 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 deﬁne cð1Þ by  1   1 T 1 T ð1Þ T ð0Þ F þ L GL k c ¼ Aþ L L   and iteratively deﬁne  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 brieﬂy discussed in [12]. A convergence result was pointed out. That is, kc  ckþ1 k 6 Ckc  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-coeﬃcient 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 ﬁnd a spline function s 2 S satisfying the boundary conditions exactly, we can ﬁnd 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-coeﬃcients 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 ﬁnd 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 ﬁnding 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-coeﬃcients, 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 deﬁnite 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-coeﬃcient vector cp , its integral over a triangle is simply a weighted sum of its Bcoeﬃcients (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-coeﬃcients c of a spline over X to the B-coeﬃcients 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 stiﬀness 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 deﬁnite 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 suﬃciently large or h is suﬃciently 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-coeﬃcient vectors of spline functions u; w 2 S. The diﬀerence with the previous section is the presence of the nonlinear term, the trilinear form q deﬁned 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 ﬁrst get acquainted with the trilinear form (4.2). Let e encode the B-coeﬃcient 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 ﬁrst note that qðh; u; wÞ ¼ qðh; w; uÞ. In terms of B-coeﬃcient 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 coeﬃcient vector c of the spline approximation of the weak solution satisﬁes Hc ¼ 0 and Bc ¼ b and mcT Ad þ cT QðcÞd ¼ hT Md

ð4:4Þ

for all coeﬃcient 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 coeﬃcient 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; . . ., deﬁne ð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 ﬁnd the new iterates cðnþ1Þ of Algorithm 4.4. Note that A is symmetric and positive deﬁnite 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 deﬁned 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 deﬁne 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; . . . deﬁne ð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 deﬁnite with respect to ½H ; B for m suﬃciently large. Indeed, symmetric and mA þ Qðc since A is positive deﬁnite 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 deﬁnite with respect to ½H ; B. Therefore, for m suﬃciently 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 ; BT 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 ; BT 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 ; BT 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 ﬂuid ﬂows In this section we ﬁrst present our numerical experiments with three benchmark ﬂow examples: the driven cavity ﬂow, the backward facing step ﬂow, and the ﬂow around a circular object. We show the streamlines, cavity and pressure of the driven cavity ﬂows, the streamlines for other two kinds of ﬂows. Then we present some other examples of ﬂuid ﬂows. All the ﬁgures 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 ﬂows Example 5.1. Let us consider the standard cavity ﬂow 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 ﬂow are shown in Figs. 2–5. We use several diﬀerent 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 diﬀerent 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 ﬂow (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 ﬂow (Reynolds number ¼ 1000 and 5000).

Fig. 4. The stream lines of driven cavity ﬂow (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 ﬁnite diﬀerence 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 ﬁnite diﬀerence method. Also, the vorticity of the ﬂows 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 ﬁnite diﬀerence method. Next we compute the pressure of the ﬂows 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 stiﬀ matrix subject to the smoothness conditions, boundary conditions, and a normalization condition (cf. [2] for multivariate spline method for numerical solution of partial diﬀerential equations). The pressures for ﬂows with diﬀerent Reynolds numbers are given in Figs. 8–10.

Fig. 5. The contour of the vorticity of the driven cavity ﬂow (Reynolds number ¼ 10 and 100).

Fig. 6. The contour of the vorticity of the driven cavity ﬂow (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 ﬂow (Reynolds number ¼ 10,000 and 20,000).

Fig. 8. The contour of the pressure of the driven cavity ﬂow (Reynolds number ¼ 10 and 100).

5.2. Computational experiments on backward facing step ﬂows Next we consider the well-known backward facing step ﬂows. The inﬂow condition at the left inlet is the parabolic ﬂow with maximum value 1. The outﬂow 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 ﬂow (Reynolds number ¼ 1000 and 5000).

Fig. 10. The contour of the pressure of the driven cavity ﬂow (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 ﬂows. For detailed explanation of these ﬂows, see, e.g. [11] and the reference therein. The streamlines of the ﬂows from our numerical experiments match with those in [11]. 5.3. Computational experiments of ﬂows around circular obstacle In addition to the no-slip and outﬂow boundary conditions in the previous sections, we need the following freeslip conditions for simulating the ﬂows 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 ﬂow (Reynolds number ¼ 100).

Fig. 19. Contour of the vorticity and pressure of the driven cavity ﬂow (Reynolds number ¼ 100).

Fig. 20. The stream lines of the triangular driven cavity ﬂow (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 ﬂow (Reynolds number ¼ 1000).

Fig. 22. The stream lines of the triangular driven cavity ﬂow (Reynolds number ¼ 5000).

Fig. 23. Contour of the vorticity and pressure of the triangular driven cavity ﬂow (Reynolds number ¼ 5000).

1070

M.-J. Lai, P. Wenston / Computers & Fluids 33 (2004) 1047–1073

Fig. 24. The ﬂows 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 ﬂuid 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 ﬂows around a circular disk for various Reynolds numbers (Figs. 11–17). The in-ﬂows are at a constant velocity at the left vertical boundary, i.e., u ¼ 1 and v ¼ 0. These ﬂows can be compared with the ones in [11,20]. 5.4. Computational experiments on other interested ﬂows We ﬁrst show driven cavity ﬂows over a triangular domain. Example 5.2. Let us consider the cavity ﬂow over triangular domain. We uniformly reﬁne 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 ﬂow 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 inﬂow and outﬂow conditions. Example 5.3. In this example, we compare with the ﬂows around a boxy car model and streamlined car model as in Fig. 24. The ﬂows are similar to the real simulations as compiled in [17]. The pressures of two ﬂows 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 ﬂow through sudden contraction and enlargement. We compare with the ﬂows with two levels of contraction and enlargement. The ﬂows are constant at the left vertical boundary and are the same for the two ﬂows in Fig. 26, but result two diﬀerent ﬂow 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 ﬂows 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 diﬀerential equations and scattered data ﬁtting, 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 ﬁnite elements. J Approx Theory 1990;60:245–343. [6] Evans L. Partial Diﬀerential 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 ﬁnite diﬀerence schemes. J Comput Phys 1996;124:368–82. [8] Weinan E, Liu J-G. Finite diﬀerence schemes for incompressible ﬂows in the velocity impulse density formulation. J Comput Phys 1997;130:67–76. [9] Farin G. Triangular Bernstein–Bezier 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, Neunhoeﬀer T. Numerical simulation in ﬂuid dynamics. Philadelphia: SIAM Publications; 1998. [12] Gunzburger M. Finite element methods for viscous incompressible ﬂows. 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 ﬂow. 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 ﬂuid motion. Stanford, CA: Parabolic Press; 1982.