- Email: [email protected]
IB-WENO method for incompressible flow with elastic boundaries
Accepted Manuscript IB-WENO method for incompressible flow with elastic boundaries Ziqiang Cheng, Yuan Liu, Mengping Zhang, Jin Wang
PII: DOI: Reference:
S0377-0427(18)30635-6 https://doi.org/10.1016/j.cam.2018.10.028 CAM 11973
To appear in:
Journal of Computational and Applied Mathematics
Received date : 12 January 2018 Revised date : 5 July 2018 Please cite this article as: Z. Cheng, et al., IB-WENO method for incompressible flow with elastic boundaries, Journal of Computational and Applied Mathematics (2018), https://doi.org/10.1016/j.cam.2018.10.028 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
IB-WENO method for incompressible flow with elastic boundaries Ziqiang Cheng a
Yuan Liu b a School
Mengping Zhang a
Jin Wang c,∗
of Mathematical Sciences,
University of Science and Technology of China, Hefei, Anhui 230026, China b Department
of Mathematics and Statistics,
Mississippi State University, Mississippi State, MS 39762, USA c Department
of Mathematics, University of Tennessee at Chattanooga, Chattanooga, TN 37403, USA ∗ Email:
[email protected]
Abstract We present a new approach to compute the interaction between an elastic boundary and an incompressible viscous fluid, by integrating the immersed boundary method and a finite difference WENO scheme. The hybrid algorithm maintains simple and efficient treatment of the moving interface based on the immersed boundary formulation, and achieves high-order accuracy, with the WENO method and the Hodge decomposition technique employed, in the majority of the fluid domain that does not contain the interface. We demonstrate the performance of the proposed algorithm and its various numerical components through a number of computational examples. Key Words: Immersed boundary formulation; Hodge decomposition; WENO method; Incompressible flow.
1
Introduction
The study of fluid-structure interaction (FSI), an emerging field concerned with the interaction between viscous fluid flows and immersed solid structures, has been fast growing in recent decades. FSI often involves strongly nonlinear, multi-physics phenomena, with application to a wide range of scientific and engineering disciplines [4, 15]. Owing to the difficulties in mathematical analysis and limitations in physical experiments for these complex problems, FSI research and development largely rely on scientific computing. 1
There have been a large number of numerical methods developed for the computation of FSI problems. These methods can be broadly classified into the monolithic approach [17, 29], where a unified system represents the entire FSI problem and a single algorithm solves the fluid and structural dynamics simultaneously, and the partitioned approach [34, 35] that treats the fluid and the structure as two systems which can be computed separately in their respective domains and then connected at the interface. Compared to the monolithic counterpart, a partitioned scheme typically has a higher degree of flexibility in algorithm development and implementation [15]. Another common classification of FSI methods is based on the types of meshes employed. Conforming-mesh algorithms, which include Lagrangian methods, arbitrary Lagrangian-Eulerian methods, reduced-order modeling techniques, etc., use meshes that fit the body of the moving and deformable structure. Consequently, a meshupdating procedure is normally required at each time step for these methods. In contrast, algorithms based on non-conforming meshes employ Cartesian background meshes that remain fixed throughout the course of computation, offering advantages such as improved efficiency, ease of implementation, and natural handling of large structural deformation. Among these numerical techniques, the immersed type methods, originated from the immersed boundary (IB) method [27] invented by Peskin in 1970s, have achieved much success and popularity. These methods belong to the non-conforming mesh category, and are typically based on the partitioned approach. An immersed method normally represents the effects of the immersed solid structure by adding a forcing term to the right-hand side of the fluid momentum equation. As a result, the fluid motion can be computed on fixed grids in the entire computational domain, eliminating the need of costly re-meshing procedures. Since Peskin’s pioneering work, many variations and improvements have been made to the family of immersed methods; significant examples include the immersed interface method [20, 21], the immersed finite element method [22,37], and the direct-forcing method [6,25]. Still, the original IB method remains a very popular option in FSI computation, due to its simple formulation, straightforward implementation, and applicability to a wide range of FSI problems. Nevertheless, the IB method has several drawbacks. To succeed, it has to utilize a discrete delta function and represent the sharp interface by a band of cells with non-zero thickness. Thus the method is inherently first-order accurate in space. Moreover, it has been observed that the immersed boundary method may suffer from leakage problem [20, 28], although mass conservation can be improved via the use 2
of divergence-free finite-difference operators [28]. Common practice of using the IB method in FSI computation is to apply the Fast Fourier Transform (FFT) in the periodic space dimension(s), and utilize standard finite difference approximations in those space dimensions that are not periodic. In general, the spatial accuracy is reduced to first order in regions near the interface, due to the IB treatment. Additionally, it is common to employ a projection method for the temporal discretization of the incompressible Navier-Stokes equations with first- or second-order accuracy in time. For a FSI problem involving the interaction between an elastic structure and a viscous fluid, the velocity is continuous at the fluid-solid interface due to the no-slip condition. The pressure and the normal derivatives of the velocity and pressure, however, have jump discontinuities across the interface. On the other hand, the solution of the FSI problem usually possesses a higher degree of regularity (smoothness) away from the fluid-structure interface, and so it appears that a numerical solver does not need to be restricted by the low-order accuracy that results primarily from the interface treatment. Thus, on top of the IB method to treat the interface. it seems natural to employ a high-order/high-resolution method for the majority of the domain that does not contain the interface, especially for those space dimensions that are not periodic where common spectral methods such as FFT are not applicable. The obvious advantage of such a hybrid algorithm is that we maintain the simple interface treatment through the IB method, meanwhile improve the accuracy for the solution near and away from the interface. Popular high-order numerical methods for time-dependent flow problems include the finite-difference or finite-volume Weighted Essentially Non-Oscillatory (WENO) methods [9, 18], the finite-element discontinuous Galerkin methods [2, 3], the spectral methods [14], and so on. In particular, WENO methods are a class of variations to the Essentially Non-Oscillatory (ENO) method [7–12]. Different from the ENO reconstruction, which is based on the information of the smoothest stencil, WENO reconstruction uses a convex combination of all candidate stencils. The first WENO method was introduced by Liu, Osher and Chan in [23], in which a third-order accurate finite-volume WENO scheme was proposed. Later, Jiang and Shu [18] designed a general framework to construct arbitrary-order accurate finite-difference WENO schemes which can be easily applied to multi-dimensional cases. Since then, WENO methods are further developed on both structured and unstructured meshes with applications in many physical problems [1, 5, 13, 16, 26, 30]. The main advantage of a WENO method is that it can attain uniform high-order accuracy in smooth regions of the solution, while maintain a sharp and essentially monotone solution transition 3
in regions with discontinuities. In this study, we propose a new hybrid algorithm, referred to as IB-WENO, by integrating a finite-difference WENO method into the IB framework for FSI computation. The remainder of this paper is organized as follows. In section 2, we present the general framework of the IB formulation. In section 3, we describe the Hodge decomposition which efficiently handles the divergence-free condition in solving incompressible fluid equations. In section 4, we provide details of the spatial discretization based on a third-order WENO scheme. Step-by-step implementation of the IB-WENO algorithm is provided in section 5. Several numerical examples are presented in section 6 to demonstrate the performance of the proposed algorithm and its various components. Finally, concluding remark and some discussion are given in section 7.
2
Immersed boundary formulation
We first describe the general formulation based on the IB method. We consider the motion of an immersed elastic boundary, which represents a solid structure, inside an incompressible viscous fluid. The boundary exert a force, referred to as the FSI force, on the fluid, and moves with the local fluid at the same velocity; i.e., the no-slip condition is assumed. We also assume that the boundary is thin such that its volume can be (approximately) treated as zero. Such an elastic boundary, for example, may refer to a fiber in two-dimensional (2D) space, and a membrane in two-dimensional or three-dimensional (3D) space. Let Ω denote the fluid domain and Γ denote the elastic boundary. The motion of the fluid in Ω is governed by the incompressible Navier-Stokes equations ( ρ ∂u + u · ∇u + ∇p = µ4u + f , ∂t ∇ · u = 0.
(2.1)
Here ρ is the density, u is the velocity, p is the pressure, µ is the (dynamic) viscosity, and f is the FSI force which acts on the immersed boundary Γ. We assume that both ρ and µ are constant. In the IB framework, the fluid motion is handled by the Eulerian approach; that is, equation (2.1) is solved on a fixed Cartesian mesh. For ease of discussion, we focus on 2D space in what follows. We let x = (x, y) denote the Eulerian coordinates for the fluid domain. Meanwhile, the motion of the elastic boundary is tracked by the Lagrangian approach. Let X(s, t) represent the position of the boundary at time t, where s is a variable to parameterize the boundary Γ. Let also F(s, t) be the elastic force density 4
at time t. If we denote the elastic energy functional as E(s, t), then we have F=−
℘E ; ℘X
(2.2)
i.e., −F is the Frechet derivative of E with the configuration X [27]. If F is known, the FSI force f (x, t) on the Eulerian coordinates can be evaluated by Z f (x, t) = F(s, t)δ x − X(s, t) ds ,
(2.3)
Γ
where δ denotes the Dirac delta function. Essentially, equation (2.3) shows that the
FSI force f (x, t) is singular; it is zero everywhere in the fluid domain except on the boundary Γ (which has zero volume). The IB method smooths out this singularity by replacing the the Dirac delta function with a discrete delta function. In 2D space, a typical discrete delta function is given by (see [27]): δh (x, y) =
1 φ(x/h)φ(y/h), h2
where h = ∆x = ∆y, and 0, m ≤ −2 √ 1 2 , 5 + 2m − −7 − 12m − 4m −2 ≤ m ≤ −1 8 √ 1 3 + 2m + 1 − 4m − 4m2 , −1 ≤ m ≤ 0 φ(m) = 8 √ 1 2 , 3 − 2m + 1 + 4m − 4m 0≤m≤1 8 1 √ 2 , 5 − 2m − −7 + 12m − 4m 1≤m≤2 8 0, 2 ≤ m
(2.4)
(2.5)
Numerically, the use of a discrete delta function represents an interpolation of the elastic force from the boundary Γ to the fluid domain Ω. With the force f computed, the fluid equation (2.1) can then be solved in the entire domain Ω using the Cartesian mesh. Once the fluid velocity is determined, the velocity of the solid boundary can be calculated using the no-slip condition. As a result, we have Z u(X(s, t), t) = u(x, t)δ(x − X(s, t)) dx.
(2.6)
Ω
Consequently, the location of the boundary can be updated by ∂X (s, t) = u(X(s, t), t). ∂t 5
(2.7)
For the implementation of the IB method, the Dirac delta function δ in equation (2.6) is again replaced by a discrete delta function, say δh , representing an interpolation of the fluid velocity from the domain Ω to the boundary Γ. As a short summary, for each cycle of computation, the IB method typically starts from the calculation of the elastic force based on equation (2.2), followed by an interpolation of the force from the Lagrangian points to the Eulerian grids using equation (2.3) and an appropriate discrete delta function. The fluid equation (2.1) is then solved in the entire domain, and the computed fluid velocity is interpolated back to the elastic boundary using equation (2.6), which subsequently updates the position of the boundary with equation (2.7).
3
Hodge decomposition
To facilitate the incorporation of a WENO method into the IB framework, we employ the Hodge decomposition to treat the incompressibility condition in the fluid equations. Let us rewrite the 2D Navier-Stokes equations in (2.1) as ( ρ(ut + F(u)x + G(u)y ) + ∇p = µ4u + f , ∇ · u = 0, where
u u= , v
u2 F(u) = , uv
We rearrange (3.1) into the following form
uv G(u) = 2 , v
f f= 1 . f2
1 1 ut + ∇ p = −F(u)x − G(u)y + (µ4u + f ). ρ ρ
(3.1)
(3.2)
(3.3)
Using the Hodge decomposition, (3.3) can be written as ut = P(R(u, f )),
(3.4)
where
1 R(u, f ) = −F(u)x − G(u)y + (µ4u + f ) (3.5) ρ and P denotes the Hodge projection operator. Details of the Hodge decomposition are provided below. The Hodge decomposition describes the decomposition of a flow field into its divergence-free and curl-free components. Let a be a vector field on a bounded domain. Then it can be uniquely decomposed in the form a = b + ∇c 6
(3.6)
for some vector b and scalar c , where b satisfies the divergence-free condition. Alternatively, We call b the Hodge projection of a and denote b = P(a) if the following two conditions are satisfied: • ∇ · b = 0, • ∇ × b = ∇ × a. If a is given, b can be computed by an appropriate numerical method. For instance, if periodic boundary conditions are assumed on a domain with uniform mesh, we can apply the Fast Fourier Transform (FFT). As a simple illustration, we assume a = (a1 , a2 ), b = (b1 , b2 ) are 2π-periodic vector fields, and a is known. Through Fourier expansions, we have X
Ny /2
X
a ˆm,kl eI(kx+ly) ,
X
X
ˆbm,kl eI(kx+ly) ,
Nx /2
am (x, y) =
(3.7)
k=−Nx /2 l=−Ny /2 Ny /2
Nx /2
bm (x, y) =
k=−Nx /2 l=−Ny /2
for m = 1, 2, where I denotes the imaginary unit, and Nx and Ny denote the numbers of points on the x and y directions, respectively. In order to find ˆbm,kl , we start with the component form of the relation between a and b, ( (b1 )x + (b2 )y = 0, (b2 )x − (b1 )y = (a2 )x − (a1 )y ,
(3.8)
and discretize (3.8) by a fourth-order central difference method. We denote the difference operators in the x and y directions by Dx and Dy , respectively. We then obtain a discrete system, ( Dx (b1 )ij + Dy (b2 )ij = 0, Dx (b2 )ij − Dy (b1 )ij = Dx (a2 )ij − Dy (a1 )ij .
(3.9)
Meanwhile, perform Fourier transform to the operators Dx and Dy to obtain I 4 1 ˆ x (k) = D sin k4x − sin 2k4x , (3.10) 4x 3 6 and
ˆ y (l) = I D 4y
4 1 sin l4y − sin 2l4y . 3 6 7
(3.11)
Equation (3.9) is then equivalent to the following system in the frequency domain, ( ˆ x (k)ˆb1,kl + D ˆ y (l)ˆb2,kl = 0, D (3.12) ˆ x (k)ˆb2,kl − D ˆ y (l)ˆb1,kl = D ˆ x (k)ˆ ˆ y (l)ˆ D a2,kl − D a1,kl , from which ˆb1,kl and ˆb2,kl can be easily calculated. To ensure the stability, an exponential filter [19, 24] is employed in our study. With the computation on the Hodge decomposition, equation (3.4) leads to a semi-discrete system, ut = L (u) ,
(3.13)
where the spatial discretization of the right-hand side of equation (3.4) is accomplished by using a WENO method (details provided in Sec. 4). The semi-discrete system (3.13) is further discretized by a TVD Runge-Kutta scheme. For instance, in our numerical simulation we have employed the second-order and third-order TVD RungeKutta methods [33] which, respectively, take the following forms: u(1) = un + ∆tL (un ), 1 1 1 un+1 = un + u(1) + ∆tL (u(1) ); 2 2 2
(3.14)
u(1) = un + ∆tL (un ), 3 1 1 u(2) = un + u(1) + ∆tL (u(1) ), 4 4 4 2 2 1 un+1 = un + u(2) + ∆tL (u(2) ). 3 3 3
(3.15)
and
4
Spatial discretization
In this section, we focus on the numerical evaluation of R(u, f ) in (3.5), which leads to the spatial discretization of the fluid equations. Notice that both the advection and diffusion terms are included. We discretize the advection terms by a high-order WENO method, and the diffusion terms by a high-order central difference method. As a result, we obtain R(uij , fij ) = −
ˆ i+1/2,j − F ˆ i−1/2,j ˆ i,j+1/2 − G ˆ i,j−1/2 1 F G − + (µDuij + fij ), 4x 4y ρ
where D = Dx + Dy denotes the fourth-order central difference operator, and 1 1 4 5 4 1 Dx uij = − ui−2,j + ui−1,j − uij + ui+1,j − ui+2,j . 4x2 12 3 2 3 12 8
(4.1)
(4.2)
ˆ i+1/2,j is consistent with the physical flux, and is Lipschitz conThe numerical flux F tinuous with respect to all its variables. Notice that ∂F 1 b b 1 ), (4.3) ≈ (F 1 − F i− 2 ,j ∂x ∆x i+ 2 ,j ∂G 1 b b (4.4) ≈ (G 1 − G i,j− 12 ). ∂y ∆y i,j+ 2 Lax-Friedrichs flux splitting is used for constructing upwinding fluxes 1 1 F+ (u) = (F(u) + α1 u), F− (u) = (F(u) − α1 u), (4.5) 2 2 1 1 G+ (u) = (G(u) + α2 u), G− (u) = (G(u) − α2 u), (4.6) 2 2 where the maximum advection speeds α1 and α2 are approximated by α1 = max |F0 (u)|(xi ,yj ) ,
α2 = max |G0 (u)|(xi ,yj ) .
i,j
(4.7)
i,j
Numerical fluxes in (4.3) and (4.4) are computed by b 1 =F b+ 1 + F b− 1 , F i+ ,j i+ ,j i+ ,j 2
2
b b+ b− G i,j+ 1 = Gi,j+ 1 + Gi,j+ 1 . 2
2
2
(4.8)
2
b+ 1 , F b− 1 , G b+ 1, The third-order WENO interpolation is utilized to approximate F i+ 2 ,j i+ 2 ,j i,j+ 2 b − 1 . For example, to compute F b + 1 , we employ F+ (ui−1,j ), F+ (ui,j ), F+ (ui+1,j ) and G i,j+ 2
i+ 2 ,j
to obtain
b + 1 = ω0 ( 1 F+ (ui,j ) + 1 F+ (ui+1,j )) + ω1 (− 1 F+ (ui−1,j ) + 3 F+ (ui,j )), F i+ 2 ,j 2 2 2 2 where α0 α1 , ω1 = , ω0 = α0 + α1 α0 + α1 2 3
α0 =
, ( + β0 )2 β0 = (F+ (ui+1,j ) − F+ (ui,j ))2 ,
1 3
α1 =
, ( + β1 )2 β1 = (F+ (ui,j ) − F+ (ui−1,j ))2 .
b − 1 , we utilize F− (ui,j ), F− (ui+1,j ), F− (ui+2,j ) and obtain To compute F i+ ,j
(4.9)
(4.10)
2
3 1 1 1 b− 1 = ω e0 ( F− (ui+1,j ) − F− (ui+2,j )) + ω e1 ( F− (ui,j ) + F− (ui+1,j )), F i+ 2 ,j 2 2 2 2 where α e0 α e1 , ω e1 = , ω e0 = α e0 + α e1 α e0 + α e1 α e0 =
1 3
( + βe0 )2
α e1 =
,
βe0 = (F− (ui+2,j ) − D− (ui+1,j ))2 ,
2 3
( + βe1 )2
,
βe1 = (F− (ui+1,j ) − F− (ui,j ))2 .
(4.11)
(4.12)
Similarly, we use numerical values G+ (ui,j−1 ), G+ (ui,j ) and G+ (ui,j+1 ) to compute b + 1 , and use G− (ui,j ), G− (ui,j+1 ) and G− (ui,j+2 ) to compute G b− 1. G i,j+ 2
i,j+ 2
9
5
Algorithm implementation
We now describe the implementation procedure of the aforementioned numerical methods that form a new algorithm for FSI computation. In general, we consider a solid structure represented by an elastic boundary immersed in an incompressible viscous fluid, which occupies a rectangular domain [xmin , xmax ] × [ymin , ymax ]. The fluid domain is partitioned by uniform Eulerian mesh xi,j = (xi , yj ), and uni,j denotes the numerical approximation of u(xi , yj , tn ), where 0 ≤ i ≤ Nx and 0 ≤ j ≤ Ny .
Let also the mesh size be h in both the x and y directions. Meanwhile, the elastic boundary is represented by Lagrangian coordinates, with X(ks , tn ) = (Xkns , Ykns ) at time tn , which is tracked from the initial position (Xk0s , Yk0s ), where 0 ≤ ks < Ns . The following steps advance the solution from t = tn to t = tn+1 .
Step 1: We first update the location of the Lagrangian points on the boundary by 1 half time-step; i.e., at t = tn+ 2 , we have n+1/2
Xks and
n+1/2
Yks
4t X 2 0≤q≤N
= Xkns +
= Ykns +
4t X 2 0≤q≤N
x
x
X
unqr δh (xq − Xkns , yr − Ykns )h2 ,
(5.1)
X
n vqr δh (xq − Xkns , yr − Ykns )h2 .
(5.2)
0≤r≤Ny
0≤r≤Ny
These correspond to the discretization of equations (2.6) and (2.7) by a half time-step. Step 2: We calculate the elastic force n+ 1 Fks 2
n+1/2
=
n+1/2
n+1/2
n+1/2
Tks +1/2 τks +1/2 − Tks −1/2 τks −1/2 4sks
n+1/2
n+1/2
(5.3)
n+1/2
at each Lagrangian point Xks = (Xks , Yks ), based on equation (2.2). Here n+1/2 n+1/2 we assume that each pair of adjacent Lagrangian points Xks and Xks +1 form a n+1/2
linar elastic fiber. Using Hooke’s law, the fiber tension Tks +1/2 is computed by n+1/2
n+1/2 Tks +1/2
= T0
n+1/2
|Xks +1 − Xks 4sks +1/2
|
!
−1 ,
(5.4)
where T0 is the elastic (or, spring) constant, and n+1/2 τks +1/2
=
n+1/2
n+1/2
n+1/2
n+1/2
Xks +1 − Xks
|Xks +1 − Xks
n+1/2
n+1/2
|
(5.5)
is the unit tangent of the fiber. Tks −1/2 and τks −1/2 are calculated in a similar way. 10
Step 3: We distribute the elastic force from the Lagrangian points to the Eulerian grids by X n+1/2 n+1/2 fi,j = Fks δh (xi − Xkns , yj − Ykns )∆sks (5.6) 0≤ks
based on equation (2.3). Step 4: We use the Hodge decomposition and WENO reconstruction to solve the incompressible Navier-Stokes equations (2.1) on the Eulerian mesh. A second- or third-order TVD Runge-Kutta method is employed for the temporal discretization (see equations 3.14 and 3.15) to generate the fluid velocity un+1 at t = tn+1 . Step 5: We then update the location of the elastic boundary at t = tn+1 by Xkn+1 = Xkns + 4t s
X
X
n+1/2
uij
0≤i≤Nx 0≤j≤Ny
n+1/2
δh (xi − Xks
n+1/2
, y j − Yk s
)h2
(5.7)
)h2 .
(5.8)
and Ykn+1 = Ykns + 4t s
6
X
X
n+1/2
vij
0≤i≤Nx 0≤j≤Ny
n+1/2
δh (xi − Xks
n+1/2
, yj − Yks
Numerical examples
Our proposed algorithm involves several components based on different numerical methods. In what follows, we provide a number of simulation examples to thoroughly test these numerical components, as well as the entire algorithm. Example 6.1. We first test the accuracy of our Hodge projection described in Section 3. We assume that the vector a(x, y) in (3.6) takes the form of sin2 x cos y + 2 sin (x + y) cos (x + y) cos y a= −2 sin x cos x sin y + 2 sin (x + y) cos (x + y) cos y − sin2 (x + y) sin y on the domain [0, 2π] × [0, 2π]. It is easy to verify that the exact solution for the
vector b is given by
sin2 x cos y b= . −2 sin x cos x sin y
We numerically solve b through (3.7)-(3.12), and the errors of b1 and b2 in L∞ norm (among all the grid points) are presented in Table 6.1, where N = Nx = Ny . We clearly observe that uniform fourth order is attained, consistent with the fourthorder finite difference method involved in the discretization. Similar results in L1 and L2 norms are also obtained (not shown here). 11
N L∞ error of b1 64 × 64 8.26E-5 128 × 128 5.21E-6 256 × 256 3.26E-7 512 × 512 2.04E-8 1024 × 1024 1.28E-9 2048 × 2048 7.97E-11
Order L∞ error of b2 6.04E-5 3.99 3.81E-6 4.00 2.39E-7 4.00 1.49E-8 4.00 9.34E-10 4.00 5.83E-11
Order 3.99 4.00 4.00 4.00 4.00
Table 6.1: Accuracy test for the Hodge decomposition. Example 6.2. Next, we verify the accuracy of our WENO reconstruction when solving the following compressible Navier-Stokes equations ( ρ(ut + F(u)x + G(u)y ) = µ4u + f , u(x, y, 0) = u0 (x, y)
(6.1)
on the computational domain [0, 2π] × [0, 2π] with periodic boundary conditions. The initial velocity is taken as
sin2 (x + 2y) cos(3y) u0 (x, y) = . cos3 (2x − 3y) sin (x + y)
The density is set as ρ = 1 and the viscosity is µ = 0.05. The exact solution to (6.1) can be found as sin2 (x + 2y) cos(3y)e−2µt u(x, y, t) = . 2 cos3 (2x − 3y + t) sin (x + y)e−3µt Here we apply the third-order WENO scheme for the convection term and fourthorder central difference scheme for the diffusion term. The temporal discretization is based on the third-order TVD Runge-Kutta method. Numerical results with the third-order WENO reconstruction at the end time T = 2 are given in Table 6.2, for both components of the velocity, u and v, in L∞ norm (among all the grid points). It is clear to see that, as the mesh size if refined, third-order accuracy is attained. Similar patterns in L1 and L2 norms are also observed (not presented here). Example 6.3. We proceed to check the accuracy of our combined Hodge decomposition and WENO method, using the Taylor-Green vortex test problem based on the incompressible Navier-Stokes equations (3.1). On the domain [0, 2π] × [0, 2π], an exact solution to the problem is given by 2t − −2k Re , u(x, y, t) = − cos(kx) sin(ke)e 2k2 t (6.2) v(x, y, t) = sin(kx) cos(ky)e− Re , 2t 4k p(x, y, t) = − 1 (cos(kx) + cos(ky))e− Re , 4
12
N L∞ error of u 80 × 80 5.36E-2 160 × 160 9.02E-3 320 × 320 1.33E-3 640 × 640 1.63E-4 1280 × 1280 1.54E-5
Order L∞ error of v 6.18E-2 2.57 8.57E-3 2.76 1.12E-3 3.04 1.37E-4 3.39 1.48E-5
Order 2.85 2.93 3.04 3.20
Table 6.2: Accuracy test for the WENO reconstruction. where k = 4 and µ = 0.05. As before, we use the third-order TVD Runge-Kutta method for the time marching, the third-order WENO scheme for the convective term, and the fourth-order central difference scheme for the viscous term. Numerical results at T = 2 are presented in Table 6.3 for both components of the velocity in L∞ norm. Again we observe that as the mesh size is refined, third-order accuracy is achieved, consistent with the third-order WENO scheme. Similar patterns hold for the errors in L1 and L2 norms. N L∞ error of u 80 × 80 3.80E-3 160 × 160 6.24E-4 320 × 320 7.90E-5 640 × 640 9.07E-6 1280 × 1280 9.49E-7
Order L∞ error of v 3.80E-3 2.60 6.24E-4 2.98 7.90E-5 3.12 9.07E-6 3.26 9.49E-7
Order 2.60 2.98 3.12 3.26
Table 6.3: Accuracy test for the combined Hodge-WENO method. Example 6.4. We further test the performance of our combined Hodge-WENO method when the solution is not smooth. We consider the double shear layer problem with the governing equation (3.1) and computational domain [0, 2π] × [0, 2π]. The initial condition is
u(x, y, 0) =
(
tanh ((y − π/2)/ρ),
y≤π
tanh ((3π/2 − y)/ρ),
y>π
and v(x, y, 0) = δ sin x. The boundary condition is periodic. The parameters are taken as µ = 0, ρ = π/15, δ = 0.05. Numerical approximations for the vorticity contours at T = 8 are plotted in Figure (6.1) with grid points of 128 × 128 and 512 × 512, respectively. The results agree well with those in literature (e.g., [32]). 13
Figure 6.1: Vorticity contours in Example 6.4. Left: N = 128. Right: N = 512. Example 6.5. Now we proceed to test our entire IB-WENO algorithm that integrates the IB method, the Hodge decomposition, and the WENO scheme, using the canonical problem of a stretched rubber band relaxing in an incompressible viscous fluid (see Figure 6.2). At its resting state, the rubber band is a circle with a radius r0 = 0.15. Initially, the rubber band is stretched to an ellipse, centered at (Ox , Oy ) = (0.5, 0.5), whose major semi-axis and minor semi-axis are a = 0.4 and b = 0.15, respectively. Both the fluid and the elastic boundary (i.e., the rubber band) are stationary at time t = 0. Once released, the boundary is set into motion due to the restoring elastic force and, since the fluid is incompressible, the boundary converges to a final circle which minimizes the energy and which has the same area as that of the initial ellipse. The √ radius of the final steady-state circle is re = ab ≈ 0.245. We set the fluid density ρ = 1 and the viscosity µ = 0.01. In our numerical simulation, the initial Lagrangian points on the boundary are prescribed as
2πks 0 , Xks = Ox + a cos Ns (6.3) 2πks 0 Yks = Oy + b sin , Ns where 0 ≤ ks ≤ Ns − 1 and Ns is the number of Lagrangian points. To reduce the
possible leakage of the volume (i.e., the area in 2D case), we impose the restriction that |X0ks − X0ks −1 | < 14
h , 2
where h = ∆x = ∆y is the fluid mesh size.
Figure 6.2: The stretched rubber band problem in Example 6.5. The solid, dashed, and dash-dot lines represent the initial, final, and resting positions of the rubber band, respectively. We have applied both the standard IB method (implemented with FFT) and the proposed IB-WENO algorithm to this problem. We find that both algorithms achieve first-order accuracy. This is expected, since the errors at the moving interface (where discontinuity occurs) dominate for the computation of the entire FSI problem, and the IB treatment of the interface reduces the overall accuracy to first order. In order to better compare the two algorithms, we look into the location and shape of the interface profile, as well as the area enclosed, throughout the process. Due to numerical errors, the computed interface does not approach a perfect, final circle. Thus, to quantify such errors, we calculate the minimum and maximum distances between the center (Ox , Oy ) and all the computed boundary points, at each time step t = tn , by q n rmin = min (Xkns − Ox )2 + (Ykns − Oy )2 , 0≤k≤Ns −1 q (6.4) n n n 2 2 rmax = max (Xks − Ox ) + (Yks − Oy ) , 0≤k≤Ns −1
and compare with the exact value re ≈ 0.245. Figure 6.3 shows the results by using the standard IB method and the IB-WENO method, where we observe that the values
of rmin and rmax generated by the standard IB method appear to diverge from the exact solution re over time, which leads to the “volume leakage” [20,28]. In contrast, 15
the results produced by the IB-WENO algorithm have better accuracy and stability over time.
Figure 6.3: Values of rmin and rmax versus time. Solid lines (red and green) represent the results generated by the standard IB method, and dashed lines (blue and black) represent the results generated by the IB-WENO method. The dash-dot line (pink) represents the exact solution. In addition, if we regard the computed interface as an “approximate circle” at each time tn , we may calculate the radius of this approximate circle by the average distance between the center (Ox , Oy ) and all the computed boundary points, r¯n =
Ns −1 q 1 X (Xkns − Ox )2 + (Ykns − Oy )2 . Ns k =0
(6.5)
s
Consequently, the area enclosed by the computed interface at t = tn can be approximately determined as S n = π(¯ rn )2 , and we compare this value with the exact area Se = πab (which is the area enclosed by the initial ellipse). Figure 6.4 displays the results by the two methods: the standard IB, and the IB-WENO. We observe that the standard IB method generates smaller error initially, but the error keeps increasing over time (another indication of the volume leakage). In contrast, the IB-WENO 16
method seems to yield more stable and consistent result, with better accuracy eventually.
Figure 6.4: |S n − Se | versus time. The (red) solid line represents the result of the IB-WENO method, and the (green) dashed line represents the result of the standard IB method.
7
Discussion
We have presented a new approach that integrates the IB method, the Hodge decomposition, and a finite-difference WENO scheme for FSI computation. With this hybrid algorithm, named IB-WENO, we maintain the simple treatment of the fluidstructure interface, while achieve a high-order accuracy for regions near and away from the interface where the solution is sufficiently smooth. This work represents a pilot study in combining the IB and WENO methods. The combined algorithm is easy to implement, applicable to both periodic and non-periodic domains, and our numerical results demonstrate the performance of the algorithm and its various components. In particular, its overall accuracy is comparable to or better than that of the standard combination of the IB method and FFT, when implemented in periodic 17
space dimensions. Additionally, our results in the stretched rubber band test problem shows a better volume conservation than that of the standard IB method. In a recent study [36], the authors proposed a method that combines a WENO scheme and the ghost-cell method (a variation of the original IB method) to compute shock-structure interaction, where the fluid flow is described by the Euler equations; i.e., fluid viscosity is neglected. The strength of their method is that it may handle high-speed flow with strong discontinuities. The formulation and implementation of their method, however, are much more complex than ours. Currently, our integrated algorithm achieves first-order accuracy at the interface because of the IB treatment, and third-order accuracy in other parts of the domain because of the third-order WENO scheme employed. There is an apparent mismatch of the orders of accuracy for these two components in the hybrid algorithm. On one hand, the degradation of accuracy at the interface is inevitable due to the jump discontinuity of the solution. The majority of the immersed type methods achieve first-order accuracy, with only a few exceptions (such as the immersed interface method [20,21]) that require additional, highly nontrivial procedures to boost the accuracy at the interface. On the other hand, a different direction to alleviate the order “mismatch” and improve the accuracy at the interface is to employ a fine mesh in the small region that contains the interface, and use a coarse mesh for the remaining parts of the domain. In other words, the low-order method for the interface treatment (such as our IB method) is based on smaller grid size, whereas the high-order method for the majority of the domain (such as our WENO scheme) is based on larger grid size, so as to reach a balance for the accuracy (not the order of accuracy, but rather the magnitude of the numerical error) between the two methods. Consequently, a local mesh refinement procedure [31] would be a natural choice to achieve this goal. We are currently working along this line and implementing the immersed boundary method with adaptive mesh refinement.
Acknowledgments Y.L. acknowledges partial support from the Simons Foundation (No. 426993). M.Z. acknowledges partial support from the National Natural Science Foundation of China (No. 11471305). The work of J.W. was partially supported by the National Science Foundation (Nos. 1319078 and 1520672), and by the Tennessee Higher Education Commission through the CEACSE program. The authors thank the editor for handling this paper and the reviewer for helpful comments. 18
References [1] Marcos Castro, Bruno Costa, and Wai Sun Don. High order weighted essentially non-oscillatory WENO-Z schemes for hyperbolic conservation laws. Journal of Computational Physics, 230(5):1766–1792, 2011. [2] Bernardo Cockburn and Chi-Wang Shu. TVB runge-kutta local projection discontinuous Galerkin finite element method for conservation laws. ii. general framework. Mathematics of computation, 52(186):411–435, 1989. [3] Bernardo Cockburn and Chi-Wang Shu. The runge–kutta discontinuous galerkin method for conservation laws v: multidimensional systems. Journal of Computational Physics, 141(2):199–224, 1998. [4] Earl H Dowell and Kenneth C Hall. Modeling of fluid-structure interaction. Annual Review of Fluid Mechanics, 33(1):445–490, 2001. [5] Travis C Fisher, Mark H Carpenter, Nail K Yamaleev, and Steven H Frankel. Boundary closures for fourth-order energy stable weighted essentially non-oscillatory finite-difference schemes. Journal of Computational Physics, 230(10):3727–3752, 2011. [6] Robert D Guy and David A Hartenstine. On the accuracy of direct forcing immersed boundary methods with projection methods. Journal of Computational Physics, 229(7):2479–2496, 2010. [7] Ami Harten. High resolution schemes for hyperbolic conservation laws. Journal of computational physics, 49(3):357–393, 1983. [8] Ami Harten. Preliminary results on the extension of ENO schemes to twodimensional problems. Nonlinear Hyperbolic Problems, C. Carasso, P.-A. Raviart, and D. Serre, eds., LNM, 1270:23–40, 1987. [9] Ami Harten. ENO schemes with subcell resolution. Journal of Computational Physics, 83(1):148–184, 1989. [10] Ami Harten, Bjorn Engquist, Stanley Osher, and Sukumar R Chakravarthy. Uniformly high order accurate essentially non-oscillatory schemes, III. Journal of computational physics, 71(2):231–303, 1987. [11] Ami Harten and Stanley Osher. Uniformly high-order accurate nonoscillatory schemes. I. SIAM Journal on Numerical Analysis, 24(2):279–309, 1987. 19
[12] Ami Harten, Stanley Osher, Bj¨orn Engquist, and Sukumar R Chakravarthy. Some results on uniformly high-order accurate essentially nonoscillatory schemes. Applied Numerical Mathematics, 2(3-5):347–377, 1986. [13] Andrew K Henrick, Tariq D Aslam, and Joseph M Powers. Mapped weighted essentially non-oscillatory schemes: achieving optimal order near critical points. Journal of Computational Physics, 207(2):542–567, 2005. [14] Jan S Hesthaven, Sigal Gottlieb, and David Gottlieb. Spectral methods for time-dependent problems. Cambridge University Press, 21, 2007. [15] Gene Hou, Jin Wang, and Anita Layton. Numerical methods for fluid-structure interaction—a review. Communications in Computational Physics, 12(2):337– 377, 2012. [16] Changqing Hu and Chi-Wang Shu. Weighted essentially non-oscillatory schemes on triangular meshes. Journal of Computational Physics, 150(1):97–127, 1999. [17] Bj¨orn H¨ ubner, Elmar Walhorn, and Dieter Dinkler. A monolithic approach to fluid–structure interaction using space–time finite elements. Computer methods in applied mechanics and engineering, 193(23):2087–2104, 2004. [18] Guang-Shan Jiang and Chi-Wang Shu. Efficient implementation of weighted ENO schemes. Journal of computational physics, 126(1):202–228, 1996. [19] David A Kopriva. A practical assessment of spectral accuracy for hyperbolic problems with discontinuities. Journal of Scientific Computing, 2(3):249–262, 1987. [20] Randall J LeVeque and Zhilin Li. Immersed interface methods for stokes flow with elastic boundaries or surface tension. SIAM Journal on Scientific Computing, 18(3):709–735, 1997. [21] Zhilin Li. An overview of the immersed interface method and its applications. Taiwanese journal of mathematics, pages 1–49, 2003. [22] Wing Kam Liu and Shaoqiang Tang. Mathematical foundations of the immersed finite element method. Computational Mechanics, 39(3):211–222, 2007. [23] Xu-Dong Liu, Stanley Osher, and Tony Chan.
Weighted essentially non-
oscillatory schemes. Journal of computational physics, 115(1):200–212, 1994. 20
[24] Andrew Majda, James McDonough, and Stanley Osher. The Fourier method for nonsmooth initial data. Mathematics of Computation, 32(144):1041–1081, 1978. [25] J Mohd-Yusof. Combined immersed boundary/b-spline methods for simulations of flow in complex geometries. Annual Research Briefs, Center for Turbulence Research, Stanford University, pages 317–327, 1999. [26] Mohammad Motamed, Colin B Macdonald, and Steven J Ruuth. On the linear stability of the fifth-order weno discretization. Journal of Scientific Computing, 47(2):127–149, 2011. [27] Charles S Peskin. The immersed boundary method. Acta numerica, 11:479–517, 2002. [28] Charles S Peskin and Beth Feller Printz. Improved volume conservation in the computation of flows with immersed elastic boundaries. Journal of computational physics, 105(1):33–46, 1993. [29] PB Ryzhakov, Riccardo Rossi, SR Idelsohn, and E O˜ nate. A monolithic lagrangian approach for fluid–structure interaction problems. Computational mechanics, 46(6):883–899, 2010. [30] Susana Serna and Antonio Marquina.
Power ENO methods: a fifth-order
accurate weighted power ENO method.
Journal of Computational Physics,
194(2):632–658, 2004. [31] Chaopeng Shen, Jing-Mei Qiu, and Andrew Christlieb. Adaptive mesh refinement based on high order finite difference weno scheme for multi-scale simulations. Journal of Computational Physics, 230(10):3780–3802, 2011. [32] Chi-Wang Shu.
Essentially non-oscillatory and weighted essentially non-
oscillatory schemes for hyperbolic conservation laws. Advanced numerical approximation of nonlinear hyperbolic equations, pages 325–432, 1998. [33] Chi-Wang Shu and Stanley Osher.
Efficient implementation of essentially
non-oscillatory shock-capturing schemes. Journal of Computational Physics, 77(2):439–471, 1988. [34] J Steindorf and HG Matthies. Numerical efficiency of different partitioned methods for fluid-structure interaction. ZAMM-Journal of Applied Mathematics and
21
Mechanics/Zeitschrift f¨ ur Angewandte Mathematik und Mechanik, 80(S2):557– 558, 2000. [35] Jan Vierendeels, Kris Dumont, and PR Verdonck. A partitioned strongly coupled fluid-structure interaction method to model heart valve dynamics. Journal of Computational and Applied Mathematics, 215(2):602–609, 2008. [36] Min Xu, Tao Yang, and Mingjun Wei. Implementation of immersed boundary method in weno scheme to simulate shock-structure interaction. ASME 2017 Fluids Engineering Division Summer Meeting, pages V01BT11A013– V01BT11A013, 2017. [37] LT Zhang and M Gay. Immersed finite element method for fluid-structure interactions. Journal of Fluids and Structures, 23(6):839–857, 2007.
22
PII: DOI: Reference:
S0377-0427(18)30635-6 https://doi.org/10.1016/j.cam.2018.10.028 CAM 11973
To appear in:
Journal of Computational and Applied Mathematics
Received date : 12 January 2018 Revised date : 5 July 2018 Please cite this article as: Z. Cheng, et al., IB-WENO method for incompressible flow with elastic boundaries, Journal of Computational and Applied Mathematics (2018), https://doi.org/10.1016/j.cam.2018.10.028 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
IB-WENO method for incompressible flow with elastic boundaries Ziqiang Cheng a
Yuan Liu b a School
Mengping Zhang a
Jin Wang c,∗
of Mathematical Sciences,
University of Science and Technology of China, Hefei, Anhui 230026, China b Department
of Mathematics and Statistics,
Mississippi State University, Mississippi State, MS 39762, USA c Department
of Mathematics, University of Tennessee at Chattanooga, Chattanooga, TN 37403, USA ∗ Email:
[email protected]
Abstract We present a new approach to compute the interaction between an elastic boundary and an incompressible viscous fluid, by integrating the immersed boundary method and a finite difference WENO scheme. The hybrid algorithm maintains simple and efficient treatment of the moving interface based on the immersed boundary formulation, and achieves high-order accuracy, with the WENO method and the Hodge decomposition technique employed, in the majority of the fluid domain that does not contain the interface. We demonstrate the performance of the proposed algorithm and its various numerical components through a number of computational examples. Key Words: Immersed boundary formulation; Hodge decomposition; WENO method; Incompressible flow.
1
Introduction
The study of fluid-structure interaction (FSI), an emerging field concerned with the interaction between viscous fluid flows and immersed solid structures, has been fast growing in recent decades. FSI often involves strongly nonlinear, multi-physics phenomena, with application to a wide range of scientific and engineering disciplines [4, 15]. Owing to the difficulties in mathematical analysis and limitations in physical experiments for these complex problems, FSI research and development largely rely on scientific computing. 1
There have been a large number of numerical methods developed for the computation of FSI problems. These methods can be broadly classified into the monolithic approach [17, 29], where a unified system represents the entire FSI problem and a single algorithm solves the fluid and structural dynamics simultaneously, and the partitioned approach [34, 35] that treats the fluid and the structure as two systems which can be computed separately in their respective domains and then connected at the interface. Compared to the monolithic counterpart, a partitioned scheme typically has a higher degree of flexibility in algorithm development and implementation [15]. Another common classification of FSI methods is based on the types of meshes employed. Conforming-mesh algorithms, which include Lagrangian methods, arbitrary Lagrangian-Eulerian methods, reduced-order modeling techniques, etc., use meshes that fit the body of the moving and deformable structure. Consequently, a meshupdating procedure is normally required at each time step for these methods. In contrast, algorithms based on non-conforming meshes employ Cartesian background meshes that remain fixed throughout the course of computation, offering advantages such as improved efficiency, ease of implementation, and natural handling of large structural deformation. Among these numerical techniques, the immersed type methods, originated from the immersed boundary (IB) method [27] invented by Peskin in 1970s, have achieved much success and popularity. These methods belong to the non-conforming mesh category, and are typically based on the partitioned approach. An immersed method normally represents the effects of the immersed solid structure by adding a forcing term to the right-hand side of the fluid momentum equation. As a result, the fluid motion can be computed on fixed grids in the entire computational domain, eliminating the need of costly re-meshing procedures. Since Peskin’s pioneering work, many variations and improvements have been made to the family of immersed methods; significant examples include the immersed interface method [20, 21], the immersed finite element method [22,37], and the direct-forcing method [6,25]. Still, the original IB method remains a very popular option in FSI computation, due to its simple formulation, straightforward implementation, and applicability to a wide range of FSI problems. Nevertheless, the IB method has several drawbacks. To succeed, it has to utilize a discrete delta function and represent the sharp interface by a band of cells with non-zero thickness. Thus the method is inherently first-order accurate in space. Moreover, it has been observed that the immersed boundary method may suffer from leakage problem [20, 28], although mass conservation can be improved via the use 2
of divergence-free finite-difference operators [28]. Common practice of using the IB method in FSI computation is to apply the Fast Fourier Transform (FFT) in the periodic space dimension(s), and utilize standard finite difference approximations in those space dimensions that are not periodic. In general, the spatial accuracy is reduced to first order in regions near the interface, due to the IB treatment. Additionally, it is common to employ a projection method for the temporal discretization of the incompressible Navier-Stokes equations with first- or second-order accuracy in time. For a FSI problem involving the interaction between an elastic structure and a viscous fluid, the velocity is continuous at the fluid-solid interface due to the no-slip condition. The pressure and the normal derivatives of the velocity and pressure, however, have jump discontinuities across the interface. On the other hand, the solution of the FSI problem usually possesses a higher degree of regularity (smoothness) away from the fluid-structure interface, and so it appears that a numerical solver does not need to be restricted by the low-order accuracy that results primarily from the interface treatment. Thus, on top of the IB method to treat the interface. it seems natural to employ a high-order/high-resolution method for the majority of the domain that does not contain the interface, especially for those space dimensions that are not periodic where common spectral methods such as FFT are not applicable. The obvious advantage of such a hybrid algorithm is that we maintain the simple interface treatment through the IB method, meanwhile improve the accuracy for the solution near and away from the interface. Popular high-order numerical methods for time-dependent flow problems include the finite-difference or finite-volume Weighted Essentially Non-Oscillatory (WENO) methods [9, 18], the finite-element discontinuous Galerkin methods [2, 3], the spectral methods [14], and so on. In particular, WENO methods are a class of variations to the Essentially Non-Oscillatory (ENO) method [7–12]. Different from the ENO reconstruction, which is based on the information of the smoothest stencil, WENO reconstruction uses a convex combination of all candidate stencils. The first WENO method was introduced by Liu, Osher and Chan in [23], in which a third-order accurate finite-volume WENO scheme was proposed. Later, Jiang and Shu [18] designed a general framework to construct arbitrary-order accurate finite-difference WENO schemes which can be easily applied to multi-dimensional cases. Since then, WENO methods are further developed on both structured and unstructured meshes with applications in many physical problems [1, 5, 13, 16, 26, 30]. The main advantage of a WENO method is that it can attain uniform high-order accuracy in smooth regions of the solution, while maintain a sharp and essentially monotone solution transition 3
in regions with discontinuities. In this study, we propose a new hybrid algorithm, referred to as IB-WENO, by integrating a finite-difference WENO method into the IB framework for FSI computation. The remainder of this paper is organized as follows. In section 2, we present the general framework of the IB formulation. In section 3, we describe the Hodge decomposition which efficiently handles the divergence-free condition in solving incompressible fluid equations. In section 4, we provide details of the spatial discretization based on a third-order WENO scheme. Step-by-step implementation of the IB-WENO algorithm is provided in section 5. Several numerical examples are presented in section 6 to demonstrate the performance of the proposed algorithm and its various components. Finally, concluding remark and some discussion are given in section 7.
2
Immersed boundary formulation
We first describe the general formulation based on the IB method. We consider the motion of an immersed elastic boundary, which represents a solid structure, inside an incompressible viscous fluid. The boundary exert a force, referred to as the FSI force, on the fluid, and moves with the local fluid at the same velocity; i.e., the no-slip condition is assumed. We also assume that the boundary is thin such that its volume can be (approximately) treated as zero. Such an elastic boundary, for example, may refer to a fiber in two-dimensional (2D) space, and a membrane in two-dimensional or three-dimensional (3D) space. Let Ω denote the fluid domain and Γ denote the elastic boundary. The motion of the fluid in Ω is governed by the incompressible Navier-Stokes equations ( ρ ∂u + u · ∇u + ∇p = µ4u + f , ∂t ∇ · u = 0.
(2.1)
Here ρ is the density, u is the velocity, p is the pressure, µ is the (dynamic) viscosity, and f is the FSI force which acts on the immersed boundary Γ. We assume that both ρ and µ are constant. In the IB framework, the fluid motion is handled by the Eulerian approach; that is, equation (2.1) is solved on a fixed Cartesian mesh. For ease of discussion, we focus on 2D space in what follows. We let x = (x, y) denote the Eulerian coordinates for the fluid domain. Meanwhile, the motion of the elastic boundary is tracked by the Lagrangian approach. Let X(s, t) represent the position of the boundary at time t, where s is a variable to parameterize the boundary Γ. Let also F(s, t) be the elastic force density 4
at time t. If we denote the elastic energy functional as E(s, t), then we have F=−
℘E ; ℘X
(2.2)
i.e., −F is the Frechet derivative of E with the configuration X [27]. If F is known, the FSI force f (x, t) on the Eulerian coordinates can be evaluated by Z f (x, t) = F(s, t)δ x − X(s, t) ds ,
(2.3)
Γ
where δ denotes the Dirac delta function. Essentially, equation (2.3) shows that the
FSI force f (x, t) is singular; it is zero everywhere in the fluid domain except on the boundary Γ (which has zero volume). The IB method smooths out this singularity by replacing the the Dirac delta function with a discrete delta function. In 2D space, a typical discrete delta function is given by (see [27]): δh (x, y) =
1 φ(x/h)φ(y/h), h2
where h = ∆x = ∆y, and 0, m ≤ −2 √ 1 2 , 5 + 2m − −7 − 12m − 4m −2 ≤ m ≤ −1 8 √ 1 3 + 2m + 1 − 4m − 4m2 , −1 ≤ m ≤ 0 φ(m) = 8 √ 1 2 , 3 − 2m + 1 + 4m − 4m 0≤m≤1 8 1 √ 2 , 5 − 2m − −7 + 12m − 4m 1≤m≤2 8 0, 2 ≤ m
(2.4)
(2.5)
Numerically, the use of a discrete delta function represents an interpolation of the elastic force from the boundary Γ to the fluid domain Ω. With the force f computed, the fluid equation (2.1) can then be solved in the entire domain Ω using the Cartesian mesh. Once the fluid velocity is determined, the velocity of the solid boundary can be calculated using the no-slip condition. As a result, we have Z u(X(s, t), t) = u(x, t)δ(x − X(s, t)) dx.
(2.6)
Ω
Consequently, the location of the boundary can be updated by ∂X (s, t) = u(X(s, t), t). ∂t 5
(2.7)
For the implementation of the IB method, the Dirac delta function δ in equation (2.6) is again replaced by a discrete delta function, say δh , representing an interpolation of the fluid velocity from the domain Ω to the boundary Γ. As a short summary, for each cycle of computation, the IB method typically starts from the calculation of the elastic force based on equation (2.2), followed by an interpolation of the force from the Lagrangian points to the Eulerian grids using equation (2.3) and an appropriate discrete delta function. The fluid equation (2.1) is then solved in the entire domain, and the computed fluid velocity is interpolated back to the elastic boundary using equation (2.6), which subsequently updates the position of the boundary with equation (2.7).
3
Hodge decomposition
To facilitate the incorporation of a WENO method into the IB framework, we employ the Hodge decomposition to treat the incompressibility condition in the fluid equations. Let us rewrite the 2D Navier-Stokes equations in (2.1) as ( ρ(ut + F(u)x + G(u)y ) + ∇p = µ4u + f , ∇ · u = 0, where
u u= , v
u2 F(u) = , uv
We rearrange (3.1) into the following form
uv G(u) = 2 , v
f f= 1 . f2
1 1 ut + ∇ p = −F(u)x − G(u)y + (µ4u + f ). ρ ρ
(3.1)
(3.2)
(3.3)
Using the Hodge decomposition, (3.3) can be written as ut = P(R(u, f )),
(3.4)
where
1 R(u, f ) = −F(u)x − G(u)y + (µ4u + f ) (3.5) ρ and P denotes the Hodge projection operator. Details of the Hodge decomposition are provided below. The Hodge decomposition describes the decomposition of a flow field into its divergence-free and curl-free components. Let a be a vector field on a bounded domain. Then it can be uniquely decomposed in the form a = b + ∇c 6
(3.6)
for some vector b and scalar c , where b satisfies the divergence-free condition. Alternatively, We call b the Hodge projection of a and denote b = P(a) if the following two conditions are satisfied: • ∇ · b = 0, • ∇ × b = ∇ × a. If a is given, b can be computed by an appropriate numerical method. For instance, if periodic boundary conditions are assumed on a domain with uniform mesh, we can apply the Fast Fourier Transform (FFT). As a simple illustration, we assume a = (a1 , a2 ), b = (b1 , b2 ) are 2π-periodic vector fields, and a is known. Through Fourier expansions, we have X
Ny /2
X
a ˆm,kl eI(kx+ly) ,
X
X
ˆbm,kl eI(kx+ly) ,
Nx /2
am (x, y) =
(3.7)
k=−Nx /2 l=−Ny /2 Ny /2
Nx /2
bm (x, y) =
k=−Nx /2 l=−Ny /2
for m = 1, 2, where I denotes the imaginary unit, and Nx and Ny denote the numbers of points on the x and y directions, respectively. In order to find ˆbm,kl , we start with the component form of the relation between a and b, ( (b1 )x + (b2 )y = 0, (b2 )x − (b1 )y = (a2 )x − (a1 )y ,
(3.8)
and discretize (3.8) by a fourth-order central difference method. We denote the difference operators in the x and y directions by Dx and Dy , respectively. We then obtain a discrete system, ( Dx (b1 )ij + Dy (b2 )ij = 0, Dx (b2 )ij − Dy (b1 )ij = Dx (a2 )ij − Dy (a1 )ij .
(3.9)
Meanwhile, perform Fourier transform to the operators Dx and Dy to obtain I 4 1 ˆ x (k) = D sin k4x − sin 2k4x , (3.10) 4x 3 6 and
ˆ y (l) = I D 4y
4 1 sin l4y − sin 2l4y . 3 6 7
(3.11)
Equation (3.9) is then equivalent to the following system in the frequency domain, ( ˆ x (k)ˆb1,kl + D ˆ y (l)ˆb2,kl = 0, D (3.12) ˆ x (k)ˆb2,kl − D ˆ y (l)ˆb1,kl = D ˆ x (k)ˆ ˆ y (l)ˆ D a2,kl − D a1,kl , from which ˆb1,kl and ˆb2,kl can be easily calculated. To ensure the stability, an exponential filter [19, 24] is employed in our study. With the computation on the Hodge decomposition, equation (3.4) leads to a semi-discrete system, ut = L (u) ,
(3.13)
where the spatial discretization of the right-hand side of equation (3.4) is accomplished by using a WENO method (details provided in Sec. 4). The semi-discrete system (3.13) is further discretized by a TVD Runge-Kutta scheme. For instance, in our numerical simulation we have employed the second-order and third-order TVD RungeKutta methods [33] which, respectively, take the following forms: u(1) = un + ∆tL (un ), 1 1 1 un+1 = un + u(1) + ∆tL (u(1) ); 2 2 2
(3.14)
u(1) = un + ∆tL (un ), 3 1 1 u(2) = un + u(1) + ∆tL (u(1) ), 4 4 4 2 2 1 un+1 = un + u(2) + ∆tL (u(2) ). 3 3 3
(3.15)
and
4
Spatial discretization
In this section, we focus on the numerical evaluation of R(u, f ) in (3.5), which leads to the spatial discretization of the fluid equations. Notice that both the advection and diffusion terms are included. We discretize the advection terms by a high-order WENO method, and the diffusion terms by a high-order central difference method. As a result, we obtain R(uij , fij ) = −
ˆ i+1/2,j − F ˆ i−1/2,j ˆ i,j+1/2 − G ˆ i,j−1/2 1 F G − + (µDuij + fij ), 4x 4y ρ
where D = Dx + Dy denotes the fourth-order central difference operator, and 1 1 4 5 4 1 Dx uij = − ui−2,j + ui−1,j − uij + ui+1,j − ui+2,j . 4x2 12 3 2 3 12 8
(4.1)
(4.2)
ˆ i+1/2,j is consistent with the physical flux, and is Lipschitz conThe numerical flux F tinuous with respect to all its variables. Notice that ∂F 1 b b 1 ), (4.3) ≈ (F 1 − F i− 2 ,j ∂x ∆x i+ 2 ,j ∂G 1 b b (4.4) ≈ (G 1 − G i,j− 12 ). ∂y ∆y i,j+ 2 Lax-Friedrichs flux splitting is used for constructing upwinding fluxes 1 1 F+ (u) = (F(u) + α1 u), F− (u) = (F(u) − α1 u), (4.5) 2 2 1 1 G+ (u) = (G(u) + α2 u), G− (u) = (G(u) − α2 u), (4.6) 2 2 where the maximum advection speeds α1 and α2 are approximated by α1 = max |F0 (u)|(xi ,yj ) ,
α2 = max |G0 (u)|(xi ,yj ) .
i,j
(4.7)
i,j
Numerical fluxes in (4.3) and (4.4) are computed by b 1 =F b+ 1 + F b− 1 , F i+ ,j i+ ,j i+ ,j 2
2
b b+ b− G i,j+ 1 = Gi,j+ 1 + Gi,j+ 1 . 2
2
2
(4.8)
2
b+ 1 , F b− 1 , G b+ 1, The third-order WENO interpolation is utilized to approximate F i+ 2 ,j i+ 2 ,j i,j+ 2 b − 1 . For example, to compute F b + 1 , we employ F+ (ui−1,j ), F+ (ui,j ), F+ (ui+1,j ) and G i,j+ 2
i+ 2 ,j
to obtain
b + 1 = ω0 ( 1 F+ (ui,j ) + 1 F+ (ui+1,j )) + ω1 (− 1 F+ (ui−1,j ) + 3 F+ (ui,j )), F i+ 2 ,j 2 2 2 2 where α0 α1 , ω1 = , ω0 = α0 + α1 α0 + α1 2 3
α0 =
, ( + β0 )2 β0 = (F+ (ui+1,j ) − F+ (ui,j ))2 ,
1 3
α1 =
, ( + β1 )2 β1 = (F+ (ui,j ) − F+ (ui−1,j ))2 .
b − 1 , we utilize F− (ui,j ), F− (ui+1,j ), F− (ui+2,j ) and obtain To compute F i+ ,j
(4.9)
(4.10)
2
3 1 1 1 b− 1 = ω e0 ( F− (ui+1,j ) − F− (ui+2,j )) + ω e1 ( F− (ui,j ) + F− (ui+1,j )), F i+ 2 ,j 2 2 2 2 where α e0 α e1 , ω e1 = , ω e0 = α e0 + α e1 α e0 + α e1 α e0 =
1 3
( + βe0 )2
α e1 =
,
βe0 = (F− (ui+2,j ) − D− (ui+1,j ))2 ,
2 3
( + βe1 )2
,
βe1 = (F− (ui+1,j ) − F− (ui,j ))2 .
(4.11)
(4.12)
Similarly, we use numerical values G+ (ui,j−1 ), G+ (ui,j ) and G+ (ui,j+1 ) to compute b + 1 , and use G− (ui,j ), G− (ui,j+1 ) and G− (ui,j+2 ) to compute G b− 1. G i,j+ 2
i,j+ 2
9
5
Algorithm implementation
We now describe the implementation procedure of the aforementioned numerical methods that form a new algorithm for FSI computation. In general, we consider a solid structure represented by an elastic boundary immersed in an incompressible viscous fluid, which occupies a rectangular domain [xmin , xmax ] × [ymin , ymax ]. The fluid domain is partitioned by uniform Eulerian mesh xi,j = (xi , yj ), and uni,j denotes the numerical approximation of u(xi , yj , tn ), where 0 ≤ i ≤ Nx and 0 ≤ j ≤ Ny .
Let also the mesh size be h in both the x and y directions. Meanwhile, the elastic boundary is represented by Lagrangian coordinates, with X(ks , tn ) = (Xkns , Ykns ) at time tn , which is tracked from the initial position (Xk0s , Yk0s ), where 0 ≤ ks < Ns . The following steps advance the solution from t = tn to t = tn+1 .
Step 1: We first update the location of the Lagrangian points on the boundary by 1 half time-step; i.e., at t = tn+ 2 , we have n+1/2
Xks and
n+1/2
Yks
4t X 2 0≤q≤N
= Xkns +
= Ykns +
4t X 2 0≤q≤N
x
x
X
unqr δh (xq − Xkns , yr − Ykns )h2 ,
(5.1)
X
n vqr δh (xq − Xkns , yr − Ykns )h2 .
(5.2)
0≤r≤Ny
0≤r≤Ny
These correspond to the discretization of equations (2.6) and (2.7) by a half time-step. Step 2: We calculate the elastic force n+ 1 Fks 2
n+1/2
=
n+1/2
n+1/2
n+1/2
Tks +1/2 τks +1/2 − Tks −1/2 τks −1/2 4sks
n+1/2
n+1/2
(5.3)
n+1/2
at each Lagrangian point Xks = (Xks , Yks ), based on equation (2.2). Here n+1/2 n+1/2 we assume that each pair of adjacent Lagrangian points Xks and Xks +1 form a n+1/2
linar elastic fiber. Using Hooke’s law, the fiber tension Tks +1/2 is computed by n+1/2
n+1/2 Tks +1/2
= T0
n+1/2
|Xks +1 − Xks 4sks +1/2
|
!
−1 ,
(5.4)
where T0 is the elastic (or, spring) constant, and n+1/2 τks +1/2
=
n+1/2
n+1/2
n+1/2
n+1/2
Xks +1 − Xks
|Xks +1 − Xks
n+1/2
n+1/2
|
(5.5)
is the unit tangent of the fiber. Tks −1/2 and τks −1/2 are calculated in a similar way. 10
Step 3: We distribute the elastic force from the Lagrangian points to the Eulerian grids by X n+1/2 n+1/2 fi,j = Fks δh (xi − Xkns , yj − Ykns )∆sks (5.6) 0≤ks
based on equation (2.3). Step 4: We use the Hodge decomposition and WENO reconstruction to solve the incompressible Navier-Stokes equations (2.1) on the Eulerian mesh. A second- or third-order TVD Runge-Kutta method is employed for the temporal discretization (see equations 3.14 and 3.15) to generate the fluid velocity un+1 at t = tn+1 . Step 5: We then update the location of the elastic boundary at t = tn+1 by Xkn+1 = Xkns + 4t s
X
X
n+1/2
uij
0≤i≤Nx 0≤j≤Ny
n+1/2
δh (xi − Xks
n+1/2
, y j − Yk s
)h2
(5.7)
)h2 .
(5.8)
and Ykn+1 = Ykns + 4t s
6
X
X
n+1/2
vij
0≤i≤Nx 0≤j≤Ny
n+1/2
δh (xi − Xks
n+1/2
, yj − Yks
Numerical examples
Our proposed algorithm involves several components based on different numerical methods. In what follows, we provide a number of simulation examples to thoroughly test these numerical components, as well as the entire algorithm. Example 6.1. We first test the accuracy of our Hodge projection described in Section 3. We assume that the vector a(x, y) in (3.6) takes the form of sin2 x cos y + 2 sin (x + y) cos (x + y) cos y a= −2 sin x cos x sin y + 2 sin (x + y) cos (x + y) cos y − sin2 (x + y) sin y on the domain [0, 2π] × [0, 2π]. It is easy to verify that the exact solution for the
vector b is given by
sin2 x cos y b= . −2 sin x cos x sin y
We numerically solve b through (3.7)-(3.12), and the errors of b1 and b2 in L∞ norm (among all the grid points) are presented in Table 6.1, where N = Nx = Ny . We clearly observe that uniform fourth order is attained, consistent with the fourthorder finite difference method involved in the discretization. Similar results in L1 and L2 norms are also obtained (not shown here). 11
N L∞ error of b1 64 × 64 8.26E-5 128 × 128 5.21E-6 256 × 256 3.26E-7 512 × 512 2.04E-8 1024 × 1024 1.28E-9 2048 × 2048 7.97E-11
Order L∞ error of b2 6.04E-5 3.99 3.81E-6 4.00 2.39E-7 4.00 1.49E-8 4.00 9.34E-10 4.00 5.83E-11
Order 3.99 4.00 4.00 4.00 4.00
Table 6.1: Accuracy test for the Hodge decomposition. Example 6.2. Next, we verify the accuracy of our WENO reconstruction when solving the following compressible Navier-Stokes equations ( ρ(ut + F(u)x + G(u)y ) = µ4u + f , u(x, y, 0) = u0 (x, y)
(6.1)
on the computational domain [0, 2π] × [0, 2π] with periodic boundary conditions. The initial velocity is taken as
sin2 (x + 2y) cos(3y) u0 (x, y) = . cos3 (2x − 3y) sin (x + y)
The density is set as ρ = 1 and the viscosity is µ = 0.05. The exact solution to (6.1) can be found as sin2 (x + 2y) cos(3y)e−2µt u(x, y, t) = . 2 cos3 (2x − 3y + t) sin (x + y)e−3µt Here we apply the third-order WENO scheme for the convection term and fourthorder central difference scheme for the diffusion term. The temporal discretization is based on the third-order TVD Runge-Kutta method. Numerical results with the third-order WENO reconstruction at the end time T = 2 are given in Table 6.2, for both components of the velocity, u and v, in L∞ norm (among all the grid points). It is clear to see that, as the mesh size if refined, third-order accuracy is attained. Similar patterns in L1 and L2 norms are also observed (not presented here). Example 6.3. We proceed to check the accuracy of our combined Hodge decomposition and WENO method, using the Taylor-Green vortex test problem based on the incompressible Navier-Stokes equations (3.1). On the domain [0, 2π] × [0, 2π], an exact solution to the problem is given by 2t − −2k Re , u(x, y, t) = − cos(kx) sin(ke)e 2k2 t (6.2) v(x, y, t) = sin(kx) cos(ky)e− Re , 2t 4k p(x, y, t) = − 1 (cos(kx) + cos(ky))e− Re , 4
12
N L∞ error of u 80 × 80 5.36E-2 160 × 160 9.02E-3 320 × 320 1.33E-3 640 × 640 1.63E-4 1280 × 1280 1.54E-5
Order L∞ error of v 6.18E-2 2.57 8.57E-3 2.76 1.12E-3 3.04 1.37E-4 3.39 1.48E-5
Order 2.85 2.93 3.04 3.20
Table 6.2: Accuracy test for the WENO reconstruction. where k = 4 and µ = 0.05. As before, we use the third-order TVD Runge-Kutta method for the time marching, the third-order WENO scheme for the convective term, and the fourth-order central difference scheme for the viscous term. Numerical results at T = 2 are presented in Table 6.3 for both components of the velocity in L∞ norm. Again we observe that as the mesh size is refined, third-order accuracy is achieved, consistent with the third-order WENO scheme. Similar patterns hold for the errors in L1 and L2 norms. N L∞ error of u 80 × 80 3.80E-3 160 × 160 6.24E-4 320 × 320 7.90E-5 640 × 640 9.07E-6 1280 × 1280 9.49E-7
Order L∞ error of v 3.80E-3 2.60 6.24E-4 2.98 7.90E-5 3.12 9.07E-6 3.26 9.49E-7
Order 2.60 2.98 3.12 3.26
Table 6.3: Accuracy test for the combined Hodge-WENO method. Example 6.4. We further test the performance of our combined Hodge-WENO method when the solution is not smooth. We consider the double shear layer problem with the governing equation (3.1) and computational domain [0, 2π] × [0, 2π]. The initial condition is
u(x, y, 0) =
(
tanh ((y − π/2)/ρ),
y≤π
tanh ((3π/2 − y)/ρ),
y>π
and v(x, y, 0) = δ sin x. The boundary condition is periodic. The parameters are taken as µ = 0, ρ = π/15, δ = 0.05. Numerical approximations for the vorticity contours at T = 8 are plotted in Figure (6.1) with grid points of 128 × 128 and 512 × 512, respectively. The results agree well with those in literature (e.g., [32]). 13
Figure 6.1: Vorticity contours in Example 6.4. Left: N = 128. Right: N = 512. Example 6.5. Now we proceed to test our entire IB-WENO algorithm that integrates the IB method, the Hodge decomposition, and the WENO scheme, using the canonical problem of a stretched rubber band relaxing in an incompressible viscous fluid (see Figure 6.2). At its resting state, the rubber band is a circle with a radius r0 = 0.15. Initially, the rubber band is stretched to an ellipse, centered at (Ox , Oy ) = (0.5, 0.5), whose major semi-axis and minor semi-axis are a = 0.4 and b = 0.15, respectively. Both the fluid and the elastic boundary (i.e., the rubber band) are stationary at time t = 0. Once released, the boundary is set into motion due to the restoring elastic force and, since the fluid is incompressible, the boundary converges to a final circle which minimizes the energy and which has the same area as that of the initial ellipse. The √ radius of the final steady-state circle is re = ab ≈ 0.245. We set the fluid density ρ = 1 and the viscosity µ = 0.01. In our numerical simulation, the initial Lagrangian points on the boundary are prescribed as
2πks 0 , Xks = Ox + a cos Ns (6.3) 2πks 0 Yks = Oy + b sin , Ns where 0 ≤ ks ≤ Ns − 1 and Ns is the number of Lagrangian points. To reduce the
possible leakage of the volume (i.e., the area in 2D case), we impose the restriction that |X0ks − X0ks −1 | < 14
h , 2
where h = ∆x = ∆y is the fluid mesh size.
Figure 6.2: The stretched rubber band problem in Example 6.5. The solid, dashed, and dash-dot lines represent the initial, final, and resting positions of the rubber band, respectively. We have applied both the standard IB method (implemented with FFT) and the proposed IB-WENO algorithm to this problem. We find that both algorithms achieve first-order accuracy. This is expected, since the errors at the moving interface (where discontinuity occurs) dominate for the computation of the entire FSI problem, and the IB treatment of the interface reduces the overall accuracy to first order. In order to better compare the two algorithms, we look into the location and shape of the interface profile, as well as the area enclosed, throughout the process. Due to numerical errors, the computed interface does not approach a perfect, final circle. Thus, to quantify such errors, we calculate the minimum and maximum distances between the center (Ox , Oy ) and all the computed boundary points, at each time step t = tn , by q n rmin = min (Xkns − Ox )2 + (Ykns − Oy )2 , 0≤k≤Ns −1 q (6.4) n n n 2 2 rmax = max (Xks − Ox ) + (Yks − Oy ) , 0≤k≤Ns −1
and compare with the exact value re ≈ 0.245. Figure 6.3 shows the results by using the standard IB method and the IB-WENO method, where we observe that the values
of rmin and rmax generated by the standard IB method appear to diverge from the exact solution re over time, which leads to the “volume leakage” [20,28]. In contrast, 15
the results produced by the IB-WENO algorithm have better accuracy and stability over time.
Figure 6.3: Values of rmin and rmax versus time. Solid lines (red and green) represent the results generated by the standard IB method, and dashed lines (blue and black) represent the results generated by the IB-WENO method. The dash-dot line (pink) represents the exact solution. In addition, if we regard the computed interface as an “approximate circle” at each time tn , we may calculate the radius of this approximate circle by the average distance between the center (Ox , Oy ) and all the computed boundary points, r¯n =
Ns −1 q 1 X (Xkns − Ox )2 + (Ykns − Oy )2 . Ns k =0
(6.5)
s
Consequently, the area enclosed by the computed interface at t = tn can be approximately determined as S n = π(¯ rn )2 , and we compare this value with the exact area Se = πab (which is the area enclosed by the initial ellipse). Figure 6.4 displays the results by the two methods: the standard IB, and the IB-WENO. We observe that the standard IB method generates smaller error initially, but the error keeps increasing over time (another indication of the volume leakage). In contrast, the IB-WENO 16
method seems to yield more stable and consistent result, with better accuracy eventually.
Figure 6.4: |S n − Se | versus time. The (red) solid line represents the result of the IB-WENO method, and the (green) dashed line represents the result of the standard IB method.
7
Discussion
We have presented a new approach that integrates the IB method, the Hodge decomposition, and a finite-difference WENO scheme for FSI computation. With this hybrid algorithm, named IB-WENO, we maintain the simple treatment of the fluidstructure interface, while achieve a high-order accuracy for regions near and away from the interface where the solution is sufficiently smooth. This work represents a pilot study in combining the IB and WENO methods. The combined algorithm is easy to implement, applicable to both periodic and non-periodic domains, and our numerical results demonstrate the performance of the algorithm and its various components. In particular, its overall accuracy is comparable to or better than that of the standard combination of the IB method and FFT, when implemented in periodic 17
space dimensions. Additionally, our results in the stretched rubber band test problem shows a better volume conservation than that of the standard IB method. In a recent study [36], the authors proposed a method that combines a WENO scheme and the ghost-cell method (a variation of the original IB method) to compute shock-structure interaction, where the fluid flow is described by the Euler equations; i.e., fluid viscosity is neglected. The strength of their method is that it may handle high-speed flow with strong discontinuities. The formulation and implementation of their method, however, are much more complex than ours. Currently, our integrated algorithm achieves first-order accuracy at the interface because of the IB treatment, and third-order accuracy in other parts of the domain because of the third-order WENO scheme employed. There is an apparent mismatch of the orders of accuracy for these two components in the hybrid algorithm. On one hand, the degradation of accuracy at the interface is inevitable due to the jump discontinuity of the solution. The majority of the immersed type methods achieve first-order accuracy, with only a few exceptions (such as the immersed interface method [20,21]) that require additional, highly nontrivial procedures to boost the accuracy at the interface. On the other hand, a different direction to alleviate the order “mismatch” and improve the accuracy at the interface is to employ a fine mesh in the small region that contains the interface, and use a coarse mesh for the remaining parts of the domain. In other words, the low-order method for the interface treatment (such as our IB method) is based on smaller grid size, whereas the high-order method for the majority of the domain (such as our WENO scheme) is based on larger grid size, so as to reach a balance for the accuracy (not the order of accuracy, but rather the magnitude of the numerical error) between the two methods. Consequently, a local mesh refinement procedure [31] would be a natural choice to achieve this goal. We are currently working along this line and implementing the immersed boundary method with adaptive mesh refinement.
Acknowledgments Y.L. acknowledges partial support from the Simons Foundation (No. 426993). M.Z. acknowledges partial support from the National Natural Science Foundation of China (No. 11471305). The work of J.W. was partially supported by the National Science Foundation (Nos. 1319078 and 1520672), and by the Tennessee Higher Education Commission through the CEACSE program. The authors thank the editor for handling this paper and the reviewer for helpful comments. 18
References [1] Marcos Castro, Bruno Costa, and Wai Sun Don. High order weighted essentially non-oscillatory WENO-Z schemes for hyperbolic conservation laws. Journal of Computational Physics, 230(5):1766–1792, 2011. [2] Bernardo Cockburn and Chi-Wang Shu. TVB runge-kutta local projection discontinuous Galerkin finite element method for conservation laws. ii. general framework. Mathematics of computation, 52(186):411–435, 1989. [3] Bernardo Cockburn and Chi-Wang Shu. The runge–kutta discontinuous galerkin method for conservation laws v: multidimensional systems. Journal of Computational Physics, 141(2):199–224, 1998. [4] Earl H Dowell and Kenneth C Hall. Modeling of fluid-structure interaction. Annual Review of Fluid Mechanics, 33(1):445–490, 2001. [5] Travis C Fisher, Mark H Carpenter, Nail K Yamaleev, and Steven H Frankel. Boundary closures for fourth-order energy stable weighted essentially non-oscillatory finite-difference schemes. Journal of Computational Physics, 230(10):3727–3752, 2011. [6] Robert D Guy and David A Hartenstine. On the accuracy of direct forcing immersed boundary methods with projection methods. Journal of Computational Physics, 229(7):2479–2496, 2010. [7] Ami Harten. High resolution schemes for hyperbolic conservation laws. Journal of computational physics, 49(3):357–393, 1983. [8] Ami Harten. Preliminary results on the extension of ENO schemes to twodimensional problems. Nonlinear Hyperbolic Problems, C. Carasso, P.-A. Raviart, and D. Serre, eds., LNM, 1270:23–40, 1987. [9] Ami Harten. ENO schemes with subcell resolution. Journal of Computational Physics, 83(1):148–184, 1989. [10] Ami Harten, Bjorn Engquist, Stanley Osher, and Sukumar R Chakravarthy. Uniformly high order accurate essentially non-oscillatory schemes, III. Journal of computational physics, 71(2):231–303, 1987. [11] Ami Harten and Stanley Osher. Uniformly high-order accurate nonoscillatory schemes. I. SIAM Journal on Numerical Analysis, 24(2):279–309, 1987. 19
[12] Ami Harten, Stanley Osher, Bj¨orn Engquist, and Sukumar R Chakravarthy. Some results on uniformly high-order accurate essentially nonoscillatory schemes. Applied Numerical Mathematics, 2(3-5):347–377, 1986. [13] Andrew K Henrick, Tariq D Aslam, and Joseph M Powers. Mapped weighted essentially non-oscillatory schemes: achieving optimal order near critical points. Journal of Computational Physics, 207(2):542–567, 2005. [14] Jan S Hesthaven, Sigal Gottlieb, and David Gottlieb. Spectral methods for time-dependent problems. Cambridge University Press, 21, 2007. [15] Gene Hou, Jin Wang, and Anita Layton. Numerical methods for fluid-structure interaction—a review. Communications in Computational Physics, 12(2):337– 377, 2012. [16] Changqing Hu and Chi-Wang Shu. Weighted essentially non-oscillatory schemes on triangular meshes. Journal of Computational Physics, 150(1):97–127, 1999. [17] Bj¨orn H¨ ubner, Elmar Walhorn, and Dieter Dinkler. A monolithic approach to fluid–structure interaction using space–time finite elements. Computer methods in applied mechanics and engineering, 193(23):2087–2104, 2004. [18] Guang-Shan Jiang and Chi-Wang Shu. Efficient implementation of weighted ENO schemes. Journal of computational physics, 126(1):202–228, 1996. [19] David A Kopriva. A practical assessment of spectral accuracy for hyperbolic problems with discontinuities. Journal of Scientific Computing, 2(3):249–262, 1987. [20] Randall J LeVeque and Zhilin Li. Immersed interface methods for stokes flow with elastic boundaries or surface tension. SIAM Journal on Scientific Computing, 18(3):709–735, 1997. [21] Zhilin Li. An overview of the immersed interface method and its applications. Taiwanese journal of mathematics, pages 1–49, 2003. [22] Wing Kam Liu and Shaoqiang Tang. Mathematical foundations of the immersed finite element method. Computational Mechanics, 39(3):211–222, 2007. [23] Xu-Dong Liu, Stanley Osher, and Tony Chan.
Weighted essentially non-
oscillatory schemes. Journal of computational physics, 115(1):200–212, 1994. 20
[24] Andrew Majda, James McDonough, and Stanley Osher. The Fourier method for nonsmooth initial data. Mathematics of Computation, 32(144):1041–1081, 1978. [25] J Mohd-Yusof. Combined immersed boundary/b-spline methods for simulations of flow in complex geometries. Annual Research Briefs, Center for Turbulence Research, Stanford University, pages 317–327, 1999. [26] Mohammad Motamed, Colin B Macdonald, and Steven J Ruuth. On the linear stability of the fifth-order weno discretization. Journal of Scientific Computing, 47(2):127–149, 2011. [27] Charles S Peskin. The immersed boundary method. Acta numerica, 11:479–517, 2002. [28] Charles S Peskin and Beth Feller Printz. Improved volume conservation in the computation of flows with immersed elastic boundaries. Journal of computational physics, 105(1):33–46, 1993. [29] PB Ryzhakov, Riccardo Rossi, SR Idelsohn, and E O˜ nate. A monolithic lagrangian approach for fluid–structure interaction problems. Computational mechanics, 46(6):883–899, 2010. [30] Susana Serna and Antonio Marquina.
Power ENO methods: a fifth-order
accurate weighted power ENO method.
Journal of Computational Physics,
194(2):632–658, 2004. [31] Chaopeng Shen, Jing-Mei Qiu, and Andrew Christlieb. Adaptive mesh refinement based on high order finite difference weno scheme for multi-scale simulations. Journal of Computational Physics, 230(10):3780–3802, 2011. [32] Chi-Wang Shu.
Essentially non-oscillatory and weighted essentially non-
oscillatory schemes for hyperbolic conservation laws. Advanced numerical approximation of nonlinear hyperbolic equations, pages 325–432, 1998. [33] Chi-Wang Shu and Stanley Osher.
Efficient implementation of essentially
non-oscillatory shock-capturing schemes. Journal of Computational Physics, 77(2):439–471, 1988. [34] J Steindorf and HG Matthies. Numerical efficiency of different partitioned methods for fluid-structure interaction. ZAMM-Journal of Applied Mathematics and
21
Mechanics/Zeitschrift f¨ ur Angewandte Mathematik und Mechanik, 80(S2):557– 558, 2000. [35] Jan Vierendeels, Kris Dumont, and PR Verdonck. A partitioned strongly coupled fluid-structure interaction method to model heart valve dynamics. Journal of Computational and Applied Mathematics, 215(2):602–609, 2008. [36] Min Xu, Tao Yang, and Mingjun Wei. Implementation of immersed boundary method in weno scheme to simulate shock-structure interaction. ASME 2017 Fluids Engineering Division Summer Meeting, pages V01BT11A013– V01BT11A013, 2017. [37] LT Zhang and M Gay. Immersed finite element method for fluid-structure interactions. Journal of Fluids and Structures, 23(6):839–857, 2007.
22