1
行政院國家科學委員會專題研究計畫成果報告
邊界元素法在黏性流之應用
A boundar y element methods for incompr essible viscous flows
計畫編號：NSC 89-2611-E-002-045
執行期限：89 年 8 月 1 日至 91 年 7 月 31 日
主持人：黃維信 國立台灣大學造船及海洋工程學系
一、中文摘要 本 文 旨 在 以 邊 界 元 素 法 求 解 Navier-Stokes Equations。首先將非線性的 對流項在局部地區給予線性化處理，使得 Navier-Stokes equations 變成局部的線性方 程式。之後，將局部區域的線性方程式轉 化為穩態的橢圓型偏微分方程式，再將此 組 方 程 式 轉 換 為 modified Helmholtz equation 並求此方程式的 Green 函數。利 用此 Green 函數及 Green‘s Identity 將偏微 分方程式轉換為積分方程式，並使用邊界 元素法作為求解的工具。關鍵詞：邊界元素法、黏性流 Abstr act
A linearized Navier-Stokes equation is assumed to approximate the nonlinear steady Navier-Stokes equations locally. In this linearized procedure, the convective terms are simplified but still kept in the
whole equation. Usually, when the
convective terms are included in the equations, the formulation of the boundary element methods will contain a volume integral instead of a full boundary integral. In this proposal, we propose a new formulation to avoid the volume integral.
Keywor ds: Navier-Stokes Equation, boundary element method
二、緣由與目的
The numerical treatment of fluid dynamic problems is of increasing theoretical and practical importance in engineering and science. In the computational fluid dynamics, finite difference and finite element methods are the two most important and broadly used tools in the present stage. In the past decade, the
researchers are devoted themselves to
developing some new techniques different from the finite difference and finite element methods. One of them is the boundary element method. The advantage of the boundary element method over the other numerical
approaches is the dimension of the
computational domain decreasing by one for the linear differential equations. Unfortunately, viscous fluid flow is governed by non-linear equations and the boundary value problem cannot be converted into a boundary integral equation only. When a volume integral is involved in the integral formulation, several techniques were developed to overcome this type of problems in the boundary element method such as dual reciprocal method, Fourier series method, and multiple reciprocal method. Basically, in these methods the volume integral is contained in the integral formulation, and the volume integral is computed directly or hidden after some operation of boundary integrals.
When modeling the Navier-Stokes equation by a boundary integral formulation, the nonlinear convective term is frequently
2 simplified as a linear term. In the early development, the nonlinear convective term was considered as a pseudo-body force, then the linear differential equation might be transform into boundary integral form easily. In such cases, the fundamental solution of The Laplace equation or the Stokes equations has been utilized in which the convective term has not been taken into consideration. In most recent, the linearized convective term is considered as a part of differential operator. One of them utilizes the Oseen equation to simulate the local fluid flows for the Navier-Stokes equation.
The other employs velocity-vorticity
formulation of the conservation equation. Both of them contain volume integrals in their boundary integral equations.
In this paper, we have develop a new formulation of the boundary integral equation without a volume integral for the viscous fluid flows. In this case, the numerical completeness will become clear and simple, and the
numerical technique developed for the
traditional boundary element methods in the past can be used without too much modification.
三、結果與討論
For simplicity, the primitive variable formulation is used to describe the viscous fluid flows. The continuity equation and the
Navier-Stokes equations represent the
conservation laws of mass and momentum, respectively, for viscous incompressible fluid flow confined in a bounded domain. They are written below in an indicial notation form for a Cartesian coordinate system.
0 = ∂ ∂ i i x u (1) and ∂ ∂ + ∂ ∂ ∂ ∂ + ∂ ∂ − = i j j i j i i x u x u x x p Dt Du ν ρ 1 , (2)
where xi is the coordinate, ui the flow velocity
at the ith component, ρ the density, p the pressure, and ν the kinematic viscosity.
When a boundary element method is applied to represent a viscous flow for either steady or unsteady problems, the nonlinear terms of the Navier-Stokes equations have to be locally linearized, and the resulted
equations become diffusion-convective
equations in a sub-domain such as
j i j i j j i j i j i j i x u v x u x u x x p x u U t u ∂ ∂ + ∂ ∂ + ∂ ∂ ∂ ∂ + ∂ ∂ − = ∂ ∂ + ∂ ∂ ν ρ 1 , (3)
where Ujis the averaged velocity over the local
sub-domain and uj= Uj + vj. If the unsteady
term is descretized by a finite difference scheme, and the perturbated convective term is treated as a source term, this set of equations can be expressed by a simple mathematical form, d bu u a ui = j i j + + ∇2 2 , , (4)
where aj= Uj, bu represents the partial effect of
unsteady terms, and d represents the
combination effect of perturbated convective terms, a part of unsteady terms due to the previous time step from the finite difference, and the pressure gradient. In the following, we will condider Eq. (4) with and without the unsteady term, separately, and omit the subscript i for simplicity. At first, if the
unsteady term is neglected in Eq. (4), then
d u a ui = j i j + ∇2 2 , . (5) Let j j i i x a a a x da x f e u k k 2 ) ( − = (6)
3 By substituting the above equation into Eq. (5), a the modified Helmholtz equation can be obtained:
0 2
2 − =
∇ f k f (7)
where k2 =ajaj in a subdomain. This differential equation can furtherly be transferred into a boundary integral equation without a volume integral, if the Green formula and a suitable Green's function are adopted. The Green's function of the modified Helmholtz equation in an unbounded domain satisfies the equation:
) ( 2
2 − =δ r−ξr
∇ G k G x , (8) and is given for the two-dimensional case as
π 2 / ) ( 0 kr K G=− (9)
and for the three-dimensional case as
r e
G=− −kr/4π , (10) where r is the distance from the source point
x
r
to the field point ξr, and K0 is the modified
Bessel function of the second kind. To form the boundary integral representation for Eq. (7) a few steps are needed. Subtract Eq. (7) multiplying a factor G from Eq. (8) multiplying a factor f, integrate it over a subdomain, and transform the volume integral into the surface integral by the Gauss divergence theorem. We have
dS n f G n G f f S ] [ ∂ ∂ − ∂ ∂ =
∫
ε , (11) and∫
∂∂ = S dS n G0 ε (12)where G0 represents the function G when k =
0.
If the unsteady term is kept in Eq. (4), the transformation in Eq. (6) has to be changed slightly such as:
b d x f e u= ajxj ( )− . (13)
After substituting Eq. (13) into Eq. (4), a slightly different equation is obtained:
0 ) (
2 − + =
∇ f ajaj b f (14)
Basically, Eqs. (6) and (14) are of the same type of equations. The coefficient k2³ 0 in Eq. (6), but it is not always true for Eq. (14). Equations (9) and (10) are still the Green's function of Eq. (14), if k2 = aj aj + b.
For Eq. (5), let u = u + daixi 2ajaj , and substitute it into Eq. (11). The Green's formula can be modified such as:
dS n u G u n Ga n G e u j j S x aj j j )] ) [( ) ( ∂ ∂ − + ∂ ∂ =
∫
−ξ ε .For a boundary element method, the above equation can be discretized as a matrix form:
0 = ∂ ∂ + n u G u H
Though Eq. (4) and (5) are very similar, their basic properties are not so close. If the source term d is treated as a constant in a
subdomain, the type of homogeneous solutions of these two equations is different. The homogeneous solutions of these equations are the most natural bases in the subdomain. For a
two-dimensional case, the homogeneous
solution in Eq. (4) includes only exponential functions such as:
y c p a or a x p a a e c u0 0 ( ) ( ( ) ) 2 2 2 2 1 1± − + ± + − = m ,
where c = aj aj + b and p is an arbitray constant.
For Eq. (5), it contains more types of functions such as: 2 2 1 1 ) ) ( ( ) ( 0 0 ) ( 2 2 2 2 1 1 c y a x a c e c u a a p x a or a p y + − + = ± − + ± m + .
To solve a PDE problem uniquely, a preassigned condition on the whole boundary is necessary. In the boundary element method, the boundary data are only specified on a few nodes of the boundary, and the other part of boundary data are interpolated by low-order polynomials. From the above analysis, such
4 low-order polynomials may not match the natural behavior inside the domain of the
Navier-Stokes equations. In some
finite-difference type methods, such as finite volume methods or finite analytical methods,
exponential functions are chosen as
approximating functions. Many of their results indicated that the use of exponential functions maked the numerical scheme stable when comparing with pure polynomials. Physically,
it means that the correct choice of
approximating functions for the problem result in stable solutions.
Up to now, the users of boundary element methods follow finite element methods to use isoparamatric elements to construct their solutions. Usually, finite element methods use low-order polynomials to construct the shape of elements, and so for the physical properties. In this paper, we still approximate the shape of elements by low-order polynomials, but no longer for the physical properties. Both exponential functions and polynomials are combined as basis functions in an element, which is no longer an isoparametric element.
For simplicity, Eq. (5) is chosen as a test case and the basis functions on the boundary are chosen such as:
2 2 1 1 2 2 0 0 c e 1 2 c (a x a y) c u = a x+ a y + − + ,
where c0, c1, and c2 are undetermined
coefficients.
四、參考文獻
1. M.M. Grigoriev and A.V. Fafurin, 'A boundary element method for steady viscous fluid flow using penalty function formulation,' Int. J. Numer. Meth. Fluids, 25, 907-929 (1997).
2. M.M. Grigoriev and G.F. Dargush, 'A poly-region boundary element method for incompressible viscous fluid flows,' Int. J. Numer. Meth. Eng., 46, 1127-1158 (1999). 3. L. Skegret, M. Hribersek and G. Kuhn,
'Computational fluid dynamics by
boundary-domain integral method,' Int. J.
Numer. Meth. Eng., 46, 1291-1311 (1999).
五、成果自評
本研究內容與原計畫相符、並達成預期 目標、研究成果兼具學術及應用價值、適合 在學術期刊發表。
