# Stabilized Discretization in Spline Element Method for Solution of Two-Dimensional Navier-Stokes Problems

Authors: Neng Wan; Ke Du; Tao Chen; Sentang Zhang; Gongnan Xie
Journal: Abstract and Applied Analysis
 (2014)

---

## Abstract

In terms of the poor geometric adaptability of spline element method, a geometric precision spline method, which uses the rational Bezier patches to indicate the solution domain, is proposed for two-dimensional viscous uncompressed Navier-Stokes equation. Besides fewer pending unknowns, higher accuracy, and computation efficiency, it possesses such advantages as accurate representation of isogeometric analysis for object boundary and the unity of geometry and analysis modeling. Meanwhile, the selection of B-spline basis functions and the grid definition is studied and a stable discretization format satisfying inf-sup conditions is proposed. The degree of spline functions approaching the velocity field is one order higher than that approaching pressure field, and these functions are defined on one-time refined grid. The Dirichlet boundary conditions are imposed through the Nitsche variational principle in weak form due to the lack of interpolation properties of the B-splines functions. Finally, the validity of the proposed method is verified with some examples.


---

## Body

## 1. Introduction

Comparing with traditional finite element method, spline element method (on the basis of Galerkin Principle and Spline Function Theory) involves less calculation, higher precision, and fewer pending unknowns and it is easier to construct high-order coordination unit. Thus, it has attracted much attention, and Chinese scholars have gained much achievement [1, 2]. However, it mainly focuses on structural problems [3], such as elastic beam, shell, vibration, and dynamic response and the research on nonstructural problem like fluid is far from enough. Currently, the two main methods of the application of spline function in fluid mechanics are Collocation Method and Galerkin Method.Spline Collocation Method is similar to Chebyshev Spectrum Method in its less calculation and higher efficiency. So Botella [4] applied it to calculate the incompressible Navier-Stokes. Aiming at solving false oscillation by suppressing the pressure value, Botella [5] proposed a staggered grid collocation scheme and achieved stable numerical results. Comparing to the collocation method, Galerkin Method has higher numerical precision and maturer error analysis theory. However, an element type (e.g., Taylor-Hood Element) that satisfies inf-sup stability condition [6] needs to be constructed when Galerkin Method is applied for Navier-Stokes Equations. In the field of spline element, Kumar et al. [7] had adopted weighted extended B-splines (WEB-spline) to compute Stokes. Then, they extended to Navier-Stokes equations [8] containing nonlinear convection term and constructed stable grid discrete format. The basic idea was that the degree of spline function approaching velocity field was one order higher than that approaching pressure field while only two kinds of discrete formats, namely, linear-constant and quadratic-linear, are designed. Meanwhile, WEB-spline Method is a meshless method. It avoided cockamamie grid division by replacing unstructured grid of finite element with regular net, but boundary elements require special treatment. B-Spline Element Method was adopted by Kravchenko et al. [9, 10] to analyze turbulent flow problem to decrease the calculation amounts of large eddy simulation and direct numerical method as well as to increase resolution ratio of boundary layer by embedding partitioned grid. In addition, they adopted divergence-free B-spline to expand and eliminate pressure term in governing equations to decrease numerical disturbance of calculating results.Moreover, two disadvantages of Spline Element Method have been noticed. One is that it is restricted by “low geometric versatility” and is only appropriate for the solving domain of specially simple geometric shape (e.g., rectangular or those that can be converted into rectangular). The other is that B-spline function has no interpolation property, and the function value is in the convex hull. So Dirichlet boundary condition cannot be imposed directly at the junction. Mingquan [13] solved the first problem by converting quadrilateral area into rectangular region through double linear coordinate transformation. Ronglin et al. [14] calculated boundary value of the one with arc boundary through polar mapping. However, these attempts have failed to fully address this issue. Hughes et al. [15] put forward Isogeometric Analysis Method to bridge geometric modeling and finite element analysis. It can be applied in any complex geometric area, but the primary function needs to be rational function and it is less efficient than finite element and spline element.This paper aims at solving incompressible Navier-Stokes equation and the main idea is as follows.(
1
) On the basis of Isogeometric Analysis Method, solution domain can be precisely represented by making rational Bezier patches as geometric mapping and the spline element can be ascertained with the geometry that evens the B-spline function approaching physical field. Appropriate function space can be more flexibly chosen by separating the expressions of geometric solving domain and physical field; (
2
) construct discrete format of stable spline element that meets inf-sup conditions; (
3
) impose essential boundary condition through Nitsche variational principle for B-spline function’s lack of interpolation property.


## 2. Flow of Navier-Stokes

Assume that the boundary∂
Ω of a 2D closed connected region Ω
∈
R
2 satisfies Lipschitz succession. The incompressible Navier-Stokes flow equation in dimensionless form is
(1)
-
μ
Δ
u
+
(
u
·
∇
)
u
+
∇
p
=
f
,
∇
u
=
0
,
in
Ω
,
u
=
(
u
,
v
) refers to velocity vector of fluid, p refers to the pressure, and f
=
(
f
1
,
f
2
) refers to volume force source term. μ
=
1
/
Re stands for scale-free viscosity coefficient, in which Reynolds number Re
=
ρ
U
L
/
ν is a dimensionless number of representational fluid property. Nonlinear term of convection form u
·
∇
u
=
u
j
(
∂
u
i
/
∂
x
j
) is adopted in this paper, mainly because of its simple format and numerical stability of high Reynolds number. Additional boundary condition and distribution constraint of pressure field should be added to solve the above equation:
(2)
u
=
g
,
on
⁡
∂
Ω
,
∫
Ω
p
d
x
d
y
=
0
,
in
Ω
,
where g
=
(
g
1
,
g
2
) refers to Dirichlet boundary condition of speed on boundary ∂
Ω. The second equation means that the average pressure is zero.Assume that a function spaceW
=
{
p
∈
L
2
(
Ω
)
:
∫
Ω
p
d
x
d
y
=
0
}, and a vector function space S
=
{
u
∈
H
1
(
Ω
)
×
H
1
(
Ω
)
,
u
=
g
}, V
=
{
u
∈
H
1
(
Ω
)
×
H
1
(
Ω
)
,
u
=
0
} exist; then the weak form solution of constant Navier-Stokes equation can be expressed as search function (
u
,
p
)
∈
S
×
W, and it satisfies
(3)
∫
Ω
μ
∇
w
:
∇
u
d
x
+
∫
Ω
w
(
u
·
∇
u
)
d
x
-
∫
Ω
(
∇
·
w
)
p
d
x
=
∫
Ω
w
f
d
x
-
∫
Ω
θ
∇
·
u
d
x
=
0
.
In the equation, arbitrary function (
w
,
θ
)
∈
V
×
W, Galerkin discretization assumes to project physical quantities into a finite dimensional subspace, and the speed and pressure are approximately expressed as u
h
=
∑
i
=
1
n
u
N
i
u
(
x
)
{
u
i
v
i
} and p
h
=
∑
i
=
1
n
p
N
i
p
(
x
)
p
i, respectively. N
i
u and N
i
p, respectively, stand for the primary functions of finite element space of speed and pressure field. n
u and n
p are numbers of primary functions. This study is different from traditional finite element method in adopting B-spline function as the primary function.


