Numerical Solution of Fluid Model:

Stabilized Space-Time Finite Element in Moving Domain

In this section, we focus on the fluid model in moving domain. Solving fluid model involving moving domain is also the source of the computational challenges that we have. The spatial domain occupied by the fluid changes in time and the model to be used should be able to handle it. Here we use DSD/SST formulation that was introduced by the Team for Advanced Flow Simulation and Modeling (T?AFSM) in 1991 as a general-purpose interface-tracking technique for computation of flow problems with moving boundarieas and interfaces [13, 14, 39, 40].

3.2.1 Space-Time Finite Element Method

Deforming-Spatial Domain/Stabilized Space-Time is a method based on space-time finite element. In space-time finite element method, the dis-cretization is applied not only in space but also in time. Consequently, the spatial deformation is taken into account automatically and the dimension of the problem increases by one. In our case, a 2D problem becomes a 3D problem including the time dimension.

As will be shown in the next subsection, the formulation of DSD/SST is written over a sequence ofntspacetime slabsQn, whereQnis the slice of the space-time domain between the time levelst_nandt_n+1. In order to construct the finite element function spaces, we partition the time interval [0, T] into subintervalsTn= (tn, tn+1), wheren= 0, ..., nt−1,t0= 0 andtn+1 =T. Let Ω_nbe the spatial domain at time levelt_nwith its boundary Γ_n. We denote by Pn, the surface described by the boundary Γn asttraversesTn. Surface Pn

can be decomposed into (Pn)g and (Pn)_h where the Dirichlet and Neumann-type boundary conditions (2.3) are enforced, respectively. We define the

space-time slab Qn as the domain enclosed by the surfaces Ωn,Ωn+1 and Pn. Each Qn is decomposed into elements Q^e_n, where e= 1,2, ..,(ne)n, and the number of space-time elements may different for each space-time slab.

In each space-time slabQn, we define the finite element trial function spaces (S_u^h)nfor velocity and (S_p^h)n for pressure and the test function spaces (V_u^h)n

and (V_p^h)n= (S_p^h)n as follows

(S^h_u)n={u^h∈[H^1h(Qn)]²,u^h =g^hon (Pn)g}, (V_u^h)n={u^h∈[H^1h(Qn)]²,u^h =0on (Pn)g},

(S^h_p)_n= (V_p^h)_n={p^h∈[H^1h(Q_n)]}.

Here, H^1h(Q_n) represent the finite dimensional function space over the do-main Q_n,

H^1h(Q_n) ={u^h:Q_n→R|u^h,∂u^h

∂x ,∂u^h

∂y ,∂u^h

∂t ∈L²(Q_n),

u^h|_Q^e_n = polynomial inx, y, t} (3.34)

Figure 3.5: Space-time slabsQn−1 and Qn.

3.2.2 DSD/SST-SUPS Formulation

Deforming-Spatial Domain/Stabilized Space-Time SUPG-PSPG (DSD/SST-SUPS) is DSD/SST formulation based on Streamline-Upwind/Petrov-Galerkin (SUPG) stabilization and Pressure-Stabilizing/Petrov-Galerkin (PSPG) sta-blization. These stabilization terms assure the numerical stability of the

computations in advection-dominated flows and when using equal order interpolation functions for velocity and pressure, which simplifies the im-plementation. Moreover, we also use Least-Square on Incompressible Con-straint (LSIC) stabilization. We will review these stabilization terms in the Appendix 3.2.8.

DSD/SST-SUPS Formulation for incompressible flow can be written as follows,

Given (u^h)⁻_n, find u^h ∈ (S_u^h)n and p^h ∈(S_p^h)n such that ∀w^h ∈ (V_u^h)n and

∀q^h∈(V_p^h)n: Z

Qn

w^h·ρ(∂u^h

∂t +u^h· ∇u^h−f^h) dQ+ Z

Qn

ε(w^h) :σ(p^h,u^h) dQ

− Z

(Pn)h

w^hh^hdp+ Z

Qn

q^h∇ ·u^hdQ+ Z

Ωn

(w^h)⁺_nρ((u^h)⁺_n −(u^h)⁻_n) dΩ +

(n_el)n

X

e=1

Z

Q^e_n

1 ρ

τ₁ρ(∂w^h

∂t +u^h· ∇w^h) +τ₂∇q^h

·

ρ(∂u^h

∂t +u^h· ∇u^h−f^h)− ∇ ·σ(p^h,u^h)

dQ +

(nel)n

X

e=1

Z

Q^e_n

τ3∇ ·w^hρ∇ ·u^hdQ= 0, (3.35)

where

(u^h)^±_n = lim

→0u(t_n±).

The first four integrals are the Galerkin formulation of the problem.

Here,ε(w^h) :σ(p^h,u^h) represents component-wise scalar product between strain-rate tensor ε(w^h) and stress tensorσ(p^h,u^h). The fifth integral en-forces, weakly, the temporal continuity of the velocity field since the basis functions are discontinous from one space-time slab to another. The remain-ing terms are the stabilization terms. Note that this stabilization leads to a consistent formulation, in the sense that an exact solution still satisfies the stabilized formulation. In this formula, τ₁, τ₂, τ₃ are the stabilization parameters,τSU P G, τP SP G, τLSIC, respectively.

After space-time discretization, we obtain a nonlinear system of

equa-tions

K(U)U=F (3.36)





K₁ K₂ K₃ K4 K5 K6

K7 K8 K9









¯ u

¯ v

¯ p



=



 F₁ F2

F3



.

Here, ¯u,¯v,p¯are the approximate solutions and matrixKconsists of 9 block matrices. Each block matrix K_i, i = 1, ...,9 represents the corresponding terms in (3.35). The size of each block matrix is (2np)×(2np) where npis the number of nodes at each time level in the space-time slab Qn. We will write each term of DSD/SST formulation and use the index ito emphasize that it is related to the position in block Ki. We calculate the integrals for elementQ^e_n of space-time slab and then assemble for all elements to get the global system matrix.

MatrixK can be written in more detail as K=





K₁ K₂ K₃ K4 K5 K6

K7 K8 K9



, with

K₁ =T D₁+B₁+C₁(U) +M₁+S1₁(U) +S5₁ K₂ =B₂+S1₂(U) +S5₂

K₃ =G^T₃ +S2₃

K4 =B4+S14(U) +S54

K5 =T D5+B5+C5(U) +M5+S15(U) +S55

K6 =G^T₆ +S26

K₇ =G₇+S4₇(U) K₈ =G₈+S4₈(U) K₉ =S3₉

The individual terms above have the following meaning:

• Time-dependent term (T D): R

Q^e_nρw^{h ∂u}_∂t^hdQ

• Component-wise scalar product:R

Q^e_nε(w^h) :σdQ^e_n

This term consist of Gradient matrix for pressure (G^T) and Diffusion term (B).

• Convection term (C(U)):R

Q^e_nρw^h(u^h.∇u^h) dQ

• Jump term (M):R

Ω^l_nρ(w^h)⁺_n(u^h)⁺_ndΩ where Ω^l_n is spatial domain of l-th triangle element at time levelt_n.

• Gradient term (G);R

Q^e_nq^h.∇u^hdQ

• Stabilization terms (S(U)) :

(nel)n

X

e=1

Z

Q^e_n

1 ρ

τ₁ ρ(∂w^h

∂t +u^h· ∇w^h) +τ₂∇q^h

·

ρ(∂u^h

∂t +u^h· ∇u^h−f^h)− ∇ ·σ(p^h,u^h)

dQ +

(nel)n

X

e=1

Z

Q^e_n

τ3∇ ·w^hρ∇ ·u^hdQ

Stabilization terms consist of five parts: S1(U), S2, S3, S4(U), S5.as in (3.48) - (3.52). We discuss more about the stabilization parameters τ₁, τ₂, τ₃ in the stabilization section.

The right-hand side is written as

F =



 F₁ F2

F3



=





F_f1 +F_s11+F_h1+F_M1 Ff2 +Fs12+Fh2+FM2

Fs43



.

The vectors F1, F2, F3 have size (2np). The detailed of each term will be explained in the next subsections.

3.2.3 Solution of Discrete Problem

In order to solve the nonlinear system (3.36), we use Newton-Raphson method. We compute a correction ∆U of a current solution U^l at each iterationl, which yields a linear system

J(U^l) ∆U^l =F−K(U^l)U^l, (3.37) where J is the Jacobian matrix. We solve (3.37) by using GMRES with diagonal preconditioner.

As in matrixK, matrixJalso consists of nine block matrices.

J=





J1 J2 J3

J₄ J₅ J₆ J7 J8 J9



 with

J₁ =T D₁+B₁+C₁(U^∗) +CC₁(U^∗) +M₁+S1₁(U^∗) +S5₁+SS1₁(U^∗), J2 =B2+S12(U^∗) +S52+SS12(U^∗) +CC2(U^∗),

J3 =G^T₃ +S23,

J₄ =B₄+S1₄(U^∗) +S5₄+SS1₄(U^∗) +CC₄(U^∗),

J₅ =T D₅+B₅+C₅(U^∗) +CC₅(U) +M₅+S1₅(U^∗) +S5₅+SS1₅(U^∗), J6 =G^T₆ +S26,

J₇ =G₇+S4₇(U^∗), J₈ =G₈+S4₈(U^∗), J9 =S39,

where U^∗ = U^l, the current U at each iteration, CC(U^∗) and SS1(U^∗) be given by

• CC(U^∗) : ρR

Q^k_nρw^h(∆u^h· ∇u^∗) dQ,

• SS1(U^∗) :τ₁ρR

Q^k_n(u^∗· ∇w^h)(∆u· ∇u^∗), and similarly for other terms as in matrix K.

The right-hand sideFis slightly different with the previousF, F=



 F1

F₂ F3



=





F_f1 +Fss11 +F_h1+FM1

F_f2 +F_ss12 +F_h2+F_M2 Fs43



, (3.38)

whereFss1 is the right-hand side that comes from the derivative of stabiliza-tion S1, thus depend on U^∗ and since U^∗ is known, F_ss1 be a component of the right-hand side. We compute Fss1 for each Newton’s iteration; while other component of F are computed once in the beginning of Newton’s it-eration since they are independent of U^∗. The detailed computation of all terms will be explained in the next subsections.

3.2.4 Linear Finite Element in Global and Local Description The construction of finite element basis functions begins at the level of individual element. Firstly, we construct the element-level interpolation

functions and then put together all elements into a finite element mesh. In this way, we define the global basis functions.

In physical domain, elements in the mesh may have different size and shape. However, every element can be considered as an image of parent ele-ment, a simple geometrical shape which is defined in the parametric domain.

Generally, there are two points of view when dealing with linear finite elements: global point of view and element or local point of view [30]. Let us look at it in our case. We will work in two space dimension and use triangular elements. In space-time finite element formulation, it becomes a 3D problem. In this case, every space-time element has six nodes, i.e., three triangular nodes for each time level. In global point of view, the basis functions are considered to be defined everywhere on the domain of the boundary-value problem. Here, we have the quantities,

• Physical domain : Ω×[t_n, t_n+1].

