A Multirate Approach for Fluid-Structure Interaction Computation with Decoupled Methods
AA Multirate Approach for Fluid-Structure Interaction Computationwith Decoupled Methods
Lian Zhang , Mingchao Cai Mo Mu , Abstract
We investigate a multirate time step approach applied to decoupled methods in fluid andstructure interaction(FSI) computation, where two different time steps are used for fluid andstructure respectively. For illustration, the multirate technique is tested by the decoupled β scheme.Numerical experiments show that the proposed approach is stable and retains the same orderaccuracy as the original single time step schemes, while with much less computational expense. Key words:
Fluid and structure interaction, decoupled methods, multirate time step, stability, β scheme.
1. Introduction
Fluid structure interaction (FSI) problems are extremely important because they appear inmany scientific and engineering applications [3, 4, 5, 7, 14, 21, 20]. In the literature, for solvingFSI problems, both fully implicit and decoupled approaches have been applied. The fully implicitdiscretization approach leads to coupled schemes [22], in which the equations of fluid dynamics,structural mechanics, and mesh moving are solved simultaneously in a fully coupled fashion. Althoughthe coupled schemes are unconditionally stable, they result in significant difficulties and inflexibilityin the design and choice of mesh generation, PDE discretization, algebraic solvers, as well as softwaredevelopment. In recent years, some decoupled approaches, called loosely coupled or partitioned, orexplicit coupling approaches, have been developed [1, 2, 13, 17]. In these decoupled approaches,existing fluid and structure solvers are used, the equations of fluid dynamics, structural mechanics,and mesh moving are solved sequentially or independently. However, the stability and convergencecould not be guaranteed if the decoupling technique is not well-designed. For instance, in an explicitdecoupling algorithm based on the Dirichlet-Neumann splitting, one solves the fluid dynamicsequations with the velocity Dirichlet boundary conditions imposed by using the extrapolated value ofstructure velocity at the interface, then solves the structural mechanics equations with the Neumannboundary condition provided by the updated fluid interface traction, and then updates the solutionof the mesh moving by using the newest structural displacements at the interface. However, thisexplicit Dirichlet-Neumann scheme is known to be unconditionally unstable due to the so-calledartificial added-mass effect [10, 12]. Nevertheless, note that the fluid and solid possess quite different E-mail address: [email protected]. Department of Mathematics, the HongKong University of Science andTechnology. E-mail address: [email protected]. Department of Mathematics, Morgan State University. E-mail address: [email protected]. Department of Mathematics, the HongKong University of Science and Technology.
Preprint submitted to Elsevier October 30, 2018 a r X i v : . [ m a t h . NA ] O c t hysical properties, such as stiffness and velocity. It is natural to treat different models in their ownphysical regions differently for various numerical considerations. Therefore, decoupled approachesare more favorable, not only for FSI problems, but also for other coupled multi-domain, multi-physicsapplications [8, 15, 16, 18, 19].In this work, we are interested in the coupling of an incompressible viscous fluid flow modelwith a thin-walled structure model. In recent years, there are two notable works: the so-calledRobin-Neumann scheme and the β scheme. In these two schemes, the interface coupling conditionsare treated and approximated as a Robin type condition (a linear combination of Dirichlet conditionand Neumann condition). In time marching of the Robin-Neumann scheme, one firstly solves thefluid model with the Robin interface condition approximated by using the data from the solid regionat the previous time step or by certain extrapolation strategies, and then solves the structure modelwith the Neumann interface condition on the interface supplemented by the latest data computedfrom the fluid region. While, in the so-called β scheme, the authors of [6] split the thin-walledstructure equation into two parts with a constant β . One of the two parts gives a Robin-typeinterface condition which is used in the fluid step while the other part can be treated as the structureequation with Neumann condition. Different from the Robin-Neumann scheme, the structure modelwith Neumann interface condition is solved firstly in the β scheme. We note that, when β = 1, thedifference between the β scheme and Robin-Neumann scheme exists only in that which model issolved firstly. They actually apply the same strategy for handling the interface conditions. Afterserous investigation and comparison, we show that the performance of the Robin-Neumann schemeis quite similar to that of β scheme (cf. Section 4).In this work, our main interest is to extend the β scheme to a multirate time-stepping algorithm.By multi-rate timestepping, we mean that different time step sizes are used in different subdomains.Such a multirate time-stepping strategy is in accordance with the physical laws because the FSIproblems are multi-scale problems in time. Particularly, for the Stokes flows coupled with thin-walledstructures, the variables in the structure subdomain vary much more rapidly than those variables inthe fluid subdomain. In the literature, a multirate time step technique was introduced in [18, 19] forcoupled fluid-porous media flow models. The whole time interval [0 , T ] is partitioned into certaincoarse time grids with the time-step size τ coarse . Within each coarse time step, the free fluid flowsolutions are computed for multiple fine time steps with the boundary information at the interfacesupplemented by the porous medium region (using the previous time step data). When it reachesthe end of current coarse time grid, the porous medium solutions are updated by using the data fromthe fluid solutions. Such a multirate method is proved to be stable and convergent with the ordersof accuracy in space and time depending on the spatial discretization order and time discretizationorder. In this work, for the Stokes flows coupled with the thin-walled structures, we choose a finertime step size for the structure model while applies a coarse time step size for the fluid flow model.Although, in a multirate time-stepping approach, one can freely choose a fine time step for eithermodel, our numerical tests show that the current choice leads to a better numerical performance.The paper is organized as follows. In Section 2, we describe a FSI model for coupling a Stokesflow with a thin-walled structure. In Section 3, a multirate β scheme is outlined for the coupled FSImodel. Numerical experiments are presented in Section 4 to show the stability and convergence ofour scheme. Conclusions are given in Section 5. 2 . A Stokes Flow Interacting with A Thin-Walled Structure In this section, we describe the model problem studied in [6, 11]. In the coupled FSI model,the fluid flow motion is governed by the Stokes equations in a d -dimensional ( d = 2 ,
3) domainΩ f and the structure is assumed to be a linear thin-solid defined on a ( d − − manifold Γ. Theboundary ∂ Ω f = Γ ∪ Γ D ∪ Γ N with Γ D and Γ N representing the boundaries imposed with Dirichletand Neumann conditions respectively. The coupled model problem reads as: finding the fluidvelocity u f : Ω f × R + → R d , the fluid pressure p f : Ω f × R + → R , and the solid displacement d : Γ × R + → R d − such that ρ f ∂ t u f − div σ f ( u f , p f ) = 0 in Ω f , div u f = 0 in Ω f , u f = 0 on Γ D , σ f ( u f , p f ) n = f N on Γ N , (1)and u f = u s = ˙ d on Γ ,ρ s (cid:15)∂ t ˙ d + L e d + L v ˙ d = − σ f ( u f , p f ) n on Γ , d = on ∂ Γ , (2)satisfying the initial conditions u f (0) = u f , d (0) = d . Here, ρ f and ρ s are the fluid density and the solid density respectively, (cid:15) is the solid thickness, ˙ d isthe solid velocity, n is the exterior unit normal vector to ∂ Ω f , ε ( u f ) = 12 ( ∇ u f + ∇ u Tf ) , σ f ( u f , p f ) = − p f I + 2 µε ( u f )with µ being the fluid dynamic viscosity, f N is a given surface force on Γ N , L e and L v stand forthe elastic and viscous contributions respectively. Here and hereafter, we use L s = L e + L v torepresent the solid tensor. In the coupled model, two interface conditions are enforced: the Dirichletcondition (2) guarantees the continuity between the fluid velocity and the structure velocity at Γ;the Neumann condition (2) ensures the continuity of the stresses at Γ. We comment here that theequation (2) is not only an interface coupling condition but also the structure governing equationin the coupled model.Let V f and V s be the H local spaces, associated with the appropriate global Dirichlet conditions,for the fluid and structure regions, respectively. Let Q f be the L pressure space for Stokes model. V f = H (Ω f ) = (cid:8) u f ∈ ( H (Ω f )) d | u f = 0 on Γ D (cid:9) , V s = H (Ω s ) = (cid:8) u s ∈ ( H (Γ)) d − | u s = 0 on ∂ Γ (cid:9) ,Q f = L (Ω f ) . Defining V ≡ { ( u f , u s ) ∈ V f × V s | u f | Γ = u s | Γ } , then the implicit or monolithic weak problemfor the coupled model reads as: finding ( u , p f ) ∈ V × Q f , and d ∈ V s , such that u s = ∂ d ∂t and (cid:40) ( δ t u , v ) + a Ω ( u ; u , d , v ) + b ( v , p f ) = f ( v ) , ∀ v ∈ V ,b ( u , q ) = 0 , ∀ q ∈ Q f , (3)3here u ≡ ( u f , u s ), v ≡ ( v f , v s ), δ t u ≡ ( ρ f ∂ u f ∂t , ρ s ∂ u s ∂t ), a Ω ( u ; u , d , v ) ≡ a Ω f ( u f ; u f , v f ) + a Ω s ( u s ; u s , d , v s ) represent the stress tensor parts, b ( u , q ) = − (cid:82) Ω f q div u f . Note that the Neumanninterface condition is automatically guaranteed in the weak form, while the Dirichlet interfacecondition is enforced in the definition of V , which reflects the coupling.
3. Numerical Algorithms β Scheme
Algorithm 1:
The β scheme.For m = 0 , , , ...N − u m +1 s such that (cid:40) ρ s ε ˜ u m +1 s − u ms ∆ t + L s ( d m +1 ) = − β σ f ( u mf , p mf ) n , on Γ ,d t d m +1 = ˜ u m +1 s , on Γ . (4)2. Fluid step: find u m +1 f , p m +1 f and u m +1 s such that ρ f ∆ t ( u m +1 f − u mf ) − div σ f ( u m +1 f , p m +1 f ) = , in Ω f , div u m +1 f = 0 , in Ω f ,ρ s ε u m +1 s − ˜ u m +1 s ∆ t = − σ f ( u m +1 f , p m +1 f ) n + β σ f ( u mf , p mf ) n , on Γ , u m +1 f = u m +1 s , on Γ . (5)In Algorithm 1, we describe the β scheme proposed in [6]. The key of the β scheme is that thestructure equation is split as ρ s (cid:15) u m +1 s − ˜ u m +1 s (cid:124) (cid:123)(cid:122) (cid:125) + ˜ u m +1 s − u ms ∆ t + L s ( d m +1 , ˙ d m +1 ) = − σ f ( u m +1 f , p m +1 f ) n + β σ ( u mf , p nf ) n (cid:124) (cid:123)(cid:122) (cid:125) − β σ f ( u nf , p nf ) n . (6)Here, the ” (cid:124)(cid:123)(cid:122)(cid:125) ” parts are used in the fluid step as a Robin-type interface condition whereas theother parts are computed in the structure step. 4 .2. A Multirate β Scheme
Figure 1: An illustration of a multirate time stepping technique.
Algorithm 2:
A multirate β scheme.For k = 0 , , , ...N −
1, set m k = r · k ,1. Structure steps: for m = m k , m k + 1 , m k + 2 , ..., m k +1 − (cid:40) ρ s ε ˜ u m +1 s − u ms ∆ t s + L s ( d m +1 ) = − β σ f ( u m k f , p m k f ) n , on Γ ,d t s d m +1 = ˜ u m +1 s , on Γ . (7)2. Fluid step: ρ f ∆ t f ( u m k +1 f − u m k f ) − div σ f ( u m k +1 f , p m k +1 f ) = , in Ω f , div u m k +1 f = 0 , in Ω f ,ρ s ε u mk +1 s − ˜ u mk +1 s ∆ t f = − σ f ( u m k +1 f , p m k +1 f ) n + β σ f ( u m k f , p m k f ) n , on Γ , u m k +1 f = u m k +1 s , on Γ . (8)In the β scheme, the coupled FSI system is split into fluid and structure steps in a sequentialmanner. It allows us to solve the fluid model and the structure model separately. However, in thealgorithm, both the fluid solver and the structure solver use the same time step size. We note thatthe time scale in the structure part maybe different from the time scale in the fluid part. It is notnecessary to use the same time step size in both steps. Thus, we apply a multirate time steppingtechnique to the β scheme. Intuitively, there are two possible choices of the time-stepping technique:one is to use a bigger time step size for the fluid solver, the other is using a bigger time step sizefor the structure solver. Based on our numerical observations, the algorithm which uses a bigger5ime step size for the fluid solver whereas applies a smaller time step size for the structure solver (cf.Figure 1) gives a better accuracy. We therefore describe our multirate β scheme in Algorithm 2 andthe corresponding fully discrete weak form is given in Algorithm 3. Algorithm 3:
The fully discrete weak form for the multirate β scheme.For k = 0 , , , ...N −
1, set m k = r · k
1. Structure step: for m = m k , m k + 1 , m k + 2 , ..., m k +1 − , find ˜ u m +1 sh ∈ V sh with d t s d m +1 h = ˜ u m +1 sh such that ∀ v sh ∈ V sh , there holds ρ s ε (cid:18) ˜ u m +1 sh − u msh ∆ t s , v sh (cid:19) Γ + a s ( d m +1 h , v sh ) = − β (cid:16) σ f ( u m k fh , p m k fh ) n , v sh (cid:17) Γ . (9)2. Fluid step: find ( u m k +1 fh , u m k +1 sh , p m k +1 fh ) ∈ ( V fh , V sh , Q fh ) with u m k +1 fh | Γ = u m k +1 sh such that ∀ ( v fh , v sh , q fh ) ∈ ( V fh , V sh , Q fh ) with v fh | Γ = v sh , there holds ρ f (cid:32) u m k +1 fh − u m k fh ∆ t f , v fh (cid:33) Ω + a f ( u m k +1 fh , v fh ) − b ( p m k +1 fh , v fh ) + b ( q fh , u m k +1 fh )+ ρ s ε (cid:18) u m k +1 sh − ˜ u m k +1 sh ∆ t f , v sh (cid:19) Γ = β (cid:16) σ f ( u m k fh , p m k fh ) n , v sh (cid:17) Γ . (10) Remark 1. If m = m k , we have u ms = u m k s from the fluid step. If m > m k , we take u ms = ˜ u ms from the structure step.
4. Numerical Experiments
In this section, we present numerical experiments to demonstrate the convergence and stabilityperformance of the multirate β scheme. The benchmark test is for numerically solving a 2D pressurewave interacting with a thin-walled structure. The displacements of the interface are assumed tobe infinitesimal and that the Reynolds number in the fluid is assumed to be small. The 2D fluiddomain is a rectangle Ω f = [0 , L ] × [0 , R ] with L = 6cm and R = 0 . , L ] × R . See Figure 2 for the geometry configuration. Figure 2: Geometrical configuration
The physical parameters are: ρ f = 1 . ρ s = 1 .
1, and µ = 0 . L s ( d , ˙ d ) = c ∂ x d + c d , c = E(cid:15) ν ) , c = E(cid:15)R (1 − ν ) with (cid:15) = 0 .
1, the Poisson ratio ν = 0 .
5, and theYoung modulus E = 0 . · . During T ∗ = 5 · − seconds, a pressure-wave, P ( t ) = P max (1 − cos(2 tπ/T ∗ )) / P max = 2 · ,
6s prescribed on the fluid inlet boundary, a zero traction is enforced on the fluid outlet boundary,a no-slip condition is imposed on the lower boundary y = 0. For the solid, we fix the twoendpoints by imposing d = at x = 0 and x = 6. In all the following tests, similar to [11], wegenerate a reference solution using the fully implicit scheme with a high space-time grid resolution( h = 3 . × − , ∆ t = 10 − ). Figure 3: Comparisons of the numerical results obtained by the implicit scheme, the RN scheme and the β schemeunder the setting: h = 0 .
05 and ∆ t s = 10 − . In Figure 3, we compare the numerical results obtained by using the implicit scheme, theRobin-Neumann scheme and the β scheme. The mesh size and the time step size setting for theRN and the β scheme is: h = 0 .
05 and ∆ t = 10 − . From the results, we observe that both theRobin-Neumann scheme and the β scheme give very good approximations to the solution obtainedby using the implicit scheme. Most importantly, we see clearly that the results obtained by using theRobin-Neumann scheme have little difference with those obtained by using the β scheme. Therefore,in the following, we only report the numerical results obtained by using the β scheme or the multirate β scheme.In addition, we test two different multirate strategies: choosing a bigger time step size for thefluid model or a bigger time step size for the structure model. The numerical results are presentedin Figure 4. From the figure, we see that if the time ratio r = 2, a bigger step size in the fluid modelwhile applying a smaller time size for the structure model gives more accurate numerical solutionthan that obtained by using the other strategy. Furthermore, from our experiments and experience,the multirate β scheme with a smaller time step size for the fluid part is unstable and the numericalresults will be messed up when r = 5 or 10.In order to test whether a large time step ratio will cause instability, we fix ∆ t s and h whilevary the time ratio r = 1 , , , ,
50. The numerical results for t = 0 .
015 are reported in Figure 5with the structure time step size ∆ t s = 10 − while the mesh size h = 0 . h = 0 .
01 (right).From the left part of the figure, we see that the solid displacement along interface obtained bythe multirate β scheme with r = 1 , , ,
10 are almost the same as that obtained by using theimplicit scheme. Moreover, the multirate β scheme with r = 20 ,
50 are still stable although the7 igure 4: Comparison of the β scheme and two different multirate β schemes with h = 0 . , ∆ t = 10 − and r = 2.(a) (b)Figure 5: Numerical displacements under the settings: h = 0 .
01 (left) h = 0 .
01 (right) and ∆ t s = 10 − . errors become lager because of the lager time step size for the fluid model. To further investigate thestability and the convergence of the multirate β scheme, we apply a finer mesh size h = 0 .
01 (whilekeeping ∆ t s = 10 − ). The numerical results are presented in the right part of Figure 5. From theresults, we have almost the same observations as those obtained under the setting h = 0 .
1. Therefore,the multirate β scheme is stable even the time size ratio is large. To have a good approximation,one only needs to keep the time step size ratio be not too large.In Figure 6, for t = 0 . , . , . β scheme, and the multirate β scheme. (We comment here that the multirate β scheme is nothing else but the β scheme when r = 1.) By comparing the results obtained by usingdifferent algorithms, we see that the numerical results obtained by using the multirate β scheme arevery good approximations to those obtained by using the implicit scheme.In order to examine the orders of convergence which are second order in h and first order in t ,we decrease the mesh size by a factor of two and the time step size by a factor of four at each level8 a) t = 0 .
005 (b) t = 0 .
010 (c) t = 0 . t = 0 .
005 (e) t = 0 .
010 (f) t = 0 . t = 0 .
005 (h) t = 0 .
010 (i) t = 0 . t = 0 . , . , .
015 obtained by the implicit scheme (top), the multirate β scheme with r = 1 (middle) and r = 10 (bottom) with h = 0 .
01 and ∆ t s = 0 . h = 0 .
1, ∆ t s = 0 . { h, t s } = { . · (0 . i , . · (0 . i } , i = 0 , , , , . (11)In Figure 7, we present the relative errors of the primary variables ( u f , p f and d ) at t = 0 . β scheme (with the step ratio being equal to r = 1 or 10). From thefigures, we see that the numerical error is decreased by a factor of four as the mesh size and thetime step size are refined once. (a) Relative error of u f (b) Relative error of p f (c) Relative error of d Figure 7: Relative error of primary variables with the spacing h and time t in (11).
Finally, in order to highlight the advantage of the multirate β scheme, we compare the CPUtimes for the different numerical algorithms under several settings of the time step sizes and meshsizes. We fix, ∆ t s = 10 − and vary h = , , , , . The CPU times of using differentnumerical algorithms are summarized in Table 1. From the table, we see that the multirate β schemetakes much less cost than that of the implicit scheme, in particular, when r is large. Therefore, weconclude that the multirate schemes improve the efficiency.10 able 1: CPU times (in seconds) for the implicit scheme and the multirate β scheme (with r = 1 or 10) under differentsettings of mesh sizes (∆ t s = 10 − is fixed). implicit Scheme multirate β scheme r=1 multirate β scheme r=10 h = h = h = h = h =
5. Conclusions
Fluid structure interaction problems appear in many engineering and science applications. Suchproblems are multi-domain, multi-physics problems with multiscales. In this paper, we developa multirate β scheme for solving the coupled model of Stokes flow interacting with a thin-walledstructure. We note that the incompressible fluid model and the thin-walled structure model possessdifferent time scales. It is natural to apply a multirate time strategy to solve such a model. First ofall, our algorithm is a decoupled algorithm which have many advantages as stated in the introductionpart. Moreover, extensive numerical experiments are presented to show that the proposed scheme isefficient and accurate. Compared with the coupled implicit scheme, our algorithm uses much lesscomputational cost to achieve the same order of accuracy.
6. Acknowledgement
The first and the third authors’ research is supported in part by Hong Kong RGC CompetitiveEarmarked Research Grant HKUST16301218 and NSFC (91530319,11772281). The second author’swork is supported in part by NSF Grant