## 3. Solution of Navier-Stokes Equation

### 3.1. Nitsche Type Variational Weak Form

In order to simplify the derivation of variational weak form, nonlinear termu
·
∇
u can be ignored, and the equitation is reduced to Stokes flow equation: -
μ
Δ
u
+
∇
p
=
f. Its weak form is equivalent to variational extremum: search u
∈
S and then obtained as
(4)
J
(
u
)
=
min
⁡
w
∈
S
J
(
w
)
,
J
(
w
)
=
1
2
∫
Ω
μ
∇
w
:
∇
w
d
x
-
∫
Ω
w
f
d
x
.Then the following constraints should be obeyed:(
1
) incompressible condition: ∇
·
u
=
0; (
2
) pressure constraint: ∫
Ω
p
d
x
d
y
=
0. It should be noticed that the above conclusion is made on the basis of natural variational principle, so Dirichlet boundary condition u
|
∂
Ω
=
g should be met when finite element space is constructed. But it can be known from the above analysis that B-spline has no interpolation property, so the constraints are hard to be directly imposed as to traditional polynomial unit (e.g., interpolation of Lagrangian unit).This paper obtained unconstrained functional by imposing Dirichlet boundary condition with Nitsche method [16] and incompressibility and pressure condition with Lagrangian multiplier method:
(5)
J
L
(
w
)
=
J
(
w
)
-
∫
Ω
p
(
∇
·
w
)
d
x
+
γ
∫
Ω
p
d
x
+
∫
∂
Ω
λ
(
w
-
g
)
d
s
+
β
2
∫
∂
Ω
(
w
-
g
)
2
d
s
.
In the equation, p, γ, λ
=
{
λ
1
,
λ
2
} stands for Lagrangian multipliers and constant β is the penalty factor that depends on mesh size h. After taking the stationary value of above unconstrained functional and several transformations, Lagrangian multipliers λ
1
=
-
(
∂
u
/
∂
n
)
+
p
n
1 and λ
2
=
-
(
∂
v
/
∂
n
)
+
p
n
2 can be identified, which means the normal component of velocity gradient and pressure. Here n
=
(
n
1
,
n
2
) refers to normal vector outside the unit. After substituting it into variational weak form and adding nonlinear convection term, the following equation can be obtained:
(6)
∫
μ
∇
w
:
∇
u
d
x
-
∫
∂
Ω
μ
∂
w
∂
n
u
d
s
-
∫
∂
Ω
μ
w
∂
u
∂
n
d
s
+
β
∫
∂
Ω
w
u
d
s
+
∫
Ω
w
(
u
·
∇
u
)
d
x
-
∫
Ω
(
∇
·
w
)
p
d
x
+
∫
∂
Ω
w
n
p
d
s
=
∫
Ω
w
f
d
x
-
∫
∂
Ω
μ
∂
w
∂
n
g
d
s
+
β
∫
∂
Ω
w
g
d
s
-
∫
θ
∇
·
u
d
x
+
∫
∂
Ω
θ
n
u
d
s
+
(
∫
Ω
θ
d
x
)
γ
=
∫
∂
Ω
θ
n
g
d
s
,
∫
Ω
p
d
x
=
0
.
Equation (6) is different from (29) in [17] in that Lagrangian multiplier λ
i contains pressure part n
i, which makes the first two equations of above formula contain boundary pressure term ∫
∂
Ω
w
n
p
d
s, ∫
∂
Ω
θ
n
u
d
s, and ∫
∂
Ω
θ
n
g
d
s. It should be noticed that the modified weak form equation can ensure optimal order convergence of numerical solution.Assume that the speedu
=
(
u
1
,
u
2
)
=
(
u
,
v
) and pressure field p can be approximately expressed by adopting spline function:
(7)
u
=
∑
i
=
1
n
u
N
i
u
1
(
x
,
y
)
u
i
,
v
=
∑
i
=
1
n
v
N
i
u
2
(
x
,
y
)
p
i
p
=
∑
i
=
1
n
p
N
i
p
(
x
,
y
)
p
i
.
In above equation, velocity components u, v and pressure p adopt different primary functions, namely, N
i
u, N
i
v, and N
i
p. Then nonlinear simultaneous equations can be obtained through settlement after substituting them into the formula
(8)
[
A
1
+
C
1
0
-
G
1
+
R
1
0
0
A
2
+
C
2
-
G
2
+
R
2
0
-
G
1
T
+
R
1
T
-
G
2
T
+
R
2
T
0
W
0
0
W
T
0
]
(
u
^
v
^
p
^
γ
)
=
(
B
1
B
2
Q
0
)
.
In the equation, u
^
=
{
u
i
}, v
^
=
{
v
i
}, and p
^
=
{
p
i
} vectors need to be solved. Partitioned matrix: A
k
=
K
k
-
[
H
k
+
H
k
T
]
+
β
M
k, k
=
1,2. Element of matrix K
k: K
i
j
k
=
∫
Ω
μ
∇
N
i
u
k
·
∇
N
j
u
k
d
x
d
y; element of matrix H
k: H
i
j
k
=
∫
∂
Ω
μ
(
∂
N
i
u
k
/
∂
n
)
N
j
u
k
d
s; element of matrix M
k: M
i
j
k
=
∫
∂
Ω
N
i
u
k
N
j
u
k
d
s; element of nonlinearity matrix C
k: C
i
j
k
=
∫
Ω
N
i
u
k
(
u
(
∂
N
j
u
k
/
∂
x
)
+
v
(
∂
N
j
u
k
/
∂
y
)
)
d
x
d
y. Element of matrix G
1: G
i
j
1
=
∫
Ω
(
∂
N
i
u
/
∂
x
)
N
j
p
d
x
d
y, element of matrix G
2: G
i
j
2
=
∫
Ω
(
∂
N
i
v
/
∂
y
)
N
j
p
d
x
d
y; element of matrix R
k: R
i
j
k
=
∫
∂
Ω
N
i
u
k
n
k
N
j
p
d
s. Element of vector W: W
i
=
∫
Ω
N
i
p
(
x
)
d
x
d
y. Right-hand member of term: B
k
=
F
k
-
T
k
+
β
S
k; element of vector F
k: F
i
k
=
∫
Ω
N
i
u
k
f
k
d
x
d
y; element of vector T
k: T
i
k
=
∫
∂
Ω
μ
(
∂
N
i
u
k
/
∂
n
)
g
k
d
s; element of vector S
k: S
i
k
=
∫
∂
Ω
N
i
u
k
g
k
d
s; element of vector Q: Q
i
=
∫
∂
Ω
N
i
p
(
g
·
n
)
d
s.