• Coordinate of nodes of space-time slab Q^e_n :

{xⁿ_A−1,xⁿ_A,xⁿ_A+1,xⁿ⁻¹_A+1,xⁿ⁺¹_A ,xⁿ⁺¹_A+1}.

• Degrees of freedom: {uⁿ_A−1,uⁿ_A,uⁿ_A+1,uⁿ⁺¹_A−1,uⁿ⁺¹_A ,uⁿ⁺¹_A+1}.

• Shape functions: {N¯_A−1ⁿ ,N¯_Aⁿ,N¯_A+1ⁿ ,N¯_A−1ⁿ⁺¹,N¯_Aⁿ⁺¹,N¯_A+1ⁿ⁺¹}.

• Approximate solution:

u^h(x) =uⁿ_A−1N¯_A−1ⁿ +uⁿ_AN¯_Aⁿ+uⁿ_A+1N¯_A+1ⁿ +uⁿ⁺¹_A−1N¯_A−1ⁿ⁺¹+uⁿ⁺¹_A N¯_Aⁿ⁺¹ +uⁿ⁺¹_A+1N¯_A+1ⁿ⁺¹

These quantities are in terms of global parameters, i.e., global coordi-nates, global shape functions and global node ordering. The global point of view is useful in establishing the mathematical properties of the finite element method.

In local point of view, the above quantities are in terms of local param-eters, i.e., local coordinates, local shape functions and local node ordering, as follows,

• Parametrical domain : Ωb×[θ1, θ2].

• Coordinate of nodes of space-time slab Qb^e_n : {r¹₁,r¹₂,r¹₃,r²₁,r²₂,r²₃}.

• Degrees of freedom: {u¹₁,u¹₂,u¹₃,u²₁,u²₂,u²₃}.

• Shape functions: {N₁¹, N₂¹, N₃¹, N₁², N₂², N₃²}.

• Approximate solution:

u^h(r) =u¹₁N₁¹+u¹₂N₂¹+u¹₃N₃¹+u²₁N₁²+u²₂N₂²+u²₃N₃². In our implementation, we use θ₁ = −1, θ₂ = 1,r = (r, s)^T, where r, s ∈ [0,1], r+s≤1. For the specific forms of shape functions, see the following subsection. Note that in the local point of view, the nodal numbering is represented by numbering beginning with 1 to show that this is in local description. The local point of view is useful in the computer implementation of finite element mehod. We will use this local point of view in computation of the terms of DSD/SST that has been mentioned before.

The domains of the global and local description are related by the trans-formation

ξ : Ω×[tn, tn+1]→Ωb×[θ1, θ2], (3.39) such thatξ(xⁿ_A−1) =r¹₁, ξ(xⁿ_A) =r¹₂, ξ(xⁿ_A+1) =r¹₃, ξ(xⁿ⁺¹_A−1) =r¹₁, ξ(xⁿ⁺¹_A ) = r¹₂ and ξ(xⁿ⁺¹_A+1) =r²₂, where x = (x, y)^T,r= (r, s)^T. Similarly, we also can define the inverse of ξ,

ξ⁻¹ :Ωb×[θ₁, θ₂]→Ω×[t_n, t_n+1], 3.2.5 Shape Functions of Space-Time Element

In order to get the component of systemKandJin (3.36),(3.37) and right-hand side F in (3.38), we need to understand the shape functions of space-time element and their derivatives. Here, we assume some special features of space-time slab Qn, i.e., Qn has uniform thickness in time direction and all the nodes of the space-time slab are either on its upper or lower surface (Special DSD/SST) [29]. Uniform thickness means that we take the same time-step ∆t in the program. We also assume that the mesh on the upper surface of the slab is obtained by a deformation of the mesh on the lower surface which is governed by elasticity equation. The Special DSD/SST formulation offers efficiency in computational cost, i.e., it simplifies both the computation of shape function derivatives at Gaussian quadrature points and the formation of the element-level vectors and matrices. Since for each iteration at every time step, these computations need to be done for each element of the space-time mesh, these improvements will give a significant impact on the overall computational performance in a simulation.

The Special DSD/SST formulation uses shape function having a form of the tensor product of its spatial and temporal shape functions in parametric domain,

N_a^α(r, θ) =Na(r)T^α(θ), (3.40) where a= 1,2, .., n_el, α = 1,2, with n_el is the number of nodes in spatial domain. We use triangular elements, i.e., nel = 3. In this case, the shape functionsN_a(r) read

N1(r) =r, N2(r) =s, N3(r) = 1−r−s, and their derivatives,

N1,r = 1, N2,r = 0, N3,r =−1, N_1,s = 0, N_2,s= 1, N_3,s =−1,

with r, s ∈ [0,1] and r +s <= 1. The temporal shape functions can be defined by

T¹(θ) = 1

2(1−θ), T²(θ) = 1

2(1 +θ), and their derivatives,

T_θ¹ =−1

2, T_θ² = 1 2, withθ∈[−1,1].

As the mapping (3.39), the spatial coordinate of each element in para-metrical domain can be written as

x(r, θ) =

2

X

α=1 3

X

a=1

N_a^α(r, θ)x^α_a

=

2

X

α=1 3

X

a=1

N_a(r)T^α(θ)x^α_a

=

3

X

a=1

N_a(r)(T¹(θ)x¹_a+T²(θ)x²_a)

=

3

X

a=1

N_a(r)x_a(θ),

where x^α_a = (x^α_a, y^α_a) with a= 1,2,3; α= 1,2, are the coordinate of space-time element nodes in physical domain. We list formula for each component of x(r, θ) and their derivatives:

x(r, θ) =

3

X

a=1

N_a(r)x_a

=N₁(r)x₁(θ) +N₂(r)x₂(θ) +N₃(r)x₃(θ)

=N1(r)(T¹(θ)x¹₁+T²(θ)x²₁) +N2(r)(T¹(θ)x¹₂+T²(θ)x²₂) +N3(r)(T¹(θ)x¹₃ +T²(θ)x²₃)

= r 2

(1−θ)x¹₁+ (1 +θ)x²₁ + s

2

(1−θ)x¹₂+ (1 +θ)x²₂ +1−r−s

2

(1−θ)x¹₃+ (1 +θ)x²₃ . xr(θ) = ¹₂

(1−θ)(x¹₁−x¹₃) + (1 +θ)(x²₁−x²₃) . x_s(θ) = ¹₂

(1−θ)(x¹₂−x¹₃) + (1 +θ)(x²₂−x²₃) . x_θ(r) = r(x²₁−x¹₁) +s(x²₂−x¹₂) + (1−r−s)(x²₃−x¹₃)

2 .

y(r, θ) =

3

X

a=1

N_a(r)y_a

=N₁(r)y₁(θ) +N₂(r)y₂(θ) +N₃(r)y₃(θ)

=N1(r)(T¹(θ)y¹₁+T²(θ)y₁²) +N2(r)(T¹(θ)y₂¹+T²(θ)y²₂) +N3(r)(T¹(θ)y₃¹ +T²(θ)y₃²)

= r 2

(1−θ)y₁¹+ (1 +θ)y²₁ +s

2

(1−θ)y¹₂+ (1 +θ)y₂² +1−r−s

2

(1−θ)y₃¹+ (1 +θ)y²₃ . yr(θ) = ¹₂

(1−θ)(y¹₁−y₃¹) + (1 +θ)(y₁²−y²₃) . y_s(θ) = ¹₂

(1−θ)(y¹₂−y₃¹) + (1 +θ)(y²₂−y₃²) . y_θ(r) = r(y²₁−y₁¹) +s(y²₂−y₂¹) + (1−r−s)(y₃²−y¹₃)

2 .

The temporal coordinate can be written as, t(θ) =

2

X

α=1 3

X

a=1

N_a^α(r, θ)t^α_a.

Note that in special DSD/SST formulation, t^α_a = t^α, a = 1,2,3 since the nodes of a slab are either on its upper (α= 2) or lower surface (α= 1), we have

t(θ) =

2

X

α=1 3

X

a=1

N_a^α(r, θ)t^α=

2

X

α=1 3

X

a=1

N_a(r)T^α(θ)t^α =

3

X

a=1

N_a(r)

2

X

α=1

T^α(θ)t^α

= 1(T¹(θ)t¹+T²(θ)t²) = (1−θ)t¹+ (1 +θ)t²

2 = t¹+t²+θ(t²−t¹)

2 ,

and for derivatives,

t_r= 0, ts= 0, tθ= t²−t¹

2 = ∆t 2 .

The next step is to get the derivative of shape-function of space-time element ¯N_a^α(x, t) (in physical domain) with respect to x, y and t since the components of system matrix contain these derivative as well. By using chain rule,

N¯_a^α(x, t)

∂x = N_a^α(x(r, θ), t(θ))

∂r

∂r

∂x +N_a^α(x(r, θ), t(θ))

∂s

∂s

∂x N¯_a^α(x, t)

∂y = N_a^α(x(r, θ), t(θ))

∂r

∂r

∂y +N_a^α(x(r, θ), t(θ))

∂s

∂s

∂y N¯_a^α(x, t)

∂t = N_a^α(x(r, θ), t(θ))

∂θ

∂θ

∂t,

which can be rewritten in the following matrix form,





N¯_a,x^α (x, t) N¯_a,y^α (x, t) N¯_a,t^α (x, t)



=P





N_a,r^α (r, θ) N_a,s^α (r, θ) N_a,θ^α (r, θ)



=P





N_a,rT^α(θ) N_a,sT^α(θ) Na(r)T_θ^α



 (3.41)

where

P =





rx sx θx

r_y s_y θ_y r_t s_t θ_t





Note that we have already expressions for Na,r, Na,s and Tθ. However, we do not have explicit expressions for r = r(x, y, t), s = s(x, y, t), and θ = θ(x, y, t), so matrixP cannot be computed directly. It can be obtained from P =Q⁻¹, where

Q=





x_r y_r t_r x_s y_s t_s x_θ y_θ t_θ



=





xr(θ) yr(θ) 0 x_s(θ) y_s(θ) 0 x_θ(r) y_θ(r) ^∆t₂



.

We define

J^st(θ) = detQ= ∆t 2 det

xr(θ) yr(θ) xs(θ) ys(θ)

= ∆t 2 Υ(θ),

the Jacobian of the transformation from physical domain to parametrical domain in space-time element. Here, Υ(θ) can be written as,

Υ(θ) = det

x_r(θ) y_r(θ) xs(θ) ys(θ)

. (3.42)

Computing the inverse of matrix Q, we get

Q⁻¹=P = 1 J^st(θ)













∆t

2 ys −^∆t₂ yr 0

−^∆t₂ xs ∆t

2 xr 0

−x_θys+yθxs xθyr−yθxr Υ(θ)













=























ys(θ)

Υ(θ) −yr(θ) Υ(θ) 0

−xs(θ) Υ(θ)

xr(θ)

Υ(θ) 0

V_r(r, θ) V_s(r, θ) 2

∆t





















 ,

where

V_r(r, θ) = −x_θ(r)y_s(θ) +y_θ(r)x_s(θ)

