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]
© Copyright 2024