### 3.2. Solving of Nonlinear Equation

This study adopts Newton-Raphson method to solve the nonlinear equations (9) because matrix C
k contains nonlinear term of displacement field. First, (9) is expressed as vector form L
(
a
)
≡
0, in which L
(
a
)
=
(
L
1
,
L
2
,
L
3
,
L
4
) and vector quantity a
=
(
u
^
,
v
^
,
p
^
,
γ
):
(9)
L
1
=
A
1
u
^
+
C
1
(
u
^
)
u
^
-
G
1
p
^
+
R
1
p
^
-
B
1
,
L
2
=
A
2
v
^
+
C
2
(
u
^
)
v
^
-
G
2
p
^
+
R
2
p
^
-
B
2
,
L
3
=
(
-
G
1
T
+
R
1
T
)
u
^
+
(
-
G
2
T
+
R
2
T
)
v
^
+
γ
W
,
L
4
=
W
T
p
.
Conduct Taylor expansion on vector L
(
a
) and ignore high-order term and then obtain 0
≡
L
(
a
(
n
)
+
Δ
a
(
n
)
)
≈
L
(
u
(
n
)
)
+
(
∂
L
/
∂
a
)
a
=
a
(
n
)
Δ
a
(
n
), in which Jacobian matrix is
(10)
K
T
=
∂
L
∂
a
=
[
∂
L
1
∂
u
^
∂
L
1
∂
v
^
∂
L
1
∂
p
^
∂
L
1
∂
γ
∂
L
2
∂
u
^
∂
L
2
∂
v
^
∂
L
2
∂
p
^
∂
L
2
∂
γ
∂
L
3
∂
u
^
∂
L
3
∂
v
^
∂
L
3
∂
p
^
∂
L
3
∂
γ
∂
L
4
∂
u
^
∂
L
4
∂
v
^
∂
L
4
∂
p
^
∂
L
4
∂
γ
]
.
In the equations, for matrix ∂
L
1
/
∂
u
^
=
A
1
+
C
1
+
D
1, element of matrix D
1: D
i
j
1
=
∫
Ω
N
i
u
(
∂
u
/
∂
x
)
N
j
u
d
x
d
y, matrix ∂
L
1
/
∂
v
^
=
E
1, its element: E
i
j
1
=
∫
Ω
N
i
u
(
∂
u
/
∂
y
)
N
j
v
d
x
d
y, matrix ∂
L
1
/
∂
p
^
=
-
G
1
+
R
1, vector ∂
L
1
/
∂
γ
=
0. For matrix ∂
L
2
/
∂
u
^
=
E
2, its element: E
i
j
2
=
∫
Ω
N
i
v
(
∂
v
/
∂
x
)
N
j
u
d
x
d
y, matrix ∂
L
2
/
∂
v
^
=
A
2
+
C
2
+
D
2, in which, element of matrix D
2: D
i
j
2
=
∫
Ω
N
i
v
(
∂
v
/
∂
y
)
N
j
v
d
x
d
y, ∂
L
2
/
∂
p
^
=
-
G
2
+
R
2, vector ∂
L
2
/
∂
γ
=
0. For matrix ∂
L
3
/
∂
u
^
=
-
G
1
T
+
R
1
T, matrix ∂
L
3
/
∂
v
^
=
-
G
2
T
+
R
2
T, matrix ∂
L
3
/
∂
p
^
=
0, vector ∂
L
3
/
∂
γ
=
W. Vector ∂
L
4
/
∂
u
^
=
0, vector ∂
L
4
/
∂
v
^
=
0, vector ∂
L
4
/
∂
p
^
=
W
T, scalar ∂
L
4
/
∂
γ
=
0.The equation can be solved with Newton-Raphson iteration method:(11)
K
T
Δ
a
(
n
)
=
-
L
(
a
(
n
)
)
,
a
(
n
+
1
)
=
a
(
n
)
+
ω
Δ
a
(
n
)
.
In the equations, n refers to iterative times and 0
<
ω
≤
1 refers to relaxation factor.


## 3.1. Nitsche Type Variational Weak Form