J^st(θ) , V_s(r, θ) = x_θ(r)y_r(θ)−y_θ(r)x_r(θ) J^st(θ) .

Then (3.41) can be rewritten as





N¯_a,x^α (x, t) N¯_a,y^α (x, t) N¯_a,t^α(x, t)



=























ys(θ)

Υ(θ) −yr(θ) Υ(θ) 0

−x_s(θ) Υ(θ)

x_r(θ)

Υ(θ) 0

Vr(r, θ) Vs(r, θ) 2

∆t



























N_a,rT^α(θ) Na,sT^α(θ) N_a(r)T_θ^α





=















ys(θ)Na,r−yr(θ)Na,s

Υ(θ) T^α(θ)

−x_s(θ)N_a,r+x_r(θ)N_a,s

Υ(θ) T^α(θ)

(V_r(r, θ)N_a,r+V_s(r, θ)N_a,s)T_θ^α+(−1)^α

∆t N_a(r)













 .

(3.43) 3.2.6 Computation of Component System Matrix and

Right-hand Side of The Linearized System

Using (3.43), we can compute all components of the finite element systems matrices K and J in (3.36),(3.37) and right-hand side F in (3.38). The components of system matrix and the right-hand side are decomposed from spatial and temporal shape functions due to (3.40 ). We compute the integral over temporal domain using Gaussian quadrature with Gaussian quadrature points ˜θ_i and weights W_i (Appendix 3.2.10); while the integral over spatial domain will be calculated analytically for efficiency of computation reason [29]. It is possible to perform integration over spatial domain analytically, since we use linear triangular as spatial element. Here, we use the formula in [30]

Z

4

N₁^αN₂^βN₃^γ d4= α!β!γ!

(α+β+γ+ 2)!2A whereA is the area of triangular4that can be obtained from

2A= det





1 x1 y1

1 x₂ y₂ 1 x₃ y₃



,

with (xi, yi) is the coordinate of nodeiin triangular4. In the computation of component system matrix, we change the domain of integration from physical domain into parametrical domain. In this case, A represent the area of triangle in parametrical domain with A= 1

2.

We will calculate the linear terms first and then nonlinear terms (convec-tion term, stabiliza(convec-tion stab1 term, stabiliza(convec-tion stab4 term) in increment form. We denote each term with its abbreviated name as we mentioned in the earlier part of this section, e.g., T D for time-dependent term. In T D^e_i(a, b), index irepresents the position of block matrix of this term, su-perindex e is the number of space-time slab elementQ^e_n. This symbol also emphasizes that this computation is done for each element separately. aand bare the local numbers of spatial nodes inQ^e_n,a= 1,2,3; b= 1,2,3,α and β represent time level at Q^e_n, 1 for time level tn and 2 fortn+1,Qb^e_n isQ^e_n in parametrical domain, Qb⁴n represents the triangular (spatial element) inQb^e_n. 1) Time-dependent term (T D),

Z

Q^e_n

ρw^h∂u^h

∂t dQ (3.44)

T D₁^e(a, b)

=ρ Z

Q^e_n

N¯_a^α(x, t) ¯N_b,t^β(x, t) dxdt

=ρ Z

Qb^e_n

N¯_a^α(x(r, θ), t(θ)) ¯N_b,t^β (x(r, θ), t(θ))J^st(θ) drdθ

=ρ Z

Qb⁴n

2

X

i=1

N_a^α(r,θ˜i)N_b,t^β(r,θ˜i)J^st(˜θi)Wi dr

=ρ

2

X

i=1

J^st(˜θi)Wi

Z

Qb⁴n

Na(r)T^α(˜θi)N_b,t^β (r,θ˜i) dr

=ρ

2

X

i=1

J^st(˜θi)Wi

T^α(˜θi)T^β(˜θi) Z

Qb⁴n

Na(r)(Vr(r,θ˜i)N_b,r+Vs(r,θ˜i)N_b,s)dr

+ρ

2

X

i=1

J^st(˜θi)Wi

T^α(˜θi)(−1)^β

∆t Z

Qb⁴n

Na(r)N_b(r)dr

We obtain time-dependent term,

T D^e₁(a, b)















=ρ

2

X

i=1

T^α(˜θi)2A 24 h

A1+A2+2J^st(˜θi)⁽⁻¹⁾_∆t^β

i

, if a=b

=ρ

2

X

i=1

T^α(˜θi)2A 24 h

A1+A2+J^st(˜θi)⁽⁻¹⁾_∆t^β

i

, otherwise

T D₁^e(a, b)















=ρ

2

X

i=1

T^α(˜θi)2A

24 [A1+A2+Υ(˜θi)(−1)^β], if a=b

=ρ

2

X

i=1

T^α(˜θi)2A 24

h

A1+A2+^{Υ( ˜}₂^θi)(−1)β i

, otherwise

where

A1=T^β(˜θi) (^PNc(r)(−^∆xc

2 −^∆xa

2 )ys(˜θi)+P

Nc(r)(^∆yc₂ +^∆ya₂ )xs(˜θi))Nb,r

A₂=T^β(˜θ_i) (^PNc(r)(^∆xc₂ +^∆xa₂ )yr(˜θi)−P

Nc(r)(−^∆yc

2 −^∆ya

2 )xr(˜θi))N_b,s

T D^e₅(a, b) =T D^e₁(a, b)

2) Component wise scalar product

Z

Q^e_n

ε(w^h) :σdQ^e_n (3.45)

This term consist of

• Gradient matrix for pressure (G^T) (G^T₃)^e(a, b) =−

Z

Q^e_n

N¯_a,x^α (t) ¯N_b^β(x, t) dxdt

=− Z

Qb^e_n

N¯_a,x^α (t(θ)) ¯N_b^β(x(r, θ), t(θ))J^st(θ) drdθ

=− Z

Qb⁴n

2

X

i=1

N_a,x^α (˜θ_i)N_b^β(r,θ˜_i)J^st(˜θ_i)W_i dr

=−

2

X

i=1

J^st(˜θ_i)W_i Z

Qb⁴n

N_a,x(˜θ_i)T^α(˜θ_i)N_b(r)T^β(˜θ_i) dr

=−

2

X

i=1

J^st(˜θi)WiT^α(˜θi)T^β(˜θi)Na,x(˜θi) Z

Qb⁴n

Nb(r)

=−

2

X

i=1

∆t

2 Υ(˜θi)WiT^α(˜θi)T^β(˜θi)Na,x(˜θi)2A 6

=−∆t 12

2

X

i=1

Υ(˜θi)T^α(˜θi)T^β(˜θi)Na,x(˜θi)

and similarly with (G^T₅)^e(a, b) by changingNa,x withNa,y

(G^T₆)^e(a, b) =− Z

Q^e_n

N¯_a,y^α (t) ¯N_b^β(x, t) dxdt

=−∆t 12

2

X

i=1

Υ(˜θ_i)T^α(˜θ_i)T^β(˜θ_i)N_a,y(˜θ_i)

3) Diffusion term (B):

B₁^e(a, b) =µ Z

Q^e_n

2 ¯N_a,x^α (t) ¯N_b,x^β (t) + ¯N_a,y^α (t) ¯N_b,y^β (t) dxdt

=µ Z

Qb^e_n

2 ¯N_a,x^α (t(θ)) ¯N_b,x^β (t(θ))J^st(θ) drdθ +µ

Z

Qb^e_n

N¯_a,y^α (t(θ)) ¯N_b,y^β (t(θ))J^st(θ) drdθ

B₁^e(a, b) =µ Z

Qb⁴n

2

X

i=1

2N_a,x^α (˜θi)N_b,x^β (˜θi)J^st(˜θi)Wi dr

+µ Z

Qb^e_n 2

X

i=1

N_a,y^α (˜θi)N_b,y^β (˜θi)J^st(˜θi)Widr

=µ

2

X

i=1

J^st(˜θ_i)W_i Z

Qb⁴n

2N_a,x(˜θ_i)T^α(˜θ_i)N_b,x(˜θ_i)T^β(˜θ_i) dr

+µ

2

X

i=1

J^st(˜θ_i)W_i Z

Qb⁴n

N_a,y(˜θ_i)T^α(˜θ_i)N_b,y(˜θ_i)T^β(˜θ_i) dr

=µ

2

X

i=1

J^st(˜θ_i)W_i2N_a,x(˜θ_i)T^α(˜θ_i)N_b,x(˜θ_i)T^β(˜θ_i) Z

Qb⁴n

dr

+µ

2

X

i=1

J^st(˜θ_i)W_iN_a,y(˜θ_i)T^α(˜θ_i)N_b,y(˜θ_i)T^β(˜θ_i) Z

Qb⁴n

dr

=µ

2

X

i=1

∆t

2 Υ(˜θ_i)W_i2N_a,x(˜θ_i)T^α(˜θ_i)N_b,x(˜θ_i)T^β(˜θ_i)2A 2 +µ

2

X

i=1

∆t

2 Υ(˜θ_i)W_iN_a,y(˜θ_i)T^α(˜θ_i)N_b,y(˜θ_i)T^β(˜θ_i)2A 2

= µ∆t 2

2

X

i=1

Υ(˜θ_i)N_a,x(˜θ_i)T^α(˜θ_i)N_b,x(˜θ_i)T^β(˜θ_i)

+µ∆t 4

2

X

i=1

Υ(˜θ_i)N_a,y(˜θ_i)T^α(˜θ_i)N_b,y(˜θ_i)T^β(˜θ_i).

Similarly, B₂^e(a, b) =µ

Z

Q^e_n

N¯_a,y^α (t) ¯N_b,x^β (t) dxdt

=µ Z

Qb^e_n

N¯_a,y^α (t(θ)) ¯N_b,x^β (t(θ))J^st(θ) drdθ

=µ Z

Qb⁴n

2

X

i=1

N_a,y^α (˜θi)N_b,x^β (˜θi)J^st(˜θi)Wi dr

=µ

2

X

i=1

J^st(˜θi)Wi

Z

Qb⁴n

Na,y(˜θi)T^α(˜θi)N_b,x(˜θi)T^β(˜θi) dr

B^e₂(a, b) =µ

2

X

i=1

J^st(˜θi)WiNa,y(˜θi)T^α(˜θi)N_b,x(˜θi)T^β(˜θi) Z

Qb⁴n

dr

=µ

2

X

i=1

J^st(˜θi)WiNa,y(˜θi)T^α(˜θi)N_b,x(˜θi)T^β(˜θi)2A 2

=µ

2

X

i=1

∆t

2 Υ(˜θ_i)N_a,y(˜θ_i)T^α(˜θ_i)N_b,x(˜θ_i)T^β(˜θ_i)2A 2

= µ∆t 4

2

X

i=1

Υ(˜θ_i)N_a,y(˜θ_i)T^α(˜θ_i)N_b,x(˜θ_i)T^β(˜θ_i).

B₄^e(a, b) =µ Z

Q^e_n

N¯_a,x^α (t) ¯N_b,y^β (t) dxdt

= µ∆t 4

2

X

i=1

Υ(˜θi)Na,x(˜θi)T^α(˜θi)N_b,y(˜θi)T^β(˜θi).

