c 2015 Institute for Scientific
Computing and Information
INTERNATIONAL JOURNAL OF
NUMERICAL ANALYSIS AND MODELING
Volume 12, Number 1, Pages 125–143
NUMERICAL APPROXIMATION OF VISCOELASTIC FLOWS IN
AN ELASTIC MEDIUM
HYESUK LEE AND SHUHAN XU
(Communicated by Max Gunzburger )
Abstract. We study the numerical approximation of a flow problem governed by a viscoelastic model coupled with a one-dimensional elastic structure equation. A variational formulation
for flow equations is developed based on the Arbitrary Lagrangian-Eulerian (ALE) method and
stability of the time discretized system is investigated. For decoupled numerical algorithms we
consider both an explicit scheme of the Leap-Frog type and a fully implicit scheme.
Key words. viscoelastic flow, elastic structure, finite elements.
1. Introduction
There have been extensive numerical studies on multidisciplinary problems where
a coupled mathematical model is present for more than one media. In particular,
the interaction of fluid flows with an elastic medium is of great interest for both
industrial and biological uses. Examples of such systems include bronchial air ways
[10], blood flow in arteries [4, 21] and micro-fluidic devices [25]. A fluid-elastic
system is typically highly linked through a strong coupling of governing equations,
making its mathematical and numerical study very challenging. Subsystems interact with each other in the manner that the stress of the fluid imposes forces on
the boundaries of the elastic medium, which causes continuous displacement of the
flexible boundaries, resulting in a movement of the fluid domain. Such an interaction problem is mathematically described as a system of time-dependent, coupled,
linear/nonlinear partial differential equations with moving boundaries.
Modeling and simulation of the fluid-structure interaction have been studied by
many researchers [3, 5, 7, 27]. Fluid-structure interactions involving elementary
fluids, governed by equations for a potential functions, e.g., the Laplace equation
or wave equation, were studied in [3, 11]. Studies on multi-phase systems involving
viscous, incompressible fluids and elastic solids can be found in [5, 17, 21]. One
important application of such a system is a blood flow. Interactions of blood flow
with a vessel wall are modeled by coupled fluid-structure equations with matching
conditions on velocity and traction along the interface. The vessel wall is known
to be elastic material, thus a linear/nonlinear elastic equation has been used for
the structural mechanics. The blood flow has been modeled by the Stokes or the
Navier- Stokes equations in literature [4, 13, 21, 22], although blood is known to be
non-Newtonian in general. In particular, it is well known that blood flow in small
vessels and capillaries behaves as a viscoelastic fluid [19, 28].
There are only a few numerical studies found in literature for the non-Newtonian
fluids coupled with elastic solids. Some simulation results for the interaction of
non-Newtonian fluids with deformable bodies were reported in engineering journals
[1, 26] for the purpose of model validation. Even though there are some reports [15,
Received by the editors May 30, 2013 and, in revised form, January 13, 2014.
2000 Mathematics Subject Classification. 65M60, 65M12.
This research was supported by the NSF under grants DMS-1016182 and DMS-1418960.
125
126
H. LEE AND S. XU
Figure 1. Fluid domain
18] on numerical methods for simulation of blood flows using non-Newtonian fluid
models, detailed numerical analysis such as studies on stability of time-stepping
schemes is, in general, lacking from the current literature.
Mathematical study for partial differential equations governing a viscoelastic
fluid behavior is still far behind compared to advances in computing. It is wellknown that an analytical or numerical study of viscoelastic flows is very challenging
due to complexity of governing equations. One of the difficulties in simulating
viscoelastic flows arises from the hyperbolic nature of the constitutive equation for
which one needs to use a stabilization technique such as the streamline upwinding
Petrov-Galerkin (SUPG) [24] method or the discontinuous Galerkin method [2].
In [6, 16] we introduced a modified Johnson-Segalman model, referred to as the
Oseen-viscoelastic model, in which the velocity in the nonlinear pieces of the constitutive equation was taken to be a known function. For this model we have been
able to provide a rigorous mathematical analysis of the model equations and their
approximation, and provide some additional insights into the investigated problems
(high Weissenberg number problem, domain decomposition, optimal control) in viscoelasticity. We use the Oseen model for analysis in this paper, however numerical
tests will be for the standard Jonson-Segalman model. For the structure model, we
use a one-dimensional string model introduced in [20, 21].
The remainder of this paper is organized as follows. In the next section model
equations are described and stability of the coupled system is proved. The variational formulation of the problem based on the Arbitrary Lagrangian Eulerian is
introduced in Section 3. In Section 4 the time discretized variational formulation
and its stability are studied, and finally numerical results are presented in Section
5.
2. Model equations
We consider a viscoelastic flow problem, where flow equations are coupled with
a one-dimensional elastic structure model. Let Ω(t) be a bounded domain at time t
in IR 2 with the Lipschitz continuous boundary Γ(t). Suppose Γ(t) consists of three
parts: ΓD , ΓN and Γw (t), where ΓD ∪ ΓN is a fixed boundary and Γw (t) a moving
wall boundary.
Consider the viscoelastic model equations:
∂σ
(2.1)
+ u · ∇σ + gβ (σ, ∇u) − 2 α D(u) = 0 in Ω(t) ,
σ+λ
∂t
∂u
(2.2) Re
+ u · ∇u − ∇ · σ − 2(1 − α) ∇ · D(u) + ∇p = f in Ω(t) ,
∂t
(2.3)
∇·u
=
0 in Ω(t) ,
NUMERICAL APPROXIMATION OF VISCOELASTIC FLOWS IN AN ELASTIC MEDIUM 127
where σ denotes the extra stress tensor, u the velocity vector, p the pressure of
fluid, Re the Reynolds number, and λ is the Weissenberg number defined as the
product of the relaxation time and a characteristic strain rate. In (2.1) and (2.2),
D(u) := (∇u + ∇uT )/2 is the rate of the strain tensor, α a number such that
0 < α < 1 which may be considered as the fraction of viscoelastic viscosity, and f
the body force. In (2.1), gβ (σ, ∇u) is defined by
1+β
1−β
(σ∇u + ∇uT σ) −
(∇u σ + σ∇uT )
2
2
for β ∈ [−1, 1]. Note that (2.1) reduces to the Oldroyd-B model for the case β = 1.
For the viscoelastic fluid flow problem the major difficulty in establishing existence of a solution to the continuous variation formulation is the constitutive
equation (2.1). With this equation the nonlinear operator associated with the
model is neither coercive nor monotone. One way to overcome this difficulty is to
consider a nearby problem where the gβ term is linearized with the given velocity
b(x, t)(≈ u(x, t)). Consider the modified problem with the given velocity b in the
gβ term:
∂σ
σ+λ
(2.5)
+ u · ∇σ + gβ (σ, ∇b) − 2 α D(u) = 0 in Ω(t) ,
∂t
∂u
(2.6) Re
+ u · ∇u − ∇ · σ − 2(1 − α) ∇ · D(u) + ∇p = f in Ω(t) ,
∂t
(2.4)
gβ (σ, ∇u) :=
(2.7)
∇·u
=
0 in Ω(t) .
Initial and boundary conditions for u and σ are given as follows:
(2.8)
u(x, 0)
= u0
in Ω(0) ,
(2.9)
σ(x, 0) = σ 0
in Ω(0) ,
(2.10)
u
= 0
on ΓD ,
(2.11)
(σ + 2(1 − α) D(u) − pI) · n
= 0
on ΓN ,
where n is the outward unit normal vector to Ω(t).
We make the following assumption for b:
(2.12)
b ∈ H1 (Ω(t)),
∇ · b = 0,
kbk∞ ≤ M,
k∇bk∞ ≤ M < ∞ .
To analyze the equation (2.5), we will need a small data condition on the Weissenberg number λ and/or on the bound M so that 1 − 4λM > 0. As a model
equation of elastic structure we consider the one-dimensional generalized string
model [21], which was developed to account for longitudinal action:
(2.13)
ρw s
∂ 2η
∂3η
Es
∂2η
ˆ,
−
kGh
−
γ
+
η=Φ
2
2
2
∂t
∂z
∂t∂z
(1 − ν 2 )R02
where η represents the radial displacement of the wall with respect to the rest
configuration
(2.14)
Γ0 := {(z, r) ∈ IR 2 : z ∈ (0, L), r ∈ (0, R0 )}.
In (2.13) ρw is the wall volumetric mass, s the wall thickness, k the Timoshenko
shear correction factor, G the shear modulus, and E the Young modulus. The right
ˆ is the external force in the radial direction, ν the Poisson
hand side function Φ
128
H. LEE AND S. XU
ratio, γ a viscoelastic parameter and R0 is the radius at rest. In order to simplify
expressions, we rewrite (2.13) in the form of
(2.15)
ρw
∂2η
∂2η
∂3η
ˆ
−a 2 −b
+ c η = Φ,
2
∂t
∂t
∂t∂z 2
where a, b, c are positive constants related to the physical properties of the solid
structure described above. The structure equation is accompanied with the conditions at z = 0, L:
(2.16)
η|z=0 = 0 for all t,
η|z=L = 0
for all t.
The interface conditions between the fluid and the structure are obtained by
enforcing continuity of the velocity and the stress force:
(2.17)
∂η
e |t=0
∂t
=
u0
(2.18)
∂η
e
∂t
=
u
(2.19)
Φe =
on Γw (0) ,
on Γw (t) ,
−(−pI + σ + 2(1 − α)D(u)) · n − pext n
on Γw (t) ,
where e is a unit vector in the radial direction and pext is the external pressure.
Without loss of generality we let pext = 0 and Γ0 = Γw (0) in this paper. In
ˆ on Γw (t), which takes into account the change of
(2.19) Φ is a representation of Φ
ˆ will
configuration from the reference to the moving interface. A detailed form of Φ
be discussed in the next section. See (3.11).
Remark 2.1. In general the system of viscoelastic fluid equations (2.1)-(2.3) is
completed with initial and boundary conditions (2.8)-(2.11) and a Dirichlet type
boundary condition for the stress σ along an inflow boundary of fluid domain, i.e.,
on a part of Γ(t) where u · n < 0. In a fluid-structure system an inflow part on
the moving boundary is changed from time to time due to the interface condition
(2.18), which makes numerical studies for the system extremely challenging not
only by the change of inflow boundaries but also by a lack of boundary information
on the stress. Therefore, to simplify numerical analysis we assume that the stress
is unknown on the whole boundary. The analysis results in theorems are still valid
if a stress condition is imposed. A possible way to implement a stress boundary
condition is suggested in the conclusion section as a future work.
We use the Sobolev spaces W m,p (D) with norms k · km,p,D if p < ∞, k · km,∞,D
if p = ∞. We denote the Sobolev space W m,2 by H m with the norm k · km . The
corresponding space of vector-valued or tensor-valued functions is denoted by Hm .
If D = Ω(t), D is omitted, i.e., (·, ·) = (·, ·)Ω(t) and k · k = k · kΩ(t) . For σ, τ tensors
P2
and u, v vectors, we use : and · to denote the tensor product σ : τ := i,j=1 σij τij
P2
and the vector product u · v := i=1 ui vi . For the structure equation, we will use
(·, ·), k · k to denote (·, ·)Γ0 and k · kΓ0 , respectively.
In the next theorem we show the stability of a solution satisfying the coupled
problem (2.5)-(2.11), (2.15)-(2.19).
NUMERICAL APPROXIMATION OF VISCOELASTIC FLOWS IN AN ELASTIC MEDIUM 129
Theorem 2.2. If 1 − 4λM > 0 and u · n ≥ 0 on ΓN , a solution to the system
(2.5)-(2.11), (2.15)-(2.19) satisfies the estimate
λ
∂η
∂η
kσk20 + α Re kuk20 + ρw α k k20 + aα k k20 + cα kηk20
2
∂t
∂z
Z t
2
∂ η 2
k0 + (1 − 4λM )kσk20 + 2α(1 − α)kD(u)k20 ds
+
bαk
∂t∂z
0
(2.20)
≤ C,
where C is a constant depending on the forcing term f and initial data.
Proof. Multiplying (2.15) by ∂η
∂t and integrating over Γ0 , we have that
∂2η
∂2η
∂3η
∂η
∂η
ˆ
ρw 2 − a 2 − b
(2.21)
= Φ,
.
+ c η,
∂t
∂z
∂t∂z 2
∂t Γ0
∂t Γ0
Using integration by parts and (2.16), (2.21) implies
ρw d ∂η 2 a d ∂η 2 b ∂ 2 η 2 c d
k k +
k k + k
k +
kηk20 =
2 dt ∂t 0 2 dt ∂z 0 2 ∂t∂z 0 2 dt
(2.22)
ˆ ∂η
.
Φ,
∂t Γ0
On the other hand, multiplying (2.5), (2.6), (2.7) by σ, 2α u and p, respectively,
integrating over Ω(t) and using the Green’s theorem, we have
Z
λ
λ
d
(σ : σ) dΩ + ((u · n)σ, σ)Γw(t) ∪ΓN
2 Ω(t) dt
2
Z
d
(u · u) dΩ + α Re ((u · n)u, u)Γw(t) ∪ΓN + A((σ, u), (σ, u))
+α Re
Ω(t) dt
(2.23)
= 2α (f , u) + 2α((σ + 2(1 − α) D(u) − pI) · n, u)Γw(t) ,
where A((u, σ), (v, τ )) is defined by
A((u, σ), (v, τ )) :=
(2.24)
(σ, τ ) + λ (gβ (σ, ∇b), τ ) − 2α (D(u), τ )
+2α (σ, D(v)) + 4α(1 − α) (D(u), D(v)) .
Note that, since
(2.25)
(gβ (σ, ∇b), τ ) ≤ 4k∇bk∞ kσk0 kτ k0 ≤ 4M kσk0 kτ k0 ,
if λ M is small so that 1 − 4λM > 0, A satisfies
A((u, σ), (u, σ)) ≥
(2.26)
kσk20 − 4λM kσk20 + 4α(1 − α) kD(u)k20
= (1 − 4λM ) kσk20 + 4α(1 − α) kD(u)k20 .
Using (2.18) and the Reynolds transport formula
Z
Z
Z
∂ψ
d
(2.27)
dΩ =
ψ dΩ −
(w · n) ψ dΓ ,
dt Ω(t)
Ω(t) ∂t
Γ(t)
where w is the velocity of a moving boundary, (2.23) is reduced to
d λ
kσk20 + α Re kuk20
dt 2
λ
(2.28)
+ ((u · n)σ, σ)ΓN + α Re ((u · n)u, u)ΓN + A((σ, u), (σ, u))
2
= 2α (f , u) + 2α((σ + 2(1 − α) D(u) − pI) · n, u)Γw (t) .
130
H. LEE AND S. XU
Note that by the interface boundary conditions (2.18)-(2.19), the integral along the
interface boundary in (2.28) is written as
∂η
∂η
ˆ
.
e
= − Φ,
(2.29) ((σ +2(1−α) D(u)−pI)·n, u)Γw (t) = − Φ e,
∂t
∂t Γ0
Γw (t)
If u · n ≥ 0 on ΓN , two boundary integrals along ΓN in the left hand side of (2.28)
are positive, therefore, using (2.26), (2.28) and (2.29), we have
d λ
2
2
kσk0 + α Re kuk0 + (1 − 4λM ) kσk20 + 4α(1 − α) kD(u)k20
dt 2
∂η
ˆ
(2.30)
.
≤ 2α (f , u) − 2α Φ,
∂t Γ0
Now, multiplying (2.22) by 2α, adding to (2.30) and using the Young’s and Poincar´e
inequalities, we have that
d λ
∂η 2
∂η 2
2
2
2
kσk0 + α Re kuk0 + ρw α k k0 + aα k k0 + cα kηk0
dt 2
∂t
∂z
∂2η 2
+b α k
k + (1 − 4λM ) kσk20 + 4α(1 − α) kD(u)k20
∂t∂z 0 1
(2.31)
≤ 2α
kf k2−1 + (1 − α)kD(u)k20 .
4(1 − α)
Simplifying (2.31), integrating over the time from 0 to t and using the Gronwall
Lemma [23], the energy estimate (2.20) is obtained.
3. ALE formulation
For the fluid-elastic system, the interface moves by the displacement of structure, therefore, the fluid subproblem is a moving boundary problem. Numerical
simulation for the moving boundary problem can be performed using the Arbitrary
Lagrangian Eulerian method [14, 22], where the Eulerian frame is used in the fluid
domain while, in the solid domain, the Lagrangian formulation is used. In ALE
formulation an unknown coordinate transformation is usually introduced for the
fluid domain, and the fluid equations can be rewritten for a fixed reference domain. Specifically, we define the time-dependent bijective mapping Ψt that maps
the reference domain Ω0 to the physical domain Ω(t):
(3.1)
Ψt : Ω0 → Ω(t),
Ψt (y) = x(t, y) ,
where y and x are the spatial coordinates in Ω0 and Ω(t), respectively. The coordinate y is often called the ALE coordinate. Using Ψt the weak formulation of the
flow equations in Ω(t) can be recast into a weak problem defined in the reference
domain Ω0 . Thus the model equations in the reference domain can be considered for
numerical simulation and the transformation function Ψt needs to be determined
at each time step as a part of computation.
For a function φ : Ω(t) × [0, T ] → IR , define its corresponding function φ = g ◦ Ψt
in the ALE setting:
(3.2)
φ : Ω0 → IR ,
φ(y, t) = φ(Ψt (y), t).
The time derivative in ALE frame is also defined as
(3.3)
∂φ
|y : Ω(t) × [0, T ] → IR ,
∂t
∂φ
∂φ
|y (y, t) =
(y, t).
∂t
∂t
NUMERICAL APPROXIMATION OF VISCOELASTIC FLOWS IN AN ELASTIC MEDIUM 131
Using the chain rule,
∂φ
∂φ
|y =
+ z · ∇φ,
∂t
∂t
(3.4)
∂φ
where z := ∂x
∂t |y is the domain velocity. In (3.4) ∂t |y is the so called ALE derivative of φ. The flow equations (2.5)-(2.7) can be then written in ALE formulation
as follows:
∂σ
(3.5)
|y +(u − z) · ∇σ + gβ (σ, ∇b) − 2 α D(u) = 0 in Ω(t) ,
σ+λ
∂t
∂u
Re
|y +(u − z) · ∇u − ∇ · σ − 2(1 − α) ∇ · D(u)
∂t
(3.6)
+∇p = f in Ω(t) ,
(3.7)
∇ · u = 0 in Ω(t) .
In order to define the ALE mapping Ψt , consider the boundary position function
h : ∂Ω0 → ∂Ω(t) such that
y + η e on Γw (t)
(3.8)
h(y) =
y
on ΓD ∪ ΓN ,
where η is the displacement of the moving wall. The ALE mapping may be then
determined by solving the equation
(3.9)
∆Ψt
=
0
in Ω0 ,
Ψt
=
h
on ∂Ω0 .
This method is called the harmonic extension, where the boundary position function
on the boundary is extended onto the whole domain [8, 12].
The forcing term Φ in (2.19) given in the coordinate for Γw (t) can be recast in
the reference configuration Γ0 as

 s