In order to simplify the derivation of variational weak form, nonlinear termu
·
∇
u can be ignored, and the equitation is reduced to Stokes flow equation: -
μ
Δ
u
+
∇
p
=
f. Its weak form is equivalent to variational extremum: search u
∈
S and then obtained as
(4)
J
(
u
)
=
min
⁡
w
∈
S
J
(
w
)
,
J
(
w
)
=
1
2
∫
Ω
μ
∇
w
:
∇
w
d
x
-
∫
Ω
w
f
d
x
.Then the following constraints should be obeyed:(
1
) incompressible condition: ∇
·
u
=
0; (
2
) pressure constraint: ∫
Ω
p
d
x
d
y
=
0. It should be noticed that the above conclusion is made on the basis of natural variational principle, so Dirichlet boundary condition u
|
∂
Ω
=
g should be met when finite element space is constructed. But it can be known from the above analysis that B-spline has no interpolation property, so the constraints are hard to be directly imposed as to traditional polynomial unit (e.g., interpolation of Lagrangian unit).This paper obtained unconstrained functional by imposing Dirichlet boundary condition with Nitsche method [16] and incompressibility and pressure condition with Lagrangian multiplier method:
(5)
J
L
(
w
)
=
J
(
w
)
-
∫
Ω
p
(
∇
·
w
)
d
x
+
γ
∫
Ω
p
d
x
+
∫
∂
Ω
λ
(
w
-
g
)
d
s
+
β
2
∫
∂
Ω
(
w
-
g
)
2
d
s
.
In the equation, p, γ, λ
=
{
λ
1
,
λ
2
} stands for Lagrangian multipliers and constant β is the penalty factor that depends on mesh size h. After taking the stationary value of above unconstrained functional and several transformations, Lagrangian multipliers λ
1
=
-
(
∂
u
/
∂
n
)
+
p
n
1 and λ
2
=
-
(
∂
v
/
∂
n
)
+
p
n
2 can be identified, which means the normal component of velocity gradient and pressure. Here n
=
(
n
1
,
n
2
) refers to normal vector outside the unit. After substituting it into variational weak form and adding nonlinear convection term, the following equation can be obtained:
(6)
∫
μ
∇
w
:
∇
u
d
x
-
∫
∂
Ω
μ
∂
w
∂
n
u
d
s
-
∫
∂
Ω
μ
w
∂
u
∂
n
d
s
+
β
∫
∂
Ω
w
u
d
s
+
∫
Ω
w
(
u
·
∇
u
)
d
x
-
∫
Ω
(
∇
·
w
)
p
d
x
+
∫
∂
Ω
w
n
p
d
s
=
∫
Ω
w
f
d
x
-
∫
∂
Ω
μ
∂
w
∂
n
g
d
s
+
β
∫
∂
Ω
w
g
d
s
-
∫
θ
∇
·
u
d
x
+
∫
∂
Ω
θ
n
u
d
s
+
(
∫
Ω
θ
d
x
)
γ
=
∫
∂
Ω
θ
n
g
d
s
,
∫
Ω
p
d
x
=
0
.
Equation (6) is different from (29) in [17] in that Lagrangian multiplier λ
i contains pressure part n
i, which makes the first two equations of above formula contain boundary pressure term ∫
∂
Ω
w
n
p
d
s, ∫
∂
Ω
θ
n
u
d
s, and ∫
∂
Ω
θ
n
g
d
s. It should be noticed that the modified weak form equation can ensure optimal order convergence of numerical solution.Assume that the speedu
=
(
u
1
,
u
2
)
=
(
u
,
v
) and pressure field p can be approximately expressed by adopting spline function:
(7)
u
=
∑
i
=
1
n
u
N
i
u
1
(
x
,
y
)
u
i
,
v
=
∑
i
=
1
n
v
N
i
u
2
(
x
,
y
)
p
i
p
=
∑
i
=
1
n
p
N
i
p
(
x
,
y
)
p
i
.
In above equation, velocity components u, v and pressure p adopt different primary functions, namely, N
i
u, N
i
v, and N
i
p. Then nonlinear simultaneous equations can be obtained through settlement after substituting them into the formula
(8)
[
A
1
+
C
1
0
-
G
1
+
R
1
0
0
A
2
+
C
2
-
G
2
+
R
2
0
-
G
1
T
+
R
1
T
-
G
2
T
+
R
2
T
0
W
0
0
W
T
0
]
(
u
^
v
^
p
^
γ
)
=
(
B
1
B
2
Q
0
)
.
In the equation, u
^
=
{
u
i
}, v
^
=
{
v
i
}, and p
^
=
{
p
i
} vectors need to be solved. Partitioned matrix: A
k
=
K
k
-
[
H
k
+
H
k
T
]
+
β
M
k, k
=
1,2. Element of matrix K
k: K
i
j
k
=
∫
Ω
μ
∇
N
i
u
k
·
∇
N
j
u
k
d
x
d
y; element of matrix H
k: H
i
j
k
=
∫
∂
Ω
μ
(
∂
N
i
u
k
/
∂
n
)
N
j
u
k
d
s; element of matrix M
k: M
i
j
k
=
∫
∂
Ω
N
i
u
k
N
j
u
k
d
s; element of nonlinearity matrix C
k: C
i
j
k
=
∫
Ω
N
i
u
k
(
u
(
∂
N
j
u
k
/
∂
x
)
+
v
(
∂
N
j
u
k
/
∂
y
)
)
d
x
d
y. Element of matrix G
1: G
i
j
1
=
∫
Ω
(
∂
N
i
u
/
∂
x
)
N
j
p
d
x
d
y, element of matrix G
2: G
i
j
2
=
∫
Ω
(
∂
N
i
v
/
∂
y
)
N
j
p
d
x
d
y; element of matrix R
k: R
i
j
k
=
∫
∂
Ω
N
i
u
k
n
k
N
j
p
d
s. Element of vector W: W
i
=
∫
Ω
N
i
p
(
x
)
d
x
d
y. Right-hand member of term: B
k
=
F
k
-
T
k
+
β
S
k; element of vector F
k: F
i
k
=
∫
Ω
N
i
u
k
f
k
d
x
d
y; element of vector T
k: T
i
k
=
∫
∂
Ω
μ
(
∂
N
i
u
k
/
∂
n
)
g
k
d
s; element of vector S
k: S
i
k
=
∫
∂
Ω
N
i
u
k
g
k
d
s; element of vector Q: Q
i
=
∫
∂
Ω
N
i
p
(
g
·
n
)
d
s.


## 3.2. Solving of Nonlinear Equation

This study adopts Newton-Raphson method to solve the nonlinear equations (9) because matrix C
k contains nonlinear term of displacement field. First, (9) is expressed as vector form L
(
a
)
≡
0, in which L
(
a
)
=
(
L
1
,
L
2
,
L
3
,
L
4
) and vector quantity a
=
(
u
^
,
v
^
,
p
^
,
γ
):
(9)
L
1
=
A
1
u
^
+
C
1
(
u
^
)
u
^
-
G
1
p
^
+
R
1
p
^
-
B
1
,
L
2
=
A
2
v
^
+
C
2
(
u
^
)
v
^
-
G
2
p
^
+
R
2
p
^
-
B
2
,
L
3
=
(
-
G
1
T
+
R
1
T
)
u
^
+
(
-
G
2
T
+
R
2
T
)
v
^
+
γ
W
,
L
4
=
W
T
p
.
Conduct Taylor expansion on vector L
(
a
) and ignore high-order term and then obtain 0
≡
L
(
a
(
n
)
+
Δ
a
(
n
)
)
≈
L
(
u
(
n
)
)
+
(
∂
L
/
∂
a
)
a
=
a
(
n
)
Δ
a
(
n
), in which Jacobian matrix is
(10)
K
T
=
∂
L
∂
a
=
[
∂
L
1
∂
u
^
∂
L
1
∂
v
^
∂
L
1
∂
p
^
∂
L
1
∂
γ
∂
L
2
∂
u
^
∂
L
2
∂
v
^
∂
L
2
∂
p
^
∂
L
2
∂
γ
∂
L
3
∂
u
^
∂
L
3
∂
v
^
∂
L
3
∂
p
^
∂
L
3
∂
γ
∂
L
4
∂
u
^
∂
L
4
∂
v
^
∂
L
4
∂
p
^
∂
L
4
∂
γ
]
.
In the equations, for matrix ∂
L
1
/
∂
u
^
=
A
1
+
C
1
+
D
1, element of matrix D
1: D
i
j
1
=
∫
Ω
N
i
u
(
∂
u
/
∂
x
)
N
j
u
d
x
d
y, matrix ∂
L
1
/
∂
v
^
=
E
1, its element: E
i
j
1
=
∫
Ω
N
i
u
(
∂
u
/
∂
y
)
N
j
v
d
x
d
y, matrix ∂
L
1
/
∂
p
^
=
-
G
1
+
R
1, vector ∂
L
1
/
∂
γ
=
0. For matrix ∂
L
2
/
∂
u
^
=
E
2, its element: E
i
j
2
=
∫
Ω
N
i
v
(
∂
v
/
∂
x
)
N
j
u
d
x
d
y, matrix ∂
L
2
/
∂
v
^
=
A
2
+
C
2
+
D
2, in which, element of matrix D
2: D
i
j
2
=
∫
Ω
N
i
v
(
∂
v
/
∂
y
)
N
j
v
d
x
d
y, ∂
L
2
/
∂
p
^
=
-
G
2
+
R
2, vector ∂
L
2
/
∂
γ
=
0. For matrix ∂
L
3
/
∂
u
^
=
-
G
1
T
+
R
1
T, matrix ∂
L
3
/
∂
v
^
=
-
G
2
T
+
R
2
T, matrix ∂
L
3
/
∂
p
^
=
0, vector ∂
L
3
/
∂
γ
=
W. Vector ∂
L
4
/
∂
u
^
=
0, vector ∂
L
4
/
∂
v
^
=
0, vector ∂
L
4
/
∂
p
^
=
W
T, scalar ∂
L
4
/
∂
γ
=
0.The equation can be solved with Newton-Raphson iteration method:(11)
K
T
Δ
a
(
n
)
=
-
L
(
a
(
n
)
)
,
a
(
n
+
1
)
=
a
(
n
)
+
ω
Δ
a
(
n
)
.
In the equations, n refers to iterative times and 0
<
ω
≤
1 refers to relaxation factor.