B₅^e(a, b) =µ Z

Q^e_n

N¯_a,x^α (t) ¯N_b,x^β (t) + 2 ¯N_a,y^α (t) ¯N_b,y^β (t) dxdt

= µ∆t 4

2

X

i=1

Υ(˜θi)Na,x(˜θi)T^α(˜θi)N_b,x(˜θi)T^β(˜θi)

+µ∆t 2

2

X

i=1

Υ(˜θi)Na,y(˜θi)T^α(˜θi)N_b,y(˜θi)T^β(˜θi).

4) Jump term (M),

Z

Ω^l_n

ρ(w^h)⁺_n(u^h)⁺_ndΩ (3.46)

M₁^e(a, b) =ρ Z

Q^e_n

N_a(x)N_b(x) dx=ρ Z

Q^e_n

N_a(r)N_b(r) ˜Υ dr

=

(ρΥ˜^2A₁₂ = ₁₂^ρΥ,˜ ifa=b ρΥ˜^2A₂₄ = ₂₄^ρΥ,˜ otherwise M₅^e(a, b) =M₁^e(a, b)

5) Gradient term (G),

Z

Q^e_n

q^h.∇u^hdQ (3.47)

The integration is analogous to that of (G^T) term

G^e₇(a, b) = Z

Q^e_n

N¯_a^α(x, t) ¯N_b,x^β (t) dxdt= ∆t 12

2

X

i=1

Υ(˜θ_i)T^α(˜θ_i)T^β(˜θ_i)N_b,x(˜θ_i)

G^e₈(a, b) = Z

Q^e_n

N¯_a^α(x, t) ¯N_b,y^β (t) dxdt= ∆t 12

2

X

i=1

Υ(˜θ_i)T^α(˜θ_i)T^β(˜θ_i)N_b,y(˜θ_i)

6) Stabilization terms (S(U)) :

(nel)n

X

e=1

Z

Q^e_n

1 ρ

τ1 ρ(∂w^h

∂t +u^h· ∇w^h) +τ2∇q^h

·

ρ(∂u^h

∂t +u^h· ∇u^h−f^h)− ∇ ·σ(p^h,u^h)

dQ+

(nel)n

X

e=1

Z

Q^e_n

τ₃∇·w^hρ∇·u^hdQ In order to explain the calculation of stabilization terms, we divide the above expression into five parts.

a. Stab1 : Z

Q^e_n

τ₁1 ρ

ρ(∂w^h

∂t +u^h· ∇w^h)ρ(∂u^h

∂t +u^h· ∇u^h)

dQ (3.48)

Stab1 is nonlinear term, it will be explained later.

b. Stab2

Z

Q^k_n

ρ(∂w^h

∂t +u^h· ∇w^h)τ1

ρ∇p^hdQ (3.49)

Inx-component : τ1

R

Q^e_nρ(^∂w_∂t^h+u^∂w_∂x^h +v^∂w_∂y^h)^∂p_∂x^hdQ S2^e₃(a, b)

=τ₁ Z

Q^e_n

N¯_a,t^α(x, t) + ¯N_a,x^α (x, t)u+ ¯N_a,y^α (x, t)vN¯_b,x^β (t) dxdt

=τ₁ Z

Qb^e_n

N_a,t^α(r, θ) +N_a,x^α (θ)(^P2_γ=1^P3_c=1u^∗γc Nc(r)T^γ(θ)

N_b,x^β (θ)J^st(θ) drdθ +τ₁

Z

Qb^e_n

N_a,y^α (θ)(^P2_γ=1^P3_c=1vc^∗γNc(r)T^γ(θ)

N_b,x^β (θ)J^st(θ) drdθ

S2^e₃(a, b)

=τ₁ Z

Qb⁴n

2

X

i=1

J^st(˜θ_i)W_iN_b,x(˜θ_i)T^β(˜θ_i)h

N_a,x(˜θ_i)T^α(˜θ_i)(^P2_γ=1^P3_c=1u^∗γc Nc(r)T^γ(˜θi)

i

+τ₁ Z

Qb⁴n

2

X

i=1

J^st(˜θ_i)W_iN_b,x(˜θ_i)T^β(˜θ_i)h

N_a,y(˜θ_i)T^α(˜θ_i)(^P2_γ=1^P3_c=1u^∗γc Nc(r)T^γ(˜θi)

i

+τ1

Z

Qb⁴n

2

X

i=1

J^st(˜θi)WiNb,x(˜θi)T^β(˜θi) h

(Vr(r,θ˜i)Na,r+Vs(r,θ˜i)Na,s)T^α(˜θi) i

+τ1

Z

Qb⁴n

2

X

i=1

J^st(˜θi)WiNb,x(˜θi)T^β(˜θi)(−1)^α

∆t Na(r)

=τ1 2

X

i=1

J^st(˜θi)Nb,x(˜θi)T^β(˜θi)Na,x(˜θi)T^α(˜θi) h

P3

c=1(u^∗1_c T¹(˜θi)+u^∗2_c T²(˜θi))R

Qb4 n Nc(r)dr

i

+τ1 2

X

i=1

J^st(˜θi)Nb,x(˜θi)T^β(˜θi)Na,y(˜θi)T^α(˜θi) h

P3

c=1(v^∗1_c T¹(˜θi)+v^∗2_c T²(˜θi))R

Qb4 n Nc(r)dr

i

+τ1 2

X

i=1

J^st(˜θi)Nb,x(˜θi)T^β(˜θi)T^α(˜θi) Z

Qb⁴n

−x_θ(r)ys(˜θi) +yθ(r)xs(˜θi) J^st(˜θ_i)

! Na,rdr

+τ1 2

X

i=1

J^st(˜θi)Nb,x(˜θi)T^β(˜θi)T^α(˜θi) Z

Qb⁴n

xθ(r)yr(˜θi)−yθ(r)xr(˜θi) J^st(˜θ_i)

! Na,sdr

+τ1 2

X

i=1

J^st(˜θi)Nb,x(˜θi)T^β(˜θi)(−1)^α

∆t Z

Qb⁴n

Na(r)dr

S2^e₃(a, b)

=τ₁

2

X

i=1

J^st(˜θ_i)N_b,x(˜θ_i)T^β(˜θ_i)N_a,x(˜θ_i)T^α(˜θ_i) [^P3_c=1(u^∗1_c T¹(˜θi)+u^∗2_c T²(˜θi))^det₆ ]

+τ₁

2

X

i=1

J^st(˜θ_i)N_b,x(˜θ_i)T^β(˜θ_i)N_a,y(˜θ_i)T^α(˜θ_i) [^P3_c=1(v_c^∗1T¹(˜θi)+v^∗2_c T²(˜θi))^det₆ ]

+τ1 2

X

i=1

Nb,x(˜θi)T^β(˜θi)T^α(˜θi) (A3Na,r+A4Na,s)

+τ1 2

X

i=1

J^st(˜θi)Nb,x(˜θi)T^β(˜θi)(−1)^α

∆t 2A

6

= τ1

6

2

X

i=1

J^st(˜θi)Nb,x(˜θi)T^β(˜θi)Na,x(˜θi)T^α(˜θi) [^P3c=1(u^∗1_c T¹(˜θi)+u^∗2_c T²(˜θi))]

+τ1

6

2

X

i=1

J^st(˜θi)Nb,x(˜θi)T^β(˜θi)Na,y(˜θi)T^α(˜θi) [^P3_c=1(v^∗1_c T¹(˜θi)+v_c^∗2T²(˜θi))]

+τ1

6

2

X

i=1

Nb,x(˜θi)T^β(˜θi)T^α(˜θi) (−^P∆xc₂ ys(˜θi)+P∆yc

2 xs(˜θi))Na,r

+τ1

6

2

X

i=1

Nb,x(˜θi)T^β(˜θi)T^α(˜θi) (^P∆xc₂ yr(˜θi)−P∆yc

2 xr(˜θi))Na,s

+τ1

6

2

X

i=1

J^st(˜θi)Nb,x(˜θi)T^β(˜θi)(−1)^α

∆t , where

A₃ =

−^P∆xc₂ ys(˜θi)R

Qb4 n

Nc(r)dr+P_∆yc

2 xs(˜θi)R

Qb4 n

Nc(r)dr

A4 =

P_∆xc

2 yr(˜θi)R

Qb4

n Nc(r)dr−P∆yc 2 xr(˜θi)R

Qb4 n Nc(r)dr

Similarly, in y-component : τ1

R

Q^e_nρ(^∂w_∂t^h +u^∂w_∂x^h +v^∂w_∂y^h)^∂p_∂y^hdQ (usingN_b,y instead ofN_b,x)

S2^e₆(a, b) =τ₁ Z

Q^e_n

N¯_a,t^α (x, t) + ¯N_a,x^α (x, t)u+ ¯N_a,y^α (x, t)vN¯_b,y^β (t) dxdt

S2^e₆(a, b)

= τ₁ 6

2

X

i=1

J^st(˜θ_i)N_b,y(˜θ_i)T^β(˜θ_i)N_a,x(˜θ_i)T^α(˜θ_i) [^P3_c=1(u^∗1_c T¹(˜θi)+u^∗2_c T²(˜θi))]

+ τ₁ 6

2

X

i=1

J^st(˜θ_i)N_b,y(˜θ_i)T^β(˜θ_i)N_a,y(˜θ_i)T^α(˜θ_i) [^P3_c=1(v^∗1_c T¹(˜θi)+v_c^∗2T²(˜θi))]

+ τ₁ 6

2

X

i=1

Nb,y(˜θi)T^β(˜θi)T^α(˜θi) (−^P∆xc₂ ys(˜θi)+P_∆yc

2 xs(˜θi))Na,r

+ τ1

6

2

X

i=1

Nb,y(˜θi)T^β(˜θi)T^α(˜θi) (^P∆xc₂ yr(˜θi)−P_∆yc

2 xr(˜θi))Na,s

+ τ1

6

2

X

i=1

J^st(˜θi)Nb,y(˜θi)T^β(˜θi)(−1)^α

∆t

c. Stab3,

Z

Q^e_n

τ2

ρ∇q^h· ∇p^hdQ (3.50)

S3^e₉(a, b)

= τ₂ ρ

Z

Q^e_n

∇N¯_a^α(x, t)· ∇N¯_b^β(x, t) dxdt

= τ₂ ρ

Z

Q^e_n

∇N_a^α(r, θ)·N_b^β(r, θ)J^st(θ) drdθ

= τ₂ ρ

Z

Q^e_n

N_a,x^α (θ)N_b,x^β (θ) +N_a,y^α (θ)N_b,y^β (θ)

J^st(θ) drdθ

= τ₂ ρ

Z

Q^e_n

N_a,x(θ)T^α(θ)N_b,x(θ)T^β(θ) +N_a,y(θ)T^α(θ)N_b,y(θ)T^β(θ)

J^st(θ) drdθ

= τ2

ρ

2

X

i=1

T^α(˜θ_i)T^β(˜θ_i)J^st(˜θ_i)W_i

N_a,x(˜θ_i)N_b,x(˜θ_i) +N_a,y(˜θ_i)N_b,y(˜θ_i)Z

Q⁴n

dr

= τ2

ρ

2

X

i=1

T^α(˜θ_i)T^β(˜θ_i)J^st(˜θ_i)W_i

N_a,x(˜θ_i)N_b,x(˜θ_i) +N_a,y(˜θ_i)N_b,y(˜θ_i)2A 2

S3^e₉(a, b)

= τ₂ 2ρ

2

X

i=1

T^α(˜θ_i)T^β(˜θ_i)J^st(˜θ_i)

N_a,x(˜θ_i)N_b,x(˜θ_i) +N_a,y(˜θ_i)N_b,y(˜θ_i)

= τ₂∆t 4ρ

2

X

i=1

T^α(˜θ_i)T^β(˜θ_i)Υ(˜θ_i)

N_a,x(˜θ_i)N_b,x(˜θ_i) +N_a,y(˜θ_i)N_b,y(˜θ_i) d. Stab4,

τ₂ Z

Q^e_n

∇q^h(∂u^h

∂t +u^h· ∇u^h) (3.51)

Stab4 is nonlinear term, it will be explained later.

e. Stab5,

τ₃ Z

Q^e_n

∇w^hρ∇u^hdQ (3.52)

S5^e₁(a, b) =τ₃ρ Z

Q^e_n

N¯_a,x^α (t) ¯N_b,x^β (t) dxdt

=τ₃ρ Z

Q^e_n

N_a,x^α (θ)N_b,x^β (θ)J^st(θ) drdθ

=τ3ρ

2

X

i=1

T^α(˜θi)T^β(˜θi)J^st(˜θi)WiNa,x(˜θi)N_b,x(˜θi) Z

Q⁴n

dr

=τ3ρ

2

X

i=1

T^α(˜θi)T^β(˜θi)J^st(˜θi)WiNa,x(˜θi)N_b,x(˜θi)2A 2

= τ3ρ 2

2

X

i=1

T^α(˜θi)T^β(˜θi)J^st(˜θi)Na,x(˜θi)N_b,x(˜θi)

= τ3ρ∆t 4

2

X

i=1

T^α(˜θi)T^β(˜θi)Υ(˜θi)Na,x(˜θi)N_b,x(˜θi) Similarly,

S5^e₂(a, b) =τ₃ρ Z

Q^e_n

N¯_a,x^α (t) ¯N_b,y^β (t) dxdt

= τ3ρ 2

2

X

i=1

T^α(˜θi)T^β(˜θi)J^st(˜θi)Na,x(˜θi)N_b,y(˜θi)

= τ3ρ∆t 4

2

X

i=1

T^α(˜θi)T^β(˜θi)Υ(˜θi)Na,x(˜θi)N_b,y(˜θi)

S5^e₄(a, b) =τ₃ρ Z

Q^e_n

N¯_a,y^α (t) ¯N_b,x^β (t) dxdt

= τ3ρ 2

2

X

i=1

T^α(˜θ_i)T^β(˜θ_i)J^st(˜θ_i)N_a,y(˜θ_i)N_b,x(˜θ_i)

= τ3ρ∆t 4

2

X

i=1

T^α(˜θ_i)T^β(˜θ_i)Υ(˜θ_i)N_a,y(˜θ_i)N_b,x(˜θ_i)

S5^e₅(a, b) =τ₃ρ Z

Q^e_n

N¯_a,y^α (t) ¯N_b,y^β (t) dxdt

= τ3ρ 2

2

X

i=1

T^α(˜θi)T^β(˜θi)J^st(˜θi)Na,y(˜θi)Nb,y(˜θi)

= τ3ρ∆t 4

2

X

i=1

T^α(˜θi)T^β(˜θi)Υ(˜θi)Na,y(˜θi)N_b,y(˜θi) 7)F₁ and F₂

Z

Q^e_n

ρw^hfdQ (3.53)

F₁^e(a) =ρ Z

Q^e_n

N¯_a^α(x, t)f_x dxdt

=ρ Z

Qb^e_n

N¯_a^α(x(r, θ), t(θ))f_xJ^st(θ)drdθ

=ρ Z

Qb⁴n

2

X

i=1

Na(r)T^α(˜θi)fxJ^st(˜θi)Wi dr

=ρfx 2

X

i=1

T^α(˜θi)J^st(˜θi)Wi

Z

Qb⁴n

Na(r)dr

=ρfx 2

X

i=1

T^α(˜θi)J^st(˜θi)Wi

2A 6

= ρfx

6

2

X

i=1

T^α(˜θi)J^st(˜θi)

= ρfx∆t 12

2

X

i=1

T^α(˜θi)Υ(˜θi)

Similarly forF₂^e(a),

F₂^e(a) =ρ Z

Q^e_n

N¯_a^α(x, t)fy dxdt

= ρf_y∆t 12

2

X

i=1

T^α(˜θ_i)Υ(˜θ_i)

8)F_s4

Z

Q^e_n

τ₂∇qfdQ (3.54)

F_S4^e ₃(a)

=τ₂ Z

Q^e_n

N¯_a,x^α (t)f_x+ ¯N_a,y^α (t)f_y dxdt

=τ2

"

Z

Qb^e_n

N¯_a,x^α (θ), t(θ))fx+ ¯N_a,y^α (θ), t(θ))fy