2
∂η 
|y
on Γ0 ,
(3.10)
Φ(z, t)|x = Φ 1 +
∂z
q
2
where the expression 1 + ( ∂η
∂z ) represents the change in the surface measure
passing from Γw (t) to Γ0 . Therefore, (2.15) can be rewritten on Γ0 as
s
2
3
2
∂ η
∂ η
∂η 2
∂ η
+
c
η
=
Φ
.
1
+
(3.11)
ρw 2 − a 2 − b
∂t
∂z
∂t∂z 2
∂z
Now the fluid in a moving boundary problem in ALE formulation is summarized as
solve the fluid equations (3.5)-(3.7) and the structure equation (3.11)
with the boundary and initial conditions (2.8)-(2.11), (2.16) and
(3.12)
u
=
(3.13)
(−pI + σ + 2(1 − α)D(u)) · n
=
using the ALE mapping satisfying (3.9).
∂η
e on Γw (t) ,
∂t
Φ e on Γw (t) ,
132
H. LEE AND S. XU
For the variational formulation of the flow equations (3.5)-(3.7) in ALE framework, define function spaces for the reference domain:
U0 := {v ∈ H1 (Ω) : v = 0
on ΓD } ,
Q0 := L2 (Ω) ,
Σ0 := {τ ∈ L2 (Ω) : τ ij = τ ji } .
The function spaces for Ω(t) are then defined as
for v ∈ X0 } ,
Ut := {v : Ω(t) × [0, T ] → IR 2 , v = v ◦ Ψ−1
t
Qt := {q : Ω(t) × [0, T ] → IR , q = q ◦ Ψ−1
for p ∈ Q0 } ,
t
Σt := {τ : Ω(t) × [0, T ] → IR 2×2 , τ = τ ◦ Ψ−1
for τ ∈ Σ0 } .
t
The variational formulation of (3.5)-(3.7) in ALE framework is then to find (u, p, σ) ∈
Ut × Qt × Σt such that
∂σ
|y +(u − z) · ∇σ + gβ (σ, ∇b), τ
∂t
(3.14)
−2α (D(u), τ ) = 0 ∀τ ∈ Σt ,
∂u
Re
|y +(u − z) · ∇u, v + (σ, D(v)) + 2(1 − α)(D(u), D(v))
∂t
(3.15) −(p, ∇ · v) = (f , v) + ((σ + 2(1 − α) D(u) − pI) · n, v)Γw (t) ∀v ∈ Ut ,
(σ, τ ) + λ
(3.16) (q, ∇ · u) = 0
∀q ∈ Qt .
In order to derive the conservative variational formulation [20], consider the
Reynolds transport formula
(3.17)
Z
Z
Z
d
∂φ
∂φ
|y +φ∇x · z dV =
|x +∇x φ · z + φ∇x · z dV
φ(x, t) dV =
dt V (t)
∂t
V (t) ∂t
V (t)
for any subdomain V (t) ⊂ Ωt such that V (t) = Ψt (V0 ) with V0 ⊂ Ω0 . If v is a
for v : Ω0 → IR , we have that
function from Ωt to IR and v = v ◦ Ψ−1
t
∂v
|y = 0
∂t
(3.18)
and, therefore,
d
dt
(3.19)
(3.20)
d
dt
Z
Ωt
Z
v dΩ =
v∇x · z dΩ ,
Ωt
Ωt
φv dΩ =
Z
Z
Ωt
∂φ
|y +φ∇x · z v dΩ .
∂t
NUMERICAL APPROXIMATION OF VISCOELASTIC FLOWS IN AN ELASTIC MEDIUM 133
Using (3.14)-(3.16) and (3.20), we have the following weak formulation: find
(u, p, σ) ∈ Ut × Qt × Σt such that
(σ, τ ) + λ
(3.21)
Re
d
(u, v) + Re (−u(∇ · z) + (u − z) · ∇u, v) + (σ, D(v))
dt
+2(1 − α)(D(u), D(v)) − (p, ∇ · v)
(3.22)
(3.23)
d
(σ, τ ) + λ (−σ(∇ · z) + ((u − z) · ∇)σ + gβ (σ, ∇b), τ )
dt
−2α (D(u), τ ) = 0 ∀τ ∈ Σt ,
= (f , v) + ((σ + 2(1 − α) D(u) − pI) · n, v)Γw (t)
∀v ∈ Ut ,
(q, ∇ · u) = 0 ∀q ∈ Qt .
For the structure equation define the function space:
S := H01 (0, L).
The variational formulation of (3.11) is then to find η ∈ S such that
(3.24)
s
2
2
∂η 2
∂η ∂ξ
∂ η
∂ η ∂ξ
ρw ( 2 , ξ) + a ( ,
) +b(
,
) + c (η, ξ) = (Φ 1 +
, ξ)
∂t
∂z ∂z
∂t∂z ∂z
∂z
∀ξ ∈ S.
For the coupled problem define the test function spaces for u, η by
˜ t × S˜ := {(v, ξ) ∈ Ut × S : v ◦ Ψt |Γ0 = ξe} .
(3.25)
U
Using (2.19), (3.10) and the ALE mapping, we have
ρw (
(3.26)
∂η ∂ξ
∂ 2 η ∂ξ
∂2η
, ξ) + a ( ,
) +b(
,
) + c (η, ξ)
2
∂t
∂z ∂z
∂t∂z ∂z
= −((σ + 2(1 − α) D(u) − pI) · n, (ξ ◦ Ψ−1
t )e)Γw (t) .
Thus, by (2.24), (3.12) and (3.13), the variational problem of the fluid-structure
coupled problem in ALE frame is given by: find (u, p, σ, η) ∈ Ut × Qt × Σt × S
such that
∂η ∂ξ
∂ 2 η ∂ξ
∂2η
) +b(
,
) + c (η, ξ)
2α ρw ( 2 , ξ) + a ( ,
∂t
∂z ∂z
∂t∂z ∂z
d
+λ (σ, τ ) + λ (−σ(∇ · z), τ ) + λ ((u − z) · ∇)σ, τ )
dt
d
+2α Re (u, v) + 2α Re (−u(∇ · z), v) + 2α Re ((u − z) · ∇u, v)
dt
(3.27)
+A((σ, u), (τ , v)) − 2α(p, ∇ · v) + 2α(q, ∇ · v) = 2α (f , v)
˜ t × Σt × Qt × S,
˜ where A((σ, u), (τ , v)) is defined as (2.24).
∀(v, q, τ , ξ) ∈ U
4. Discretization
We define finite element spaces for the approximattion of (u, p) in Ω0 :
Uh,0
:=
{v ∈ U0 ∩ (C 0 (Ω))2 : v|K ∈ P2 (K)2 , ∀K ∈ Th,0 } ,
Qh,0
:=
{q ∈ Q0 ∩ C 0 (Ω) : q|K ∈ P1 (K), ∀K ∈ Th,0 } ,
where Th,0 is a triangularization of Ω0 . The stress σ is approximated in the discontinuous finite element space of piecewise linear polynomials:
Σh,0 := {τ ∈ Σ0 : τ |K ∈ P1 (K)2×2 , ∀K ∈ Th,0 } .
134
H. LEE AND S. XU
Then the finite element spaces for Ωt are defined as
Uh,t
:= {vh : Ωt × [0, T ] → IR 2 , vh = vh ◦ Ψ−1
h,t for vh ∈ Uh,0 } ,
Qh,t
:= {qh : Ωt × [0, T ] → IR , qh = q h ◦ Ψ−1
h,t for q h ∈ Qh,0 } ,
Σh,t
for σ h ∈ Σh,0 } ,
:= {σ h : Ωt × [0, T ] → IR 2×2 , σ h = σ h ◦ Ψ−1
t
where Ψh,t : Ω0 → Ωt is a discrete mapping approximated by Pl Lagrangian finite
elements such that Ψh,t (y) = xh (y, t). For the discrete ALE mapping, define the
set
Xh := {x ∈ H1 (Ω0 ) : x|K ∈ Pl (K)2 , ∀K ∈ Th,0 } .
(4.1)
For the displacement we define
Sh := {ξ ∈ S ∩ C 0 ([0, L]) : ξ|E ∈ P2 (E), ∀E ∈ T h } ,
(4.2)
where ∪T h = [0, L] and T h has matching grid points with Th,0 .
Introduce the operators θ(·, ·, ·) κ(·, ·, ·) defined by
θ(u, w, v)Ωt :=
(4.3)
1
[(u · ∇w, v)Ωt − (u · ∇v, w)Ωt ] ,
2
(4.4)
1
κ(v − z, σ, τ )Ωt := (((v − z) · ∇)σ, τ )Ωt + (∇ · v σ, τ )Ωt + < σ + − σ − , τ + >h,v−z ,
2
where the last term in (4.4) accounts for the jump in the discretized σ across an
inflow edge of element associated v − z. Using the Green’s theorem and ∇ · u = 0,
1
(u · ∇w, v)Ωt = θ(u, w, v)Ωt + ((u · n)w, v)ΓN ∪Γw (t)
2
(4.5)
and
θ(u, v, v)Ωt =
(4.6)
1
((u · n)v, v)ΓN ∪Γw (t)
2
∀v ∈ Uh,t .
Also, using integration by parts, we have
κ(uh − zh , σ h , τ h )Ωt
(4.7)
1
= −(((uh − zh ) · ∇)τ h , σ h )Ωt − (∇ · uh τ h , σ h )Ωt
2
−
−
+
+ < σ h , τ h − τ h >h,uh −zh +((∇ · zh )σ h , τ h )Ωt
+(((uh − zh ) · n)σ h , τ h )ΓN ∪Γw (t) ,
therefore,
1
[((∇ · zh )σ h , σ h )Ωt
2
+(((uh − zh ) · n)σ h , σ h )ΓN ∪Γw (t) + < σ h + − σ h − , σh + − σ h − >h,uh −zh
1
(4.8)
((∇ · zh )σ h , σ h )Ωt + (((uh − zh ) · n)σ h , σ h )ΓN ∪Γw (t) .
≥
2
κ(uh − zh , σ h , σ h )Ωt =
NUMERICAL APPROXIMATION OF VISCOELASTIC FLOWS IN AN ELASTIC MEDIUM 135
The semi-discrete variational formulation of the coupled problem in ALE frame
is then written as
∂ 2 ηh
∂ηh ∂ξh
∂ 2 ηh ∂ξh
2α ρw ( 2 , ξh )Γ0 + a (
,
)Γ0 + b (
,
)Γ0 + c (ηh , ξh )Γ0
∂t
∂z ∂z
∂t∂z ∂z
d
(σ h , τ h )Ωt + κ(uh − zh , σ h , τ h )Ωt − (σ h (∇ · zh ), τ h )Ωt
+λ
dt
d
+2α Re
(uh , vh )Ωt − (uh (∇ · zh ), vh )Ωt
dt
1
+ θ(uh , uh , vh )Ωt + ((uh · n)uh , vh )ΓN ∪Γw (t) − (zh · ∇uh , vh )Ωt
2
+A((σ h , uh ), (τ h , vh )) − 2α(ph , ∇ · vh )Ωt + 2α(qh , ∇ · uh )Ωt
˜ h × Σh × Qh × S˜h .
= 2α (f , v)Ωt ∀(vh , τ h , qh , ξh ) ∈ U
(4.9)
˜ h × S˜h is a subspace of Uh × Sh , where the matching condition in (3.25)
In (4.9) U
is satisfied in a discrete sense.
Using the backward Euler for both the fluid and structure equations and applying
the geometric conservation law(GCL) [20], we have the fully discretized systems
given as follows.
∂ηhn+1 ∂ξh
η n+1 − 2ηhn + η n−1
+
a
(
,
ξ
)
,
)Γ
2α ρw ( h
h Γ0
∆t2
∂z
∂z 0
+b (
n+1
∂ηh
∂z
−
∆t
n
∂ηh
∂z

∂ξh
,
)Γ + c (ηhn+1 , ξh )Γ0 
∂z 0
i
λ h n+1
(σ h , τ h )Ωtn+1 − (σ nh , τ h )Ωtn
∆t
2α Re n+1
(uh , vh )Ωtn+1 − (unh , vh )Ωtn
+
h∆t
+
i
+λ κ(un+1
− zn+1
, σ n+1
, τ h )Ω n+ 1 − (σ n+1
(∇ · zn+1
), τ h )Ω n+ 1
h
h
h
h
h
2
2
t
t
1
· n)un+1
, vh )ΓN ∪Γ n+ 1
+2α Re θ(un+1
, un+1
, vh )Ω n+ 1 + ((un+1
h
h
h
h
2
2
2
t
t
i
n+1
n+1
n+1
n+1
−(uh (∇ · zh ), vh )Ω n+ 1 − (zh · ∇uh , vh )Ω n+ 1
2
t
t
2
+A((σ n+1
, un+1
), (τ h , vh ))Ω n+ 1 − 2α(pn+1
, ∇ · vh )Ω n+ 1 + 2α(qh , ∇ · uh )
h
h
h
2
t
(4.10)
t
2
= 2α (f , vh )Ω n+ 1 .
t
2
Stability of the fully discretized coupled system is proved in the next theorem.
We omit the subscript h in (unh , σ nh , pnh , ηhn ) to simplify our notation.
136
H. LEE AND S. XU
Theorem 4.1. If 1 − 4λM > 0, a solution to the fully discretized system (4.10)
satisfies the estimate
α Re n+1 2
λ
kσ n+1 k20,Ωtn+1 +
ku
k0,Ωtn+1
αc kη n+1 k20 +
2∆t
∆t
n+1
∂η
ρw n+1
kη
− η n k20 + a k
k2
+α
2
∆t
∂z 0
n
X
ρw i+1
∂η i+1 − ∂η i 2
α
+
kη
− 2η i + η i−1 k20 + a k
k0
2
∆t
∂z
i=0
2b ∂(η i+1 − η i ) 2
k0 + c kη i+1 − η i k20
+ k
∆t
∂z
n
X
2α(1 − α)kD(ui+1 )k20,Ω 1 + (1 − 4λM )kσ i+1 k20,Ω
+
i+ 1
2
t
i+
2
t
λ
((ui+1 · n)σ i+1 , σ i+1 )ΓN + α Re ((ui+1 · n)ui+1 , ui+1 )ΓN
2
i=0
λ
α Re
∂η0 2
ρw
2
2
≤α
kη˙ 0 k0 + a k
k0 + c kη0 k0 +
kσ 0 k20,Ω0 +
ku0 k20,Ω0
∆t
∂z
2∆t
∆t
n
X
kf i+1 k2−1,Ω 1 .
+C
+
(4.11)
i=0
n X
∆t
i+
2
t
i=0
Proof. Letting τ = σ n+1 , v = un+1 , q = pn+1 and ξ = η n+1 − η n in (4.10), we
obtain
η n+1 − 2η n + η n−1 n+1
∂η n+1 ∂η n+1 − ∂η n
n
+
a
(
,
η
−
η
)
,
)Γ0
Γ
0
∆t2
∂z
∂z
#
n
∂η n+1
− ∂η
∂η n+1 − ∂η n
n+1 n+1
n
∂z
∂z
,η
− η )Γ0
,
)Γ0 + c (η
+b (
∆t
∂z
2α ρw (
λ
2α Re n+1 2
kσ n+1 k20,Ωtn+1 +
ku
k0,Ωtn+1
∆t
∆t
h
i
+λ κ(un+1 − zn+1 , σ n+1 , σ n+1 )Ω n+ 1 − (σ n+1 (∇ · zn+1 ), σ n+1 )Ω n+ 1
2
2
t
t
1
n+1
n+1
n+1
n+1
n+1
n+1
· n)u
,u
)ΓN ∪Γ n+ 1
+2α Re θ(u
,u
,u
)Ω n+ 1 + ((u
2
2
2
t
t
i
−(un+1 (∇ · zn+1 ), un+1 )Ω n+ 1 − (zn+1 · ∇un+1 , un+1 )Ω n+ 1
+
2
t
+A((σ
n+1
,u
n+1
), (σ
n+1
,u
n+1
t
2
))Ω n+ 1
t
2
2α Re n n+1
λ
(4.12) =
(σ n , σ n+1 )Ωtn +
(u , u
)Ωtn + 2α (f n+1 , un+1 )Ω n+ 1 .
∆t
∆t
2
t
NUMERICAL APPROXIMATION OF VISCOELASTIC FLOWS IN AN ELASTIC MEDIUM 137
Consider the left hand side of the equation first. The structure terms are turned
to:
∂η n+1 ∂η n+1 − ∂η n
η n+1 − 2η n + η n−1 n+1
n
+
a
(
,
η
−
η
)
,
)Γ0
2α ρw (
Γ
0
∆t2
∂z
∂z
#
n
∂η n+1
− ∂η
∂η n+1 − ∂η n
∂z
,
)Γ0 + c (η n+1 , η n+1 − η n )Γ0
+b ( ∂z
∆t
∂z
hρ
w
=α
kη n+1 − η n k20 − kη n − η n−1 k20 + kη n+1 − 2η n + η n−1 k20
2
∆t ∂η n 2
∂η n+1 − ∂η n 2
2b ∂(η n+1 − η n ) 2
∂η n+1 2
k0 − k
k0 + k
k0 +
k
k0
+a k
∂z
∂z
∂z
∆t
∂z
(4.13)
+c kη n+1 k20 − kη n k20 + kη n+1 − η n k20 .
By (4.8),
h
i
λ κ(un+1 − zn+1 , σ n+1 , σ n+1 )Ω n+ 1 − (σ n+1 (∇ · zn+1 ), σ n+1 )Ω n+ 1
2
2
t
t
λ h n+1
(∇ · zn+1 ), σ n+1 )Ω n+ 1
≥ − (σ
2
2
t
i
n+1
n+1
n+1
(4.14)
−(((u
−z
) · n)σ
, σ n+1 )ΓN ∪Γ n+ 1 .
t
n+1
Using integration by parts and that z
(zn+1 · ∇un+1 , un+1 )Ω n+ 1
t
2
(4.15)
2
= 0 on the fixed boundary ΓN ,
1
= − ((∇ · zn+1 )un+1 , un+1 )Ω n+ 1
2
2
t
1 n+1
n+1
n+1
+ ((z
· n)u
,u
)ΓN ∪Γ n+ 1 .
2
2
t
Thus, using (4.6) and (4.15),
1
2α Re θ(un+1 , un+1 , un+1 )Ω n+ 1 + ((un+1 · n)un+1 , un+1 )ΓN ∪Γ n+ 1
2
2
2
t
t
i
n+1
n+1
n+1
n+1
n+1
n+1
· ∇u
,u
)Ω n+ 1
−(u
(∇ · z
), u
)Ω n+ 1 − (z
2
2
t
t
h
n+1
n+1
n+1
≥ α Re −(u
(∇ · z
), u
)Ω n+ 1
2
t
i
n+1
n+1
(4.16)
+(((u
−z
) · n)un+1 , un+1 )ΓN ∪Γ n+ 1 ,
2
t
and, by (2.26), (4.14) and (4.16),
Fluid termsh at left
≥ α Re −(un+1 (∇ · zn+1 ), un+1 )Ω n+ 1
t
+(((u
n+1
n+1
−z
) · n)u
2
n+1
, un+1 )ΓN ∪Γ n+ 1
t
2
+(1 − 4λM )kσn+1 k20,Ω 1 + 4α(1 − α)kD(un+1 )k20,Ω 1
n+
n+
2
2
t
t
λ h n+1
(∇ · zn+1 ), σ n+1 )Ω n+ 1
− (σ
2
2
t
i
−(((un+1 − zn+1 ) · n)σ n+1 , σ n+1 )ΓN ∪Γ n+ 1
t
(4.17)
i
λ
2α Re n+1 2
+ kσ n+1 k20,Ωtn+1 +
ku
k0,Ωtn+1 .
∆t
∆t
2
138
H. LEE AND S. XU
On the other hand, using the Schwartz and Young’s inequalities, the right hand
side of (4.12) is bounded as
RHS
(4.18)
≤
α Re
λ
kσn k20,Ωtn + kσn+1 k20,Ωtn +
kun k20,Ωtn + kun+1 k20,Ωtn
2∆t
∆t
+Ckf n+1 k2−1,Ω 1 + δkD(un+1 )k20,Ω 1 .
n+
2
t
n+
2
t
By (4.13), (4.17) and (4.18), (4.12) imply
hρ
w
kη n+1 − η n k20 + kη n+1 − 2η n + η n−1 k20
α
2
∆t
∂η n+1 − ∂η n 2
∂η n+1 2
k0 + k
k0
+a k
∂z
∂z
n+1
n
2b
∂(η
−η ) 2
+
k
k0 + c kη n+1 k20 + kη n+1 − η n k20
∆t
∂z
λ
1
+
kσ n+1 k20,Ωtn+1 − kσ n+1 k20,Ωtn
∆t
2
1
2α Re
kun+1 k20,Ωtn+1 − kun+1 k20,Ωtn
+
∆t
2
+(1 − 4λM )kσn+1 k20,Ω
+ (4α(1 − α) − δ)kD(un+1 )k20,Ω
n+ 1
2
t
n+ 1
2
t
λh
− (σ n+1 (∇ · zn+1 ), σ n+1 )Ω n+ 1
2
2
t
i
−(((un+1 − zn+1 ) · n)σ n+1 , σ n+1 )ΓN ∪Γ n+ 1
2
t
h
+α Re −(un+1 (∇ · zn+1 ), un+1 )Ω n+ 1
t
2
i
) · n)un+1 , un+1 )ΓN ∪Γ n+ 1
2
t
n
∂η
λ
ρw n
kη − η n−1 k20 + ak
k2 + ckη n k20 +
kσ n k20,Ωtn
≤α
∆t2
∂z 0
2∆t
α Re
(4.19)
+
kD(un )k20,Ωtn + C kf n+1 k2−1,Ω 1 .
n+
∆t
2
t
+(((u
n+1
n+1
−z
The time discretization scheme in (4.10) is based on the mid-point rule satisfying
GCL [20]
(4.20)
Z tn+1 Z
Z
Z
Z
vh dΩ =
vh dΩ −
vh ∇ · zh dΩ ,
vh ∇ · zh dΩ dt = ∆t
Ωtn+1
Ωtn
tn
and we have
i
λ h n+1 2
kσ
k0,Ωtn+1 − kσn+1 k20,Ωtn
2∆t
Ω
Ωt
=
=
n+ 1
2
t
λ
2
Z
kσ n+1 k20,Ω
Ω
n+ 1
2
t
n+ 1
2
t
∇ · zn+1 dΩ
λ n+1
(σ
(∇ · zn+1 ), σ n+1 )Ω n+ 1 ,
2
2
t
Z
i
α Re h n+1 2
ku
k0,Ωtn+1 − kun+1 k20,Ωtn
∆t
= α Re
(4.21)
= α Re (un+1 (∇ · zn+1 ), un+1 )Ω n+ 1 .
Ω
n+ 1
2
t
kun+1 k20,Ω
n+ 1
2
t
∇ · zn+1 dΩ
t
2
NUMERICAL APPROXIMATION OF VISCOELASTIC FLOWS IN AN ELASTIC MEDIUM 139
Using (4.21) in (4.19) and letting δ = 2α(1 − α), we obtain
hρ
w
kη n+1 − η n k20 + kη n+1 − 2η n + η n−1 k20
α
2
∆t
∂η n+1 2
∂η n+1 − ∂η n 2
+a k
k +k
k0
∂z 0
∂z
∂(η n+1 − η n ) 2
2b
n+1 2
n+1
n 2
k0 + kη
− η k0
k
k0 + c kη
+
∆t
∂z
λ
α Re n+1 2
+
kσ n+1 k20,Ωtn+1 +
ku
k0,Ωtn+1
2∆t
∆t
n+1 2
+(1 − 4λM )kσ
k0,Ω 1 + (2α(1 − α) − δ)kD(un+1 )k20,Ω 1
n+
2
t
n+
2
t
λ
+ (((un+1 − zn+1 ) · n)σ n+1 , σ n+1 )ΓN ∪Γ n+ 1
2
2
t
+α Re (((un+1 − zn+1 ) · n)un+1 , un+1 )ΓN ∪Γ n+ 1
t
2
∂η n 2
λ
ρw n
n−1 2
n 2
kη − η
k0 + ak
k + ckη k0 +
kσ n k20,Ωtn
≤α
∆t2
∂z 0
2∆t
α Re
(4.22)
+
kD(un )k20,Ωtn + C kf n+1 k2−1,Ω 1 .
n+
∆t
2
t
Summing over n in (4.22), assuming that the fluid velocity matches with the domain
velocity on Γ n+ 21 and using zn+1 = 0 on ΓN , we obtain the estimate (4.11).
t
5. Numerical results
In this section we present results of numerical experiments for the viscoelastic
fluid-structure system (2.1)-(2.3) and (2.15). The initial domain of the fluid is
the rectangle of height H = 1 and length L = 6 whose upper bound is elastic. See
Figure 2. Both the fluid and structure are initially at rest. We simulate the pressure
pulse Pin = 2000 by imposing the following Neumann boundary conditions on the
inflow and outflow sections:
πt
) − 1]n on Γin
(σ + 2(1 − α)D(u) − pI) · n = − P2in [cos( 0.0025
(σ + 2(1 − α)D(u) − pI) · n = 0
on Γout .
The parameters for the fluid equations are given as α = 0.9825, β = 0, λ = 0.9 ,
1/Re = 0.035, while for the structure equation, a = 25000, c = 0.01, b = 400000,
ρw = 1.1. Note that for higher Weissenberg number λ, the numerical stability
for the coupled system degrades due to the required condition 1 − 4λM > 0 of
Theorem 4.1. One of the main goals of numerical experiments is to compare the
viscoelastic case with the Newtonian case (λ = 0), and for this purpose, a choice of
low Reynolds number (high viscosity) can be made to improve numerical stability
when a larger λ value is used for simulations.
A conforming space discretization is applied to the fluid-structure coupled system. For the fluid, the space discretization consists of the Taylor-Hood (P2 , P1 )
finite elements for u, p and P1 discontinuous elements for σ. The structure is discretized by continuous P2 finite elements, and the ALE mapping is approximated
by P1 elements.
The homogeneous Dirichlet boundary condition was applied for the structure
equation in the previous analysis. However, in the numerical tests, we consider the
first order absorbing boundary condition instead [7], which may be more realistic
140
H. LEE AND S. XU
Figure 2. Test domain
when applied to blood flow simulations:

q
a ∂η
 ∂η −
=0
∂t
q ρw ∂z
a ∂η
 ∂η +
∂t
ρw ∂z = 0
at z = 0
at z = L .
5.1. Experiment 1. In this experiment, we approximate the system by the LeapFrog algorithm. We set the mesh size to h = 0.1 and the time step to ∆t = 10−4 .
The wall displacements of the viscoelastic fluid in different time steps are presented
in Figure 3. Observe that the bump of the wall, which is caused by fluid stress,
moves from the inflow section to the outflow section repeatedly as time goes. The
results meet our expectation based on the physical foundation. However, during
numerical tests, we found out that the explicit Leap-Frog algorithm is not always
stable. In fact, with a higher pressure input, the explicit algorithm can converge
only up to the time 0.14s.
5.2. Experiment 2. We implemented a more stable algorithm for the system
by an implicit scheme with relaxation. Specifically, a sub-iteration involving both
structure and fluid solvers is added in each time iteration [20]:
(1) find the fluid subproblem solutions un+1
, pn+1
, σ n+1
,
k
k
k
n+1
(2) find the structure subproblem solutions ηˆk ,
n+1
(3) relax the structure solver by ηkn+1 = ω ηˆkn+1 + (1 − ω)ηk−1
, 0 ≤ ω ≤ 1,
n+1
and keep iterating until |ηkn+1 − ηk−1
| < tolerance.
For this implicit algorithm, we used ω = 0.9 and pin = 20000, which is 10 times
higher than the previous experiment and the result of this case is presented in
Figure 4. Notice that the viscoelastic fluid is reduced to the Newtonian fluid in
the case of λ = 0. Since the Newtonian flow is usually used to simulate the blood
flow in many tests, we compared the wall displacements for both the viscoelastic
and the Newtonian models with all the same parameters except λ. In Figure 4,
displacements of the viscoelastic model are represented by blue curves while results
for the Newtonian case are plotted with red curved. Clearly, similar patterns are
observed from both models. A visible difference is observed as time goes, although
the difference is not quite significant at initial times. During numerical tests, we
also noticed that the viscoelastic case calls for more subiterations to converge as
expected.
6. Conclusion
Viscoelastic flow model equations coupled with the String model was considered
for stability analysis and numerical experiments. The numerical results indicate the
approach based on the ALE method can be used to simulate the viscoelastic flow in
an elastic medium. Non-negligible differences between the viscoelastic flows and the
Newtonian flows were observed under a high pressure input. Also we noticed that
NUMERICAL APPROXIMATION OF VISCOELASTIC FLOWS IN AN ELASTIC MEDIUM 141
−3
−3
5
x 10
5
2.5
2.5
0
0
−2.5
−2.5
−5
0
1
2
3
4
5
6
x 10
−5
0
1
(a) 0.02 s
5
0
0
−2.5
−5
0
1
2
3
4
5
6
6
4
5
6
4
5
6
4
5
6
4
5
6
x 10
−5
0
1
(c) 0.06 s
2
3
(d) 0.08 s
−3
−3
x 10
5
x 10
2.5
2.5
0
0
−2.5
−2.5
−5
0
1
2
3
4
5
6
−5
0
1
(e) 0.10 s
5
2.5
2.5
0
0
−2.5
−2.5
1
2
3
4
5
6
x 10
−5
0
1
(g) 0.14 s
2
3
(h) 0.16 s
−3
−3
x 10
5
2.5
2.5
0
0
−2.5
−2.5
−5
0
3
(f) 0.12 s
x 10
−5
0
2
−3
−3
5
5
2.5
−2.5
5
4
(b) 0.04 s
x 10
2.5
5
3
−3
−3
5
2
1
2
3
(i) 0.18 s
4
5
6
x 10
−5
0
1
2
3
(j) 0.20 s
Figure 3. Displacement of the wall for Pin = 2000
the non-Newtonian property affects stability of decoupled algorithms significantly,
i.e., larger Weissenberg numbers resulted in divergence of the algorithm at earlier
times. This is obviously due to the small data assumption on λ for numerical
stability of the coupled system.
In the future work the study on the viscoelastic/non-Newtonian fluid-structure
system will be extended to a 2D-2D case with linear/nonlinear structure models.
Various decoupling time discretization schemes will be investigated for better stabilized algorithms. As pointed out earlier more realistic numerical simulations include
142
H. LEE AND S. XU
0.04
0.04
0.03
0.03
0.02
0.02
0.01
0.01
0
0
−0.01
−0.01
−0.02
−0.02
−0.03
−0.03
−0.04
0
1
2
3
4
5
6
−0.04
0
1
(a) 0.02 s
0.04
0.03
0.03
0.02
0.02
0.01
0
−0.01
−0.01
−0.02
−0.02
−0.03
5
6
4
5
6
4
5
6
4
5
6
−0.03
1
2
3
4
5
6
−0.04
0
1
(c) 0.06 s
2
3
(d) 0.08 s
0.04
0.04
0.03
0.03
0.02
0.02
0.01
0.01
0
0
−0.01
−0.01
−0.02
−0.02
−0.03
−0.03
1
2
3
4
5
6
−0.04
0
1
(e) 0.10 s
2
3
(f) 0.12 s
0.04
0.04
0.03
0.03
0.02
0.02
0.01
0.01
0
0
−0.01
−0.01
−0.02
−0.02
−0.03
−0.04
0
4
0.01
0
−0.04
0
3
(b) 0.04 s
0.04
−0.04
0
2
−0.03
1
2
3
(g) 0.14 s
4
5
6
−0.04
0
1
2
3
(h) 0.16 s
Figure 4. Displacement of the wall for Pin = 20000
implementation of a stress boundary condition along inflow boundaries. One possible way would be to compute boundary value of stress using velocity information as
given in [9] on inflow boundaries, where u · n < 0 and use this as a stress condition
for the next time or subiteration.
References
[1] J. Adepua and V. Shankar, Suppression or enhancement of interfacial instability in two-layer
plane couette flow of fene-p fluids past a deformable solid layer, J. Non-Newtonian Fluid
Mech., 141 (2007) 43-58.
[2] J. Baranger and D. Sandri, Finite element approximation of viscoelastic fluid flow: Existence
of approximate solutions and error bounds. I. Discontinuous constraints, Numer. Math., 63
(1992) 13-27.
[3] J. Boujot. Mathematical formulation of fluid-structure interaction problems, RAIRO Model.
Math. Anal. Numer., 21 (1987) 239-260.
[4] S. Canic, A. Mikelic, and J. Tambaca, Two dimensional effective model describing fluidstructure inter- action in blood flow: analysis, simulation and experimental validation,
Comptes Rendus Mechanique Acad. Sci., 333 (2005) 867-883.
[5] P. Causin, J.F. Gerbeau, and F. Nobile, Added-mass effect in the design of partitioned
algorithms for fluid-structure problems, Comput. Methods. Appl. Mech. Engrg., 194 (2005)
4506-4527.
NUMERICAL APPROXIMATION OF VISCOELASTIC FLOWS IN AN ELASTIC MEDIUM 143
[6] V.J. Ervin, J.S. Howell, and H. Lee, A two-parameter defect-correction method for computation of steady-state viscoelastic fluid flow, Appl. Math. and Compt., 196 (2008) 818-834.
[7] L. Formaggia, J.F. Gerbeau, F. Nobile, and A. Quarteroni, On the coupling of 3D and 1D
Navier-Stokes equations for flow problems in compliant vessels, Comput. Methods Appl. Mech.
Engrg. 191 (2001) 561-582.
[8] L. Formaggia and F. Nobile, A stability analysis for the arbitrary Lagrangian Eulerian formulation with finite elements, East-West J. Numer. Math., 7 (1999) 105-131.
[9] A. Fortin and A. Zine, An improved GMRES method for solving viscoelastic fluid flow problem, J. of Non-Newtonian Fluid Mechanics, 42 (1992) 1-18.
[10] Y.C. Fung, Biomechanics Circulations, Splinger-Verlag, New York, 1996.
[11] L. Gastaldi, Mixed finite element methods in fluid structure system, Numer. Math., 74 (1996)
153-176.
[12] L. Gastaldi, A priori error estimates for the Arbitrary Lagrangian Eulerian formulation with
finite elements, East-West J. Numer. Math., 9 (2001) 123-156.
[13] J.J. Heys, T.A. Manteuffel, S.F. McCormick, and J.W. Ruge, First-order system least squares
(FOSLS) for coupled fluid-elastic problems, Journal of Computational Physics, 195 (2004)
560-575.
[14] T.J.R. Hughes, W.K. Liu and T.K. Zimmermann, Lagrangian-Eulerian finite element formulation for incompressible viscous flows, Comput. Methods. Appl. Mech. Engrg., 29 (1981)
329-349.
[15] J. Janela, A. Moura, and A. Sequeira, A 3D non-Newtonian uidstructure interaction model
for blood ow in arteries, J. Comput. Appl. Math., 234 (2010) 2783-2791.
[16] H. Lee and H.C. Lee, Analysis and finite element approximation of an optimal control problem
for the Oseen viscoelastic fluid flow, J. Math.Anal. Appl., 336 (2007) 1090-1160.
[17] P. LeTallec and S. Mani, Numerical analysis of a linearized fluid-structure interaction problem, Numer. Math., 87 (2000) 317-354.
[18] M. Lukacova and A. Zauskoza, Numerical modeling of shear-thinning non-Newtonian flows
in compliant vessels, Int. J. Numer. Meth. Fluids, 56 (2008) 1409-1415.
[19] D.A. McDonald. Blood Flow in Arteries. Williams & Wilkins, second edition, 1974.
[20] F. Nobile, Numerical approximation of fluid-structure interaction problems with application
to haemodynamics, Ph.D. thesis, 2001.
[21] A. Quarteroni, M. Tuveri and A. Veneziani, Computational vascular fluid dynamics: problems, models, and methods, Comput. Visual Sci., 2 (2000) 163-197.
[22] A. Quarteroni and L. Formaggia, Mathematical modelling and numerical simulation of the
cardiovascular system, Modelling of living systems, Handbook of numerical analysis series,
North-Holland, Amsterdam, 2004.
[23] A. Quarteroni and A. Valli, Numerical approximation of partial differential equations,
Springer Series in Computational Mathematics No. 23., Springer, 1994.
[24] D. Sandri, Finite element approximation of viscoelastic fluid flow: Existence of approximate
solutions and error bounds. Continuous approximation of the stress, SIAM J. Numer. Anal.
31 (1994) 362-377.
[25] K.P. Selverov and H.A. Stone. Peristaltically driven channel flows with applications toward
micro- mixing. Physics of Fluids, 13 (2000) 1837-1859.
[26] V. Shankar and S.Kumer, Instability of viscoelastic plane couette flow past a deformable wall,
J. Non-Newtonian Fluid Mech., 116 (2004) 371-393.
[27] A. Veneziani and C. Vergara, Flow rate defective boundary conditions in haemodynamics
simulations. Internat, J. Numer. Methods Fluids, 47 (2005) 803-816.
[28] X.Y. Xu, M.W. Collins, and C.J.H. Jones, Flow studies in canine aortas, ASME J. of
Biomech. Eng., 114 (1992) 504-511.
Department of Mathematical Sciences, Clemson University, Clemson, SC 29634-0975, USA
E-mail : [email protected] and [email protected]