## 4. Stable Grid Discretization

Approximate function space that satisfies LBB condition [6] (or called inf-sup condition) should be constructed when mixed finite element method is applied to solve Navier-Stokes equation
(12)
inf
⁡
p
∈
W
sup
⁡
u
∈
S
(
-
∫
Ω
p
∇
·
u
d
Ω
|
u
|
1
∥
p
∥
0
)
≥
α
>
0
.
In the equation, α refers to the constant that is independent of grid discretization. It is theoretically difficult to prove that certain unit format meets above condition. Moreover, perfect error analysis theory on spline element method has not been established and there are only a small amount of literatures for [18]. Therefore, the stability of grid discretization is verified with a kind of numerical test called “inf-sup test,” which is similar to the patch test that proves whether the nonconforming finite element is in convergence. It is an effective tool to verify unit quality.Here the method of inf-sup test mentioned in [19] is briefly introduced. The above LBB condition can be expressed as a discrete version
(13)
α
h
=
inf
⁡
u
^
h
sup
⁡
p
^
h
(
u
^
h
T
G
p
^
h
u
^
h
T
K
u
^
h
v
h
T
Q
v
h
)
,
where element of matrix Q: Q
i
j
=
∫
Ω
N
i
p
N
j
p
d
x, matrix K
=
[
K
1
0
0
K
2
], and matrix G
=
[
G
1
G
2
]. Then generalized eigenvalue problem can be approached through a series of transformations:
(14)
P
φ
=
λ
Q
φ
.
In the equation, matrix P
=
G
T
K
-
1
G. If the eigenvalue sequence of above problem is 0
=
λ
1
=
λ
2
⋯
=
λ
k
-
1
<
λ
k
≤
λ
k
+
1
⋯
λ
n, then inf-sup constant is α
h
=
λ
k, namely, the square root of the smallest nonzero eigenvalue. Inf-sup test requires that the inf-sup constant α
h should be independence of mesh size h.The approximation capability of spline function depends on the function power and grid density, so the approximation precision could be improved through the promotion of power and grid density. Assume that parameter regionD is equally divided into n shares along arbitrary coordinate direction; then the total number of units is N
2 and the mesh size is h
=
1
/
N, and initial mesh is denoted by π
0
(
h
). As shown in Figure 1, take a unit from the grid (Figure 1(a)) and then equally divide them into four small units (Figure 1(b)). After the refinement, the total number of units is (
2
N
)
2, the mesh size is h
/
2, and the new grid is denoted by π
1
(
h
/
2
).Mesh refinement.
(a)
(b)This paper adopts such a stable discrete format, namely, the power of spline function closed to velocity fieldu is one order higher than that of spline function closed to pressure field p. What is more, velocity field u adopts one grid refinement π
1
(
h
/
2
) and pressure field p adopts original grid π
0
(
h
). As a contrast, in another unstable discrete format, velocity field and pressure field share the grids with same density and same power of primary function. For the convenience of later reference, we called the former as 4/1 format and the later as 1/1 format. What is more, such a mark u
c
a
p
d
b is adopted, in which u refers to velocity field and p refers to pressure field. The superscript a (or b) refers to power of spline function and the subscript c (or d) is 0 or 1, in which value 0 refers to origin gird and value 1 refers to refined grid.Now, a test on a group of numbers is conducted using two kinds of geometry regions shown on Figure2. Equally divide them in parameter regions and adopt n
=
1,2
,
…
,
5 for the power of spline function closed to pressure field, respectively. For each group of splines, continuously defined grids N
×
N will be adopted, in which the number of units in each direction is N
=
5,10,20,40. The numerical results are shown on Figure 3: horizontal ordinate is mesh size h
=
1
/
N and vertical coordinate is constant α
h of inf-sup. It is shown in the figure that the inf-sup constant of 4/1 format is independent of mesh size h, but the inf-sup constant of 1/1 format is gradually decreased with continuous refinement of grids.Grid division.
(a)
Square area
(b)
Circular areaInf-sup constant.
(a)
Square area
(b)
Circular area


## 5. Examples of Numerical Calculation

### 5.1. Rectangular Area