#

J^st(θ)drdθ

=τ2 2

X

i=1

Z

Qb⁴n

Na,x(˜θi)T^α(˜θi)fx+Na,y(˜θi)T^α(˜θi)fy

J^st(˜θi)Widr

=τ2 2

X

i=1

J^st(˜θi)Wi

h

Na,x(˜θi)T^α(˜θi)fx+Na,y(˜θi)T^α(˜θi)fy

iZ

Qb⁴n

dr

=τ2 2

X

i=1

J^st(˜θi)Wi

h

Na,x(˜θi)T^α(˜θi)fx+Na,y(˜θi)T^α(˜θi)fy

i2A 2

=τ2

2A 2

∆t 2

2

X

i=1

Υ(˜θi)Wi

h

Na,x(˜θi)T^α(˜θi)fx+Na,y(˜θi)T^α(˜θi)fy

i

= τ2∆t 4

2

X

i=1

Υ(˜θi)T^α(˜θi)h

Na,x(˜θi)fx+Na,y(˜θi)fy

i

9) Right-hand side of jump termFM, Z

Ω^e

ρw(u⁻) dΩ (3.55)

F_M^e₁(a) =ρ Z

Ω^e

Na(x)

3

X

b=1

u⁻_b N_b(x) dx

=ρ Z

Ω^e

N_a(r)

3

X

b=1

u⁻_b N_b(r)Υ drb

=ρ2A 24

3

X

b=1

(u⁻_b +ua)Υb

= ρΥb 24

3

X

b=1

(u⁻_b +u_a) Similarly,

F_M^e₂(a) =ρ Z

Ω^e

N_a(x)

3

X

b=1

v⁻_b N_b(x) dx

= ρΥb 24

3

X

b=1

(v⁻_b +v_a)

The nonlinear terms will be written in increment form, i.e., we take the solution u in the form

u=u^∗+ ∆u,

and substitute it into the corresponding terms and calculate the resulting contribution to the left-hand side of the system.

1) Convection Z

Q^e_n

ρw(u· ∇u) dQ

=ρ Z

Q^e_n

w((u^∗+ ∆u)· ∇(u^∗+ ∆u)) dQ

=ρ

"

Z

Q^e_n

w(u^∗· ∇u^∗) dQ+ Z

Q^e_n

w(u^∗· ∇∆u) dQ+ Z

Q^e_n

w(∆u· ∇u^∗) dQ

#

+ρ

"

Z

Q^e_n

w(∆u· ∇∆u) dQ

#

The first integral contributes to right-hand side since u^∗ is known. The second integral contributes to the left-hand side and we will call it, convec-tion term (C(U^∗)). The third integral contributes to the left-hand side and

we will call it convection derivative term (CC(U^∗)). The last integral is dropped since it is quadratic in ∆u[31].

a. Convection term (C(U^∗))

ρ Z

Q^e_n

w(u^∗· ∇∆u) dQ (3.56)

C₁^e(U^∗)(a, b)

=ρ Z

Q^e_n

N¯_a^α(x, t)(u^∗· ∇N¯_b^β(t)) dxdt

=ρ Z

Q^e_n

N¯_a^α(x, t)h

u^∗N¯_b,x^β (t) +v^∗N¯_b,y^β (t)i dxdt

=ρ Z

Qb^e_n

N_a^α(x(r, θ), t(θ)) h

u^∗N_b,x^β (t(θ)) +v^∗N_b,y^β (t(θ)) i

J^st(θ) drdt

=ρ Z

Qb^e_n

Na(r)T^α(θ) h

u^∗Nb,x(θ)T^β(θ) +v^∗Nb,y(θ)T(θ) i

J^st(θ) drdθ

=ρ

2

X

i=1

J^st(˜θ_i)W_iT^α(˜θ_i)T^β(˜θ_i)N_b,x(˜θ_i) Z

Qb⁴n

N_a(r)^P2_γ=1^P3_c=1u^∗γc Nc(r)T^γ(˜θi) dr

+ρ

2

X

i=1

J^st(˜θ_i)W_iT^α(˜θ_i)T^β(˜θ_i)N_b,y(˜θ_i) Z

Qb⁴n

N_a(r)^P2_γ=1^P3_c=1vc^∗γNc(r)T^γ(˜θi) dr

=ρ

2

X

i=1

J^st(˜θ_i)W_iT^α(˜θ_i)T^β(˜θ_i)N_b,x(˜θ_i)

" ₃ X

c=1

(u^∗1_c T¹(˜θ_i) +u^∗2_c T²(˜θ_i)) Z

Qb⁴n

N_c(r)N_a(r) dr

#

+ρ

2

X

i=1

J^st(˜θ_i)W_iT^α(˜θ_i)T^β(˜θ_i)N_b,y(˜θ_i)

" ₃ X

c=1

(v_c^∗1T¹(˜θ_i) +v^∗2_c T²(˜θ_i)) Z

Qb⁴n

N_c(r)N_a(r) dr

#

=ρ

2

X

i=1

J^st(˜θ_i)W_iT^α(˜θ_i)T^β(˜θ_i)N_b,x(˜θ_i)2A

24 (A₅+A₆)

C₁^e(U^∗)(a, b) =ρ

2

X

i=1

∆t

2 Υ(˜θi)WiT^α(˜θi)T^β(˜θi)N_b,x(˜θi)2A

24 (A5+A6)

=ρ∆t 48

2

X

i=1

Υ(˜θi)T^α(˜θi)T^β(˜θi)N_b,x(˜θi) (A5+A6) where

A₅ =

" ₃ X

c=1

(u^∗1_c +u^∗1_a )T¹(˜θ_i) + (u^∗2_c +u^∗2_a )T²(˜θ_i)

#

A₆ =

" ₃ X

c=1

(v_c^∗1+v_a^∗1)T¹(˜θ_i) + (v_c^∗2+v_a^∗2)T²(˜θ_i)

#

(3.57) C₅^e(U^∗)(a, b) =C₁^e(U^∗)(a, b)

b. Convection derivative (CC(U^∗)), ρ

Z

Q^e_n

w(∆u· ∇u^∗) dQ (3.58)

ρ Z

Q^e_n

w(∆u· ∇u^∗) dQ=















CC1(U^∗) =ρR

Q^e_nw1(∆u^∂u_∂x^∗) dQ CC2(U^∗) =ρR

Q^e_nw1(∆v^∂u_∂y^∗) dQ CC₄(U^∗) =ρR

Q^e_nw₂(∆u^∂v_∂x^∗) dQ CC₁(U^∗) =ρR

Q^e_nw₂(∆v^∂v_∂y^∗) dQ CC₁^e(U^∗)(a, b)

=ρ Z

Q^e_n

N¯_a^α(x, t) ¯N_b^β(x, t)





2

X

γ=1 3

X

c=1

u^∗γ_c N¯_c,x^γ (t)



 dxdt

=ρ Z

Qb^e_n

N_a^α(r, θ)N_b^β(r, θ)





2

X

γ=1 3

X

c=1

u^∗γ_c N_c,x^γ (θ)



J^st(θ) drdθ

=ρ Z

Q⁴n

2

X

i=1

J^st(˜θ_i)W_iN_a(r)T^α(˜θ_i)N_b(r)T^β(˜θ_i)





2

X

γ=1 3

X

c=1

u^∗γ_c N_c,x(˜θ_i)T^γ(θ)



 dr

=ρ

2

X

i=1

J^st(˜θi)WiT^α(˜θi)T^β(˜θi)





2

X

γ=1 3

X

c=1

u^∗γ_c Nc,x(˜θi)T^γ(˜θi)



 Z

Q^e_n

Na(r)N_b(r) dr

CC₁^e(U^∗)(a, b)

=

























 ρ

2

X

i=1

∆t

2 Υ(˜θ_i)W_iT^α(˜θ_i)T^β(˜θ_i)





2

X

γ=1 3

X

c=1

u^∗γ_c N_c,x(˜θ_i)T^γ(˜θ_i)



 2A

12, ifa=b.

ρ

2

X

i=1

∆t

2 Υ(˜θi)WiT^α(˜θi)T^β(˜θi)





2

X

γ=1 3

X

c=1

u^∗γ_c Nc,x(˜θi)T^γ(˜θi)



 2A

24, otherwise.

=























ρ∆t 24

2

X

i=1

Υ(˜θ_i)T^α(˜θ_i)T^β(˜θ_i)

3

X

c=1

(u^∗1_c T¹(˜θ_i) +u^∗2_c T²(˜θ_i))N_c,x(˜θ_i)

! , ifa=b.

ρ∆t 48

2

X

i=1

Υ(˜θi)T^α(˜θi)T^β(˜θi)

3

X

c=1