Taking 2D Navier-Stokes flow that is defined within unit rectangular areaΩ
=
[
0,1
]
×
[
0,1
] into account, assume the analysis formulas of flow function ψ
(
x
,
y
) and pressure function p
(
x
,
y
) are
(15)
ψ
(
x
,
y
)
=
sin
⁡
(
π
x
)
sin
⁡
(
π
y
)
π
,
p
(
x
,
y
)
=
cos
⁡
(
π
x
)
2
+
cos
⁡
(
π
y
)
2
2
+
2
π
cos
⁡
(
π
x
)
cos
⁡
(
π
y
)
.
So the components of velocity field are u
=
∂
ψ
/
∂
y
=
sin
(
π
x
)
cos
⁡
(
π
y
) and v
=
-
∂
ψ
/
∂
x
=
-
cos
⁡
(
π
x
)
sin
(
π
y
), and the components of volume force source term are f
1
=
0 and f
2
=
-
4
π
2
cos
⁡
(
π
x
)
sin
(
π
y
). Adopt grid division shown on Figure 3(a) and add Dirichlet velocity boundary conditions on the assumption of viscosity coefficient being μ
=
1.After a group of tests, it is found that the unit number from each direction isN
=
5,10,20,40 when the degree of spline function closed to pressure field is n
=
1,2
,
…
5, and t. Figure 4(a) shows the convergence curve of approximate velocity field u
h under L2-norm, in which horizontal ordinate refers to mesh size h and vertical coordinate refers to error of L
2-norm. And the convergence rate value is marked beside every curve. Similarly, Figure 4(b) and Figure 5(a) show the convergence curves of approximate velocity field u
h under H
1-norm and approximate solution of pressure field p
h under L2-norm. It should be noticed that if 4/1 discrete format is adopted, then solution of Navier-Stokes equation through spline function can get optimal convergence rate (optimal convergence rate refers to convergence of p
+
1 order under L
2-norm, convergence of p order under H
1-norm.). But for unstable 1/1 discrete format, false numerical oscillation will exist on pressure field p
h (especially u
0
2
p
0
2 and u
0
3
p
0
3, as is shown as Figure 5(b)).Convergence rate of velocity of Section5.1.
(a)
L
2-norm errors
(b)
H
2-norm errorsConvergence rate of pressure of Section5.1.
(a)
L
2-norm errors of 4/1 mesh format
(b)
L
2-norm errors of 1/1 mesh format


### 5.2. Cavity Flow

Lid-driven flow within unit square areaΩ
=
[
0,1
]
×
[
0,1
] is approached in this part. As is shown in Figure 6, side wall and cavity bottom is fix as (u
=
0, v
=
0), the lid is moving with constant velocity (u
=
1, v
=
0), and the fluid in the cavity flows with the driven of surface viscous force. Under moderate condition of Re number, except primary vortex (PV) in the center of square cavity, there still exists secondary vortex (SV) on the corner of cavity bottom.Figure 6
Lid-driven cavity flow.Now, the postprocessing is presented. Due to the relationship between flow function and velocity fieldu
=
(
u
,
v
)
(16)
-
(
∂
2
ψ
∂
x
2
+
∂
2
ψ
∂
y
2
)
=
∂
v
∂
x
-
∂
u
∂
y
ψ
|
∂
Ω
=
h
.
In the equation, h
(
x
)
-
h
(
x
0
)
=
∫
(
x
0
,
x
^
)
∇
ψ
·
τ
d
s, boundary tangent vector τ
=
(
-
n
2
,
n
1
). After solving such a Dirichlet poison equation, distribution of the stream-function can be obtained through postprocessing.Distribution of the stream-functions ofRe
=
100 and Re
=
1000 is, respectively, shown in Figure 7, in which power of spline function closed to pressure field is n
=
3, unit number from each direction is N
=
20. Values of contour line in Figures 8(a) and 8(b) are −0.10, −0.09, −0.07, −0.05, −0.03, −0.01, −0.001, −0.0001, 1
×
10
-
7, 1
×
10
-
6, 1
×
10
-
5, 0.001 and −0.115, −0.11, −0.10, −0.09, −0.07, −0.05, −0.03, −0.01, −0.001, 1
×
10
-
6, 0.001, 0.002, 0.005, 0.001, 0.0015, 0.0017, respectively. After listing minimum stream function value of primary vortex and vortex center position in Table 1 and comparing with the date of literature [11, 12], it is found that the results are basically identical.Table 1
Comparison on minimum stream function value of primary vortex and position of vortex center.
Case
x
y
ψ
min
⁡
Re = 100/present
0.6160
0.7360
−0.103523
Re = 100/reference [11]
0.6172
0.7344
−0.103423
Re = 1000/present
0.5320
0.5640
−0.118683
Re = 1000/reference [11]
0.5313
0.5625
−0.117929
Re = 1000/reference [12]
0.5300
0.5650
−0.118885Streamline distribution of fluid (p
=
3, N
=
20).
(a)
Re
=
100
(b)
Re
=
1000Velocity magnitude along parallel centerline.
(a)
Re
=
100
(b)
Re
=
1000At last, in order to visually compare the results, velocity magnitudes on parallel centerline and vertical centerline are given. In Figure8, horizontal ordinate refers to x-coordinate and vertical coordinate refers to velocity component v. In Figure 9, horizontal ordinate refers to velocity component u and vertical coordinate refers to y-coordinate vertical centerline. The grid of 20
×
20 is adopted in the calculation, and power of spline function closed to pressure field is n
=
1,2
,
…
,
5. The data chosen as testing benchmark are from Literature [11, 12], which coincides with the results in this paper.Velocity magnitude along vertical centerline.
(a)
Re
=
100
(b)
Re
=
1000


### 5.3. Unit Circular Area

2D Navier-Stokes flow in unit circular area (the center is on original point) is defined. Assume the analysis formulas of flow functionψ
(
x
,
y
) and pressure function p
(
x
,
y
) are
(17)
ψ
(
x
,
y
)
=
(
1
-
x
2
-
y
2
)
2
64
,
p
(
x
,
y
)
=
x
2
y
2
(
x
2
+
y
2
-
2
)
512
+
x
2
(
3
-
3
x
2
+
x
4
)
1536
+
y
2
(
3
-
3
y
2
+
y
4
)
1536
+
x
y
2
-
1
2048
.
Velocity components are u
=
(
-
1
/
16
)
y
(
1
-
x
2
-
y
2
) and v
=
(
1
/
16
)
x
(
1
-
x
2
-
y
2
); components of volume force source term are f
1
=
0 and f
2
=
x. Assume Dirichlet boundary condition and viscosity coefficient μ
=
1 are added. Parameterization secondary rational Bezier carve surface should be adopted in unit circular area, and its control vertexes and weight factors are shown on Table 2.Table 2
Control vertexes and weight factors of Section5.3.
i
d
i
,
1
d
i
,
2
d
i
,
3
w
i
,
1
w
i
,
2
w
i
,
3
1
(0, 1)
(−1, 1)
(−1, 0)
1
1
1
2
(1, 1)
(0, 0)
(−1, 1)
2
/
2
1
2
/
2
3
(1, 0)
(1, −1)
(0, −1)
1
1
1The set of spline function power and parameter grid is the same as that of Section5.1. Figure 10 shows the convergence curves of approximate velocity field u
h under L
2-norm and H
1-norm, both of which reached the optimal convergence rate. For comparison, Figure 11 shows the convergence curves of pressure numerical solution p
h under L
2-norm calculated in 4/1 grid and 1/1 grid, respectively. It is obvious that the false numerical oscillation of unstable format causes the degradation of convergence rate (Figure 11(b)).Convergence rate of velocity of Section5.3.
(a)
L
2-norm errors
(b)
H
1-norm errorsConvergence rate of pressure of Section5.3.
(a)
L
2-norm errors of 4/1 mesh format
(b)
L
2-norm errors of 1/1 mesh format