(u^∗1_c T¹(˜θi) +u^∗2_c T²(˜θi)Nc,x(˜θi)

! , otherwise.

And similarly for other (CC(U^∗))

CC₂^e(U^∗)(a, b)

=ρ Z

Q^e_n

N¯_a^α(x, t) ¯N_b^β(x, t)





2

X

γ=1 3

X

c=1

u^∗γ_c N¯_c,y^γ (t)



 dxdt

=























ρ∆t 24

2

X

i=1

Υ(˜θ_i)T^α(˜θ_i)T^β(˜θ_i)

3

X

c=1

(u^∗1_c T¹(˜θ_i) +u^∗2_c T²(˜θ_i))N_c,y(˜θ_i)

! , ifa=b.

ρ∆t 48

2

X

i=1

Υ(˜θi)T^α(˜θi)T^β(˜θi)

3

X

c=1

(u^∗1_c T¹(˜θi) +u^∗2_c T²(˜θi)Nc,y(˜θi)

! , otherwise.

CC₄^e(U^∗)(a, b)

=ρ Z

Q^e_n

N¯_a^α(x, t) ¯N_b^β(x, t)





2

X

γ=1 3

X

c=1

v_c^∗γN¯_c,y^γ (t)



 dxdt

=























ρ∆t 24

2

X

i=1

Υ(˜θ_i)T^α(˜θ_i)T^β(˜θ_i)

3

X

c=1

(v_c^∗1T¹(˜θ_i) +v_c^∗2T²(˜θ_i))N_c,x(˜θ_i)

! , ifa=b.

ρ∆t 48

2

X

i=1

Υ(˜θi)T^α(˜θi)T^β(˜θi)

3

X

c=1

(v_c^∗1T¹(˜θi) +v_c^∗2T²(˜θi)Nc,x(˜θi)

! , otherwise.

CC₅^e(U^∗)(a, b)

=ρ Z

Q^e_n

N¯_a^α(x, t) ¯N_b^β(x, t)





2

X

γ=1 3

X

c=1

v^∗γ_c N¯_c,y^γ (t)



 dxdt

=























ρ∆t 24

2

X

i=1

Υ(˜θ_i)T^α(˜θ_i)T^β(˜θ_i)

3

X

c=1

(v_c^∗1T¹(˜θ_i) +v^∗2_c T²(˜θ_i))N_c,y(˜θ_i)

! , ifa=b

ρ∆t 48

2

X

i=1

Υ(˜θi)T^α(˜θi)T^β(˜θi)

3

X

c=1

(v_c^∗1T¹(˜θi) +v^∗2_c T²(˜θi)Nc,y(˜θi)

! , otherwise

2) Stab4 τ₂

Z

Q^e_n

∇q ∂u

∂t +u· ∇u

dQ

=τ2

"

Z

Q^e_n

∇q ∂∆u

∂t +u^∗· ∇∆u

dQ+ Z

Q^e_n

∇q

∂∆u

∂t + ∆u· ∇u^∗

dQ

#

+τ2

"

Z

Q^e_n

∇q ∂∆u

∂t + ∆u· ∇∆u

dQ+ Z

Q^e_n

∇q ∂u^∗

∂t +u^∗· ∇u^∗

dQ

#

The first integral contributes to the left-hand side and we will call it Stab4 term (S4(U^∗)). The second integral is dropped because of convergence rea-sons. The third integral is also dropped since it is quadratic in ∆u. The last integral contributes to the right-hand side since u^∗ is known.

Stab4 term (S4(U^∗)), τ₂

Z

Q^e_n

∇q ∂∆u

∂t +u^∗· ∇∆u

dQ (3.59)

Stab4 is similar with stab2 term, only we reverseawithband α withβ. S4^e₇(U^∗)(a, b)

=τ1

Z

Q^e_n

hN¯_b,t^β (x, t) + ¯N_b,x^β (x, t)u+ ¯N_b,y^β (x, t)v

iN¯_a,x^α (t) dxdt

= τ₁ 6

2

X

i=1

J^st(˜θi)Na,x(˜θi)T^α(˜θi)Nb,x(˜θi)T^β(˜θi) [^P3c=1(u^∗1_c T¹(˜θi)+u^∗2_c T²(˜θi))]

+ τ₁ 6

2

X

i=1

J^st(˜θi)Na,x(˜θi)T^α(˜θi)Nb,y(˜θi)T^β(˜θi) [^P3c=1(v_c^∗1T¹(˜θi)+v^∗2_c T²(˜θi))]

+ τ1

6

2

X

i=1

Na,x(˜θi)T^α(˜θi)T^β(˜θi) (−^P∆xc₂ ys(˜θi)+P_∆yc

2 xs(˜θi))Nb,r

+ τ1

6

2

X

i=1

Na,x(˜θi)T^α(˜θi)T^β(˜θi) (^P∆xc₂ yr(˜θi)−P∆yc

2 xr(˜θi))Nb,s

+ τ1

6

2

X

i=1

J^st(˜θi)Na,x(˜θi)T^α(˜θi)(−1)^β

∆t

S4^e₈(U^∗)(a, b)

=τ₁ Z

Q^e_n

hN¯_b,t^β (x, t) + ¯N_b,x^β (x, t)u+ ¯N_b,y^β (x, t)vi

N¯_a,y^α (t) dxdt

= τ1

6

2

X

i=1

J^st(˜θi)Na,y(˜θi)T^α(˜θi)N_b,x(˜θi)T^β(˜θi) [^P3_c=1(u^∗1_c T¹(˜θi)+u^∗2_c T²(˜θi))]

+τ1

6

2

X

i=1

J^st(˜θi)Na,y(˜θi)T^α(˜θi)N_b,y(˜θi)T^β(˜θi) [^P3_c=1(v_c^∗1T¹(˜θi)+v^∗2_c T²(˜θi))]

+τ1

6

2

X

i=1

Na,y(˜θi)T^α(˜θi)T^β(˜θi) (−^P∆xc₂ ys(˜θi)+P_∆yc

2 xs(˜θi))N_b,r

+τ1

6

2

X

i=1

Na,y(˜θi)T^α(˜θi)T^β(˜θi) (^P∆xc₂ yr(˜θi)−P_∆yc

2 xr(˜θi))N_b,s +τ1

6

2

X

i=1

J^st(˜θi)Na,y(˜θi)T^α(˜θi)(−1)^β

∆t 3) Stab1

Z

Q^e_n

τ1

1 ρ

ρ

∂w

∂t +u· ∇w

ρ ∂u

∂t +u· ∇u

dQ (3.60)

Z

Q^e_n

τ1

1 ρ

ρ

∂w

∂t +u· ∇w

ρ ∂u

∂t +u· ∇u

dQ

= Z

Q^e_n

τ₁ρ ∂w

∂t +u^∗· ∇w ∂(u^∗+ ∆u)

∂t + (u^∗+ ∆u)· ∇(u^∗+ ∆u)

dQ Dropping ∆uin the first bracket improves the convergence rate [31].

=τ1ρ

"

Z

Q^e_n

∂w

∂t +u^∗· ∇w ∂u^∗

∂t +u^∗· ∇u^∗

dQ

#

+τ1ρ

"

Z

Q^e_n

∂w

∂t +u^∗· ∇w ∂∆u

∂t +u^∗· ∇∆u

dQ

#

+τ1ρ

"

Z

Q^e_n

∂w

∂t +u^∗· ∇w

(∆u· ∇u^∗) dQ

#

+τ1ρ

"

Z

Q^e_n

∂w

∂t +u^∗· ∇w

(∆u· ∇∆u) dQ

#

The first integral contributes to the the right-hand side since u^∗ is known.

The second integral contributes to the left-hand side and we will call it stab1 term (S1(U^∗)). The third integral contributes to the left-hand side and we will call it stab1 derivative term (SS1(U^∗)). The last integral is dropped since it is quadratic in ∆u. Since stab1 term (second integral) is compli-cated, we divide it into stab1a term, stab1b term, stab1c term and stab1d term,

a. stab1a

τ₁ρ Z

Q^e_n

∂w

∂t

∂∆u

∂t dQ (3.61)

S1a^e₁(a, b)

=τ₁ρ Z

Q^e_n

N¯_a,t^α (x, t) ¯N_b,t^β(x, t) dxdt

=τ₁ρ Z

Qb^e_n

N¯_a,t^α (x(r, θ), t(θ)) ¯N_b,t^β (x(r, θ), t(θ))J^st(θ) drdθ

=τ₁ρ Z

Qb⁴n

2

X

i=1

N_a,t^α (r,θ˜_i)N_b,t^β (r,θ˜_i)J^st(˜θ_i)W_i dr

=τ₁ρ

2

X

i=1

J^st(˜θ_i)W_i Z

Qb⁴n

N_a,t^α (r,θ˜_i)N_b,t^β(r,θ˜_i)dr

=τ₁ρ

2

X

i=1

J^st(˜θ_i)W_i Z

Qb⁴n

V_r(r,θ˜_i)N_a,r+V_s(r,θ˜_i)N_a,s

T^α(˜θ_i) + (−1)^α

∆t N_a(r)

Vr(r,θ˜i)Nb,r+Vs(r,θ˜i)Nb,s

T^β(˜θi) +(−1)^β

∆t Nb(r)

dr S1a^e₅(a, b) =S1a^e₁(a, b)

b. stab1b

τ₁ρ Z

Q^e_n

∂w

∂tu^∗· ∇∆udQ (3.62) S1b^e₁(a, b)

=τ₁ρ Z

Q^e_n

N¯_a,t^α (x, t)u^∗· ∇barN_b^β(x, t) dxdt

=τ₁ρ Z

Qb^e_n

N¯_a,t^α (x(r, θ), t(θ)) ¯N_b,t^β (x(r, θ), t(θ))J^st(θ) drdθ

=τ1ρ Z

Qb^e_n

N_a,t^α (r, θ)

u^∗·N_b,x^β (r, θ) +v^∗·N_b,y^β (r, θ)

J^st(θ) drdθ

=τ₁ρ

2

X

i=1

J^st(˜θ_i)W_i Z

Qb⁴n

V_r(r,θ˜_i)N_a,r+V_s(r,θ˜_i)N_a,s

T^α(˜θ_i) + (−1)^α

∆t N_a(r)





2

X

γ=1 3

X

c=1

u^∗γ_c Nc(r)T^γ(˜θi)N_b,x(˜θi)T^β(˜θi) +

2

X

γ=1 3

X

c=1

v_c^∗γNc(r)T^γ(˜θi)N_b,y(˜θi)T^β(˜θi)



dr

S1b^e₅(a, b) =S1b^e₁(a, b) c. stab1c

τ₁ρ Z

Q^e_n

u^∗· ∇w∂∆u

∂t dQ (3.63)

S1c^e₁(a, b)

=τ1ρ Z

Q^e_n

N¯_b,t^β(x, t)u^∗· ∇N¯_a^α(x, t) dxdt

=τ1ρ Z

Qb^e_n

N¯_b,t^β(x(r, θ), t(θ)) ¯N_a,t^α (x(r, θ), t(θ))J^st(θ) drdθ

=τ1ρ Z

Qb^e_n

N_b,t^β(r, θ) u^∗·N_a,x^α (r, θ) +v^∗·N_a,y^α (r, θ)

J^st(θ) drdθ

=τ₁ρ

2

X

i=1

J^st(˜θ_i)W_i Z

Qb⁴n

V_r(r,θ˜_i)N_b,r+V_s(r,θ˜_i)N_b,s

T^β(˜θ_i) +(−1)^β

∆t N_b(r)





2

X

γ=1 3

X

c=1

u^∗γ_c N_c(r)T^γ(˜θ_i)N_a,x(˜θ_i)T^α(˜θ_i) +

2

X

γ=1 3

X

c=1

v_c^∗γN_c(r)T^γ(˜θ_i)N_a,y(˜θ_i)T^α(˜θ_i)



dr

S1c^e₅(a, b) =S1c^e₁(a, b) d. stab1d

τ₁ρ Z

Q^e_n

(u^∗· ∇w) (u^∗· ∇∆u) dQ (3.64)

S1d^e₁(a, b)

=τ₁ρ Z

Q^e_n

u^∗· ∇N¯_a^α(x, t)

u^∗· ∇N¯_b^β(x, t) dxdt

=τ₁ρ Z

Q^e_n

u^∗N¯_a,x^α (x, t) +v^∗N¯_a,y^α (x, t)

u^∗N¯_b,x^β (x, t) +v^∗N¯_b,y^β (x, t) dxdt

=τ₁ρ Z

Qb^e_n

(u^∗N_a,x(θ)T^α(θ) +v^∗N_a,y(θ)T^α(θ))

u^∗N_b,x(θ)T^β(θ) +v^∗N_b,y(θ)T^β(θ)

J^st(θ) drdθ

=τ₁ρ

2

X

i=1

T^α(˜θ_i)T^β(˜θ_i)J^st(˜θ_i)W_i Z

Qb⁴n

(u^∗Na,x(θ) +v^∗Na,y(θ)) (u^∗N_b,x(θ) +v^∗N_b,y(θ)) dr S1d^e₅(a, b) =S1d^e₁(a, b)

Stab1 derivative term (SS1(U^∗)) can be written as, τ1ρ

Z

Q^e_n

∂w

∂t +u^∗· ∇w

(∆u^∗· ∇u^∗) dQ

(3.65) or,

SS1^e₁(a, b) =τ₁ρ Z

Q^e_n

∂w

∂t +u^∗w_1,x+v^∗w_1,y

∆u u^∗_,xdQ

SS1^e₂(a, b) =τ₁ρ Z

Q^e_n

∂w

∂t +u^∗w_1,x+v^∗w_1,y

∆v u^∗_,y dQ

SS1^e₄(a, b) =τ₁ρ Z

Q^e_n

∂w

∂t +u^∗w_1,x+v^∗w_1,y

∆u v^∗_,xdQ

SS1^e₅(a, b) =τ₁ρ Z

Q^e_n

∂w

∂t +u^∗w_1,x+v^∗w_1,y

∆v v_,ydQ

with SS1^e₁(a, b)

=τ₁ρ Z

Q^e_n

N¯_a,t^α (x, t) +u^∗N¯_a,x^α (x, t) +v^∗N¯_a,y^α (x, t)N¯_b^β(x, t)

3

X

c=1

(u^∗1_c T¹+u^∗2_c T²) ¯N_c,x(θ) dQ

=τ₁ρ Z

Qb^e_n

(N_a,t(r, θ)T^α(θ) +u^∗N_a,x(θ)T^α(θ) +v^∗N_a,y(θ)T^α(θ))N_b(r)T^β(θ)

3

X

c=1

(u^∗1_c T¹+u^∗2_c T²)N_c,x(θ)J^st(θ)drdθ

SS1^e₂(a, b)

=τ1ρ Z

Q^e_n

N¯_a,t^α (x, t) +u^∗N¯_a,x^α (x, t) +v^∗N¯_a,y^α (x, t)N¯_b^β(x, t)

3

X

c=1

(u^∗1_c T¹+u^∗2_c T²) ¯N_c,y(θ) dQ

=τ1ρ Z

Qb^e_n

(Na,t(r, θ)T^α(θ) +u^∗Na,x(θ)T^α(θ) +v^∗Na,y(θ)T^α(θ))Nb(r)T^β(θ)

3

X

c=1

(u^∗1_c T¹+u^∗2_c T²)N_c,y(θ)J^st(θ)drdθ

SS1^e₄(a, b)

=τ1ρ Z

Q^e_n

N¯_a,t^α (x, t) +u^∗N¯_a,x^α (x, t) +v^∗N¯_a,y^α (x, t)N¯_b^β(x, t)

3

X

c=1

(v_c^∗1T¹+v_c^∗2T²) ¯N_c,x(θ) dQ

=τ1ρ Z

Qb^e_n

(Na,t(r, θ)T^α(θ) +u^∗Na,x(θ)T^α(θ) +v^∗Na,y(θ)T^α(θ))Nb(r)T^β(θ)

3

X

c=1

(v_c^∗1T¹+v_c^∗2T²)N_c,x(θ)J^st(θ)drdθ

SS1^e₅(a, b)

=τ1ρ Z

Q^e_n

N¯_a,t^α (x, t) +u^∗N¯_a,x^α (x, t) +v^∗N¯_a,y^α (x, t)N¯_b^β(x, t)

3

X

c=1

(v^∗1_c T¹+v_c^∗2T²) ¯Nc,y(θ) dQ

=τ1ρ Z

Qb^e_n

(Na,t(r, θ)T^α(θ) +u^∗Na,x(θ)T^α(θ) +v^∗Na,y(θ)T^α(θ))N_b(r)T^β(θ)

3

X

c=1

(v^∗1_c T¹+v_c^∗2T²)N_c,y(θ)J^st(θ)drdθ

There is another term in the right-hand sideF(3.38) which is changing for each Newton iteration, since it depend onu^∗,

Z

Q^e_n

τ₁ρ ∂w^h

∂t +u^∗· ∇w

fdQ (3.66)

F_SS1^e ₁(a)

=τ₁ρ Z

Q^e_n

N¯_a,t^α (x, t) +u^∗∇N¯_a^α(x, t)

f_xdxdt

= τ1ρfx

6

2

X

i=1

h

T^α(˜θi) (−^P∆xc₂ ys(˜θi)+P∆yc

2 xs(˜θi))Na,r

i

+ τ1ρfx

6

2

X

i=1

T^α(˜θi) (^P∆xc₂ yr(˜θi)−P∆yc

2 xr(˜θi))Na,s+ Υ(˜θi)(−1)^α 2

+ τ1ρfx∆t 12

2

X

i=1

Υ(˜θi)T^α(˜θi)Na,x(˜θi)^P3_c=1(u^∗1_c T¹(˜θi)+u^∗2_c T²(˜θi))

+ τ1ρfx∆t 12

2

X

i=1

Υ(˜θi)T^α(˜θi)Na,y(˜θi)^P3_c=1(v^∗1_c T¹(˜θi)+v_c^∗2T²(˜θi))

Similarly for F_SS1^e ₂(a), F_SS1^e ₂(a)

=τ1ρ Z

Q^e_n

N¯_a,t^α (x, t) +u^∗∇N¯_a^α(x, t)

fy dxdt

= τ1ρfy

6

2

X

i=1

h

T^α(˜θ_i) (−^P∆xc₂ ys(˜θi)+P_∆yc

2 xs(˜θi))N_a,r i + τ₁ρf_y

6

2

X

i=1

T^α(˜θ_i) (^P∆xc₂ yr(˜θi)−P_∆yc

2 xr(˜θi))N_a,s+ Υ(˜θ_i)(−1)^α 2

+ τ₁ρf_y∆t 12

2

X

i=1

Υ(˜θ_i)T^α(˜θ_i)N_a,x(˜θ_i)^P3_c=1(u^∗1_c T¹(˜θi)+u^∗2_c T²(˜θi))

+ τ₁ρf_y∆t 12

2

X

i=1

Υ(˜θ_i)T^α(˜θ_i)N_a,y(˜θ_i)^P3_c=1(v_c^∗1T¹(˜θi)+v^∗2_c T²(˜θi))

3.2.7 Appendix: The Core of Stabilization

Solving differential equations using numerical methods often encounters some severe problems such as oscillation, singular matrix, etc. In such cases, stabilization is needed to get satisfactory result. Some of the circumstances in which such problems occur are convection-dominated problems or viola-tion of the Babuska-Brezzi condiviola-tion, which, in mixed formulaviola-tion such as Navier-Stokes equation, may be caused by using equal order of basis func-tions for velocity and pressure. The underlying idea about the stabilization can be found in [32]. In this appendix, we discuss briefly about the core idea of stabilization for time-dependent advection-diffusion 1D problem in usual finite element method (not space-time finite element method). For more details, we refer to [47].

In order to understand the core idea of stabilization, we use simple equa-tion, i.e., time-dependent advection-diffusion equation 1D,

∂φ

∂t +u∂φ

∂x−ν∂²φ

∂x² =f.

Assume that we have some suitably-defined finite-dimensional function spaces for trial function S^h and test function V^h. Weak formulation for time-dependent advection-diffusion equation 1D can be written as: find φ^h ∈S^h

such that∀w^h ∈V^h: Z

Ω

w^h∂φ^h

∂t dΩ + Z

Ω

w^hu^h∂φ^h

∂x dΩ

| {z } +∂w^h

∂x ν∂φ^h

∂x dΩ

| {z }

= Z

Γ

w^hh^h+ Z

Ω

w^hf^h.

Adv Diff

We compare the convection/advection term with diffusion term (in pare-metrical domain) to know whether the problem convection-dominated,

Adv Diff =

Z 1

−1

Nau^h2 h

∂N_b

∂ξ Υdξ ν

Z 1

−1

2 h

∂Na

∂ξ ν2 h

∂Nb

∂ξ Υdξ

= u^h

Z 1

−1

Na

∂N_b

∂ξ dξ ν_h²

Z 1

−1

∂Na

∂ξ

∂Nb

∂ξ dξ .

In this case, Υ = ^h₂ is the Jacobian of the transformation from physical domain to parametrical domain. The integrand of these two integrals have no dimension. So, to know which term is dominant, we compare only the dimension part of the formulation,

u^h

2

hν = u^h

2ν, (3.67)

which represents the element Peclet numberP e_h or element Reynold num-bersRehin Navier-Stokes equation. The problem is convection-dominated if P e_h 1 and non convection-dominated ifP e_h≈1. Note that the criterion based on the element Peclet number (element Reynold numbers) and not based on the global Peclet number (global Reynold numbers). It is closely related to the element length h in the formulation. We can make the mesh finer and finer (very small h) such that the stabilization is not necessary.

However, very smallh is numerically impractical.

Consider the stabilized formulation of time-dependent advection-diffusion equation 1D,

Z

Ω

w^h∂φ^h

∂t dΩ + Z

Ω

w^hu^h∂φ^h

∂x dΩ

| {z } +∂w^h

∂x ν∂φ^h

∂x dΩ

| {z } +

Adv Diff

nel

X

e=1

Z

Ω^e

τ u^h∂w^h

∂x ∂φ

∂t +u∂φ

∂x −ν∂²φ

∂x² −f

dΩ

| {z }

= Z

Γh

w^hh^hdΓ +f Z

Ω

w^hf dΩ.

Stab

Note that the formulation is still consistent since the stabilization term is residual based formulation. By comparing ”Adv”, ”Diff”, and ”Stab”

terms, we can get some useful information,

• If we compare ”Diff” and ”Stab”, we get the dimension ofτ. We obtain thatτ(u^h)² has the same dimension withν. It means we can consider τ(u^h)² as numerical viscosity ˜ν and stabilization can be viewed as adding the numerical viscosity. We have to add the numerical viscosity as needed, otherwise, we can lost the accuracy of the solution.

• If we compare ”Adv” and ”Diff”, we get element Peclet number P e_h as in (3.67).

P e_h= u^hh 2ν ,

• If we compare ”Adv” and ”Stab”, we get numerical element Peclet number,

P e˜ _h= u^hh

2˜ν = u^hh

2τ(u^h)² = h 2u^h

1 τ.

From this information, we can avoid the numerical difficulties caused by convection-dominated problem (assume P e_h 1) by setting the stabiliza-tion parameter τ with

P e˜ _h≈1, h

2u^h 1 τ ≈1, τ = h

2u^h. (3.68)

Note that (3.68) is one of the several selections that we take such that we add stabilization by adding the numerical viscosity as needed. In the higher dimensional case, we can represent (3.68) as

τ = h 2||u^h||, where ||u^h|| is the magnitude ofu^h.

3.2.8 Appendix: Stabilization in DSD/SST-SUPS formula-tion

In this appendix, we discuss briefly about the stabilization that we use in DSD/SST-SUPS formulation in (3.35) (refer to [47] for more details).

The formula and the role for each stabilization term are as follows,

• Streamline-Upwind/Petrov-Galerkin (SUPG)

(n_el)n

X

e=1

Z

Q^e_n

τ₁ ∂w^h

∂t +u^h· ∇w^h

·

ρ ∂u^h

∂t +u^h· ∇u^h−f^h

− ∇ ·σ(p^h,u^h)

| {z }

dQ residual

To avoid numerical instability that is caused by dominating convection.

• Pressure Stabilizing/Petrov-Galerkin (PSPG)

(nel)n

X

e=1

Z

Q^e_n

1 ρτ₂∇q^h

ρ

∂u^h

∂t +u^h· ∇u^h−f^h

− ∇ ·σ(p^h,u^h)

| {z }

dQ residual

To avoid pressure oscillations that are generated due to using equal order of basis functions for velocity and pressure with the purpose of simplifying the implementation. Moreover, without PSPG, we have zero block matrix K9 in (3.36), which leads to an ill-conditioned ma-trix.

• Least Square on Incompressible Constraint (LSIC)

(nel)n

X

e=1

Z

Q^e_n

τ3(∇ ·w^h)ρ(∇ ·u^h)

| {z } dQ residual

To reduce the divergence error that can occur in low or high Reynold numbers. To handle problems with very high Reynold numbers, be-side LSIC, we also need additional stabilization terms as in DSD/SST-VMST (DSD version with the variational multiscale turbulence model).

Note that when adding the stabilizations to the formulation, we have to make sure that the resulting formulation is still consistent, i.e., the exact solution still satisfies the resulting formulation. This is one reason we use a formula with the residual as the factor, i.e., momentum equation is used as a factor in SUPG and PSPG while continuity equation is used as a factor in LSIC.

We briefly explain the calculation of the stabilization parameters τ. In the stabilized formulation, the appropriate stabilization parameter τ plays an important role in the accuracy of the formulation [14]. In the DSD/SST formulation (3.35), we use three stabilization parameters τ₁, τ₂, τ₃ which are τSU P G, τP SP G, τLSIC, respectively. The definition of each stabilization parameter can be found in [14]. Here, we show the way to compute it in the S-DSD/SST context.

• Stabilization parameter τ₁ orτ_{SU P G} τ1=

1

τ_{SU GN12}² + 1 τ_{SU GN}² ₃

−¹

2

, (3.69)

where τ_{SU GN12}

=

2

X

α=1 3

X

a=1

∂N¯_a^α(x, t)

∂t +u^h· ∇N¯_a^α(x, t)

!−1

=

2

X

α=1 3

X

a=1

∂N¯_a^α(x, t)

∂t +uN¯_a,x^α (x, t) +vN¯_a,y^α (x, t)

!⁻¹

=

2

X

α=1 3

X

a=1

∂N¯_a^α(x(r, θ), t(θ))

∂t +uNa,x(θ)T^α(θ) +vNa,y(θ)T^α(θ)

!⁻¹

with u=

3

X

c=1

u^∗1_c T¹(θ) +u^∗2_c T²(θ) N_c(r), v=

3

X

c=1

v_c^∗1T¹(θ) +v_c^∗2T²(θ) Nc(r),

and ¯N_a,t^α (x, t),N¯_a,x^α (x, t),N¯_a,y^α (x, t) as in (3.43).

τ_{SU GN3} = h²_RGN 4ν

with

hRGN = 2

2

X

α=1 3

X

a=1

z· ∇N¯_a^α(x, t)

!−1

= 2

2

X

α=1 3

X

a=1

z₁N¯_a,x^α (x, t) +z₂N¯_a,y^α (x, t)

!−1

= 2

2

X

α=1 3

X

a=1

|z₁N_a,x(θ)T^α(θ) +z₂N_a,y(θ)T^α(θ)|

!⁻¹

and the solution gradient unit vector is defined as z= ∇ kuk

k ∇ kukk

Note that this quantity is not constant. Practically, we can evaluate it either using the center of each element or at the quadrature points.

In our simulation, we use the first approach.

• Stabilization parameter τ2 orτ_{P SP G}

τ₂ =τ₁ (3.70)

• Stabilization parameter τ3 orτLSIC

τ₃ =τ₁ kuk² (3.71)

This quantity is also not constant, we evaluate it using the center of each element.

3.2.9 Appendix: GMRES

GMRES is a projection method computing approximate solutionxn to the systemAx=b as minimizers of the residual norm

kr_nk=kb−Axnk, in the Krylov subspace

K_n= span [r₀, Ar₀, A²r₀, ..., Aⁿ⁻¹r₀],

upon finding an orthonormal basis{q₁, q2, ..., qn}of Kn and denoting by ˜Qn

the matrix with columns q1, q2, ..., qn. The above minimization is equivalent to the minimization ofkH˜_ny− kr₀ke₀kwheree₀= (1,0, ...,0) and ˜H_n is the Hessenberg matrix satisfying the similarity transformation

AQ˜_n= ˜Q_n+1H˜_n, and xnis obtained from xn=x0+ ˜Qny.

Note that the corresponding problem in our case is to approximate solu-tion ∆U^l to the system (3.37) as minimizers of the residual norm

kr^lk=kE^l−J^l∆U^lk, in the Krylov subspace

K_l= span [r₀,Jr₀,J²r₀, ...,J^l−1r₀], where E^l is the right-hand side of (3.37).

The outline of the GMRES algorithm can be written as follows [35]

Algorithm : GMRES

Given initial value x₀, we have initial residual isr₀=b−Ax₀. 1 q₁ = _kr^r0

0k

2 for j= 1,2, ..., m 3 computevj =Aqj

4 fori= 1, ..., j 5 h_i,j = (v_j, q_i) 6 vj =vj −hi,jqi

7 end

8 computeh_j+1,j =kv_jk₂ and v_j+1= _h^vj

j+1,j

9 end

10 define ˜Q_m := [q₁, ..., q_m],H˜_m={h_ij}1≤i≤m+1;1≤j≤m

11 compute y_m that minimizes kH˜_my− kr₀ke₀kand x_m=x₀+ ˜Q_my_m Step 2-9 of the algorithm is well-known with Arnoldi Iteration (modified Gram-Schmidt). Arnoldi iteration as orthogonal projection onto Krylov subspaces is an algorithm for building an orthogonal basis of the Krylov subspace Kn. This algorithm is based on the similarity transformation

AQ˜= ˜QH,˜