### 5.4. Circular Couette Flow

Lastly, the typical problem of Circular Couette Flow is considered. As is shown in Figure12, there are viscous incompressible fluids between two infinite-length concentric cylinders. The radiuses of outer cylinder and inner cylinder are R
1 and R
2 and they are rotated with the constant angular velocity of Ω
1 and Ω
2. Assume that the rotating velocity is slower and the fluid is in steady laminar flow phase; then there is an analytical solution of tangential velocity:
(18)
u
θ
=
A
r
+
B
r
,
A
=
Ω
2
R
2
2
-
Ω
1
R
1
2
R
2
2
-
R
1
2
,
B
=
(
Ω
1
-
Ω
2
)
R
1
2
R
2
2
R
2
2
-
R
1
2
.
In the equations, r
=
x
2
+
y
2 refers to radial coordinate. This paper assumes that fixation of outer cylinder is Ω
2
=
0 and angular velocity of inner cylinder is Ω
1
=
1. Due to its symmetry, 1/4 is taken for analysis and the settings of geometric definition, grid, and boundary condition are shown on Figure 13. Please refer to Table 3 for control vertex and weight factor within defined geometry area.Table 3
Control vertex and weight factor of Section5.3.
i
d
i
,
1
d
i
,
2
d
i
,
3
w
i
,
1
w
i
,
2
w
i
,
3
1
(1, 0)
(1, 1)
(0, 1)
1
2
/
2
1
2
(2, 0)
(2, 2)
(0, 2)
1
2
/
2
1Figure 12
Circular Couette flow.Figure 13
Settings of grid and boundary condition.Figure14 shows the distribution of tangential velocity u
θ. The power of spline function approaching pressure field is n
=
3 and the grid is 20
×
20. Figure 15 shows the distribution of tangential velocity u
θ on radial coordinate with angle of 45°. It can be noticed that they coincide with the analytical solution.Figure 14
Distribution of tangential velocityu
θ.Figure 15
Distribution of tangential velocityu
θ on radial coordinate.


## 5.1. Rectangular Area

Taking 2D Navier-Stokes flow that is defined within unit rectangular areaΩ
=
[
0,1
]
×
[
0,1
] into account, assume the analysis formulas of flow function ψ
(
x
,
y
) and pressure function p
(
x
,
y
) are
(15)
ψ
(
x
,
y
)
=
sin
⁡
(
π
x
)
sin
⁡
(
π
y
)
π
,
p
(
x
,
y
)
=
cos
⁡
(
π
x
)
2
+
cos
⁡
(
π
y
)
2
2
+
2
π
cos
⁡
(
π
x
)
cos
⁡
(
π
y
)
.
So the components of velocity field are u
=
∂
ψ
/
∂
y
=
sin
(
π
x
)
cos
⁡
(
π
y
) and v
=
-
∂
ψ
/
∂
x
=
-
cos
⁡
(
π
x
)
sin
(
π
y
), and the components of volume force source term are f
1
=
0 and f
2
=
-
4
π
2
cos
⁡
(
π
x
)
sin
(
π
y
). Adopt grid division shown on Figure 3(a) and add Dirichlet velocity boundary conditions on the assumption of viscosity coefficient being μ
=
1.After a group of tests, it is found that the unit number from each direction isN
=
5,10,20,40 when the degree of spline function closed to pressure field is n
=
1,2
,
…
5, and t. Figure 4(a) shows the convergence curve of approximate velocity field u
h under L2-norm, in which horizontal ordinate refers to mesh size h and vertical coordinate refers to error of L
2-norm. And the convergence rate value is marked beside every curve. Similarly, Figure 4(b) and Figure 5(a) show the convergence curves of approximate velocity field u
h under H
1-norm and approximate solution of pressure field p
h under L2-norm. It should be noticed that if 4/1 discrete format is adopted, then solution of Navier-Stokes equation through spline function can get optimal convergence rate (optimal convergence rate refers to convergence of p
+
1 order under L
2-norm, convergence of p order under H
1-norm.). But for unstable 1/1 discrete format, false numerical oscillation will exist on pressure field p
h (especially u
0
2
p
0
2 and u
0
3
p
0
3, as is shown as Figure 5(b)).Convergence rate of velocity of Section5.1.
(a)
L
2-norm errors
(b)
H
2-norm errorsConvergence rate of pressure of Section5.1.
(a)
L
2-norm errors of 4/1 mesh format
(b)
L
2-norm errors of 1/1 mesh format


## 5.2. Cavity Flow

Lid-driven flow within unit square areaΩ
=
[
0,1
]
×
[
0,1
] is approached in this part. As is shown in Figure 6, side wall and cavity bottom is fix as (u
=
0, v
=
0), the lid is moving with constant velocity (u
=
1, v
=
0), and the fluid in the cavity flows with the driven of surface viscous force. Under moderate condition of Re number, except primary vortex (PV) in the center of square cavity, there still exists secondary vortex (SV) on the corner of cavity bottom.Figure 6
Lid-driven cavity flow.Now, the postprocessing is presented. Due to the relationship between flow function and velocity fieldu
=
(
u
,
v
)
(16)
-
(
∂
2
ψ
∂
x
2
+
∂
2
ψ
∂
y
2
)
=
∂
v
∂
x
-
∂
u
∂
y
ψ
|
∂
Ω
=
h
.
In the equation, h
(
x
)
-
h
(
x
0
)
=
∫
(
x
0
,
x
^
)
∇
ψ
·
τ
d
s, boundary tangent vector τ
=
(
-
n
2
,
n
1
). After solving such a Dirichlet poison equation, distribution of the stream-function can be obtained through postprocessing.Distribution of the stream-functions ofRe
=
100 and Re
=
1000 is, respectively, shown in Figure 7, in which power of spline function closed to pressure field is n
=
3, unit number from each direction is N
=
20. Values of contour line in Figures 8(a) and 8(b) are −0.10, −0.09, −0.07, −0.05, −0.03, −0.01, −0.001, −0.0001, 1
×
10
-
7, 1
×
10
-
6, 1
×
10
-
5, 0.001 and −0.115, −0.11, −0.10, −0.09, −0.07, −0.05, −0.03, −0.01, −0.001, 1
×
10
-
6, 0.001, 0.002, 0.005, 0.001, 0.0015, 0.0017, respectively. After listing minimum stream function value of primary vortex and vortex center position in Table 1 and comparing with the date of literature [11, 12], it is found that the results are basically identical.Table 1
Comparison on minimum stream function value of primary vortex and position of vortex center.
Case
x
y
ψ
min
⁡
Re = 100/present
0.6160
0.7360
−0.103523
Re = 100/reference [11]
0.6172
0.7344
−0.103423
Re = 1000/present
0.5320
0.5640
−0.118683
Re = 1000/reference [11]
0.5313
0.5625
−0.117929
Re = 1000/reference [12]
0.5300
0.5650
−0.118885Streamline distribution of fluid (p
=
3, N
=
20).
(a)
Re
=
100
(b)
Re
=
1000Velocity magnitude along parallel centerline.
(a)
Re
=
100
(b)
Re
=
1000At last, in order to visually compare the results, velocity magnitudes on parallel centerline and vertical centerline are given. In Figure8, horizontal ordinate refers to x-coordinate and vertical coordinate refers to velocity component v. In Figure 9, horizontal ordinate refers to velocity component u and vertical coordinate refers to y-coordinate vertical centerline. The grid of 20
×
20 is adopted in the calculation, and power of spline function closed to pressure field is n
=
1,2
,
…
,
5. The data chosen as testing benchmark are from Literature [11, 12], which coincides with the results in this paper.Velocity magnitude along vertical centerline.
(a)
Re
=
100
(b)
Re
=
1000


## 5.3. Unit Circular Area

2D Navier-Stokes flow in unit circular area (the center is on original point) is defined. Assume the analysis formulas of flow functionψ
(
x
,
y
) and pressure function p
(
x
,
y
) are
(17)
ψ
(
x
,
y
)
=
(
1
-
x
2
-
y
2
)
2
64
,
p
(
x
,
y
)
=
x
2
y
2
(
x
2
+
y
2
-
2
)
512
+
x
2
(
3
-
3
x
2
+
x
4
)
1536
+
y
2
(
3
-
3
y
2
+
y
4
)
1536
+
x
y
2
-
1
2048
.
Velocity components are u
=
(
-
1
/
16
)
y
(
1
-
x
2
-
y
2
) and v
=
(
1
/
16
)
x
(
1
-
x
2
-
y
2
); components of volume force source term are f
1
=
0 and f
2
=
x. Assume Dirichlet boundary condition and viscosity coefficient μ
=
1 are added. Parameterization secondary rational Bezier carve surface should be adopted in unit circular area, and its control vertexes and weight factors are shown on Table 2.Table 2
Control vertexes and weight factors of Section5.3.
i
d
i
,
1
d
i
,
2
d
i
,
3
w
i
,
1
w
i
,
2
w
i
,
3
1
(0, 1)
(−1, 1)
(−1, 0)
1
1
1
2
(1, 1)
(0, 0)
(−1, 1)
2
/
2
1
2
/
2
3
(1, 0)
(1, −1)
(0, −1)
1
1
1The set of spline function power and parameter grid is the same as that of Section5.1. Figure 10 shows the convergence curves of approximate velocity field u
h under L
2-norm and H
1-norm, both of which reached the optimal convergence rate. For comparison, Figure 11 shows the convergence curves of pressure numerical solution p
h under L
2-norm calculated in 4/1 grid and 1/1 grid, respectively. It is obvious that the false numerical oscillation of unstable format causes the degradation of convergence rate (Figure 11(b)).Convergence rate of velocity of Section5.3.
(a)
L
2-norm errors
(b)
H
1-norm errorsConvergence rate of pressure of Section5.3.
(a)
L
2-norm errors of 4/1 mesh format
(b)
L
2-norm errors of 1/1 mesh format


## 5.4. Circular Couette Flow

Lastly, the typical problem of Circular Couette Flow is considered. As is shown in Figure12, there are viscous incompressible fluids between two infinite-length concentric cylinders. The radiuses of outer cylinder and inner cylinder are R
1 and R
2 and they are rotated with the constant angular velocity of Ω
1 and Ω
2. Assume that the rotating velocity is slower and the fluid is in steady laminar flow phase; then there is an analytical solution of tangential velocity:
(18)
u
θ
=
A
r
+
B
r
,
A
=
Ω
2
R
2
2
-
Ω
1
R
1
2
R
2
2
-
R
1
2
,
B
=
(
Ω
1
-
Ω
2
)
R
1
2
R
2
2
R
2
2
-
R
1
2
.
In the equations, r
=
x
2
+
y
2 refers to radial coordinate. This paper assumes that fixation of outer cylinder is Ω
2
=
0 and angular velocity of inner cylinder is Ω
1
=
1. Due to its symmetry, 1/4 is taken for analysis and the settings of geometric definition, grid, and boundary condition are shown on Figure 13. Please refer to Table 3 for control vertex and weight factor within defined geometry area.Table 3
Control vertex and weight factor of Section5.3.
i
d
i
,
1
d
i
,
2
d
i
,
3
w
i
,
1
w
i
,
2
w
i
,
3
1
(1, 0)
(1, 1)
(0, 1)
1
2
/
2
1
2
(2, 0)
(2, 2)
(0, 2)
1
2
/
2
1Figure 12
Circular Couette flow.Figure 13
Settings of grid and boundary condition.Figure14 shows the distribution of tangential velocity u
θ. The power of spline function approaching pressure field is n
=
3 and the grid is 20
×
20. Figure 15 shows the distribution of tangential velocity u
θ on radial coordinate with angle of 45°. It can be noticed that they coincide with the analytical solution.Figure 14
Distribution of tangential velocityu
θ.Figure 15
Distribution of tangential velocityu
θ on radial coordinate.


## 6. Conclusion

This paper solved the problem of incompressible Navier-Stokes flow through geometrical precise spline element method.(
1
) This method overcame the poor geometric versatility of spline element method; adoption of rational Bezier surface patch in mapping function can accurately express complex geometry areas. (
2
) It presents a stable discrete mesh format meeting in-sup condition, which expanded the spline method into fluid.This paper only discussed 2D fluid problem, but its conclusion can be directly generalized to 3D conditions. Problems to be solved are(
1
) variation of transient problems with time; (
2
) multiarea tiling problem, which shall adapt to more complicated computational area of topology; (
3
) parameterized method of complex solution domain.



---

*Source: 350682-2014-08-27.xml*