Due to the many applications of FSI problems in our world, there is a strong need for numerical solvers to be quick and accurate. Adaptive time-stepping is especially crucial because it allows for highly complex problems to be more efficiently and accurately solved. In cases when the dynamics of a problem are unknown, choosing a proper time step is infeasible and hence, adaptive time-stepping is used to find an appropriate time step. At certain points in time in a problem's simulation, small time steps may be necessary in ensuring the accuracy of the solution. However, at other times, larger steps may be implemented, which could help to reduce computational time. In particular, we are interested in solving problems with a viscous, incompressible fluid and an elastic structure. Our mathematical model uses the Navier-Stokes equations in the ALE form to describe the fluid whereas the solid is modeled using elastodynamic equations in the Lagrangian framework. To numerically solve this set of equations, researchers may implement a monolithic solver [44,26,31,4,32,5,16,48,25,47,33], which solves the problems in one fully coupled system, or may choose a partitioned method [1-3,11,10,15,19-21, 29,40,43,45,46], which separately solves the fluid and structure subproblems. For biomedical applications, in particular, partitioned methods encounter the added mass effect, which has been associated with numerical instabilities [13], making the development of partitioned methods for FSI problems especially challenging.
Adaptive time-stepping techniques are seldom found in FSI literature. One of the main difficulties in adaptive time stepping is maintaining stability [56], since many widely used time integration methods lose their stability properties when they are used with a variable time step [12]. Obtaining unconditional stability can be particularly difficult when an adaptive time stepping method is used together with partitioning.
Adaptive time-stepping methods for monolithic FSI problems have been explored in [41,18]. A method proposed in [41] includes a posteriori error estimation. In computing the local discretization error, the authors compare the marching scheme with the auxiliary scheme of order p and p, respectively. If the local error is not below some user-given tolerance, an optimized scaling factor is used to adjust the time step size accordingly. Here, the method is numerically explored using a monolithic solver, but is also applicable to a partitioned method. In [18], the authors tested a dual-weighted residual estimator, showing more reliability than a heuristic time-step refinement, albeit with more rigorous implementation efforts. This is numerically verified using a monolithic solver. We also mention the work in [57], where an extended variable timestepping scheme based on the Adams-Bashforth-Moulton method is presented for obtaining a solution of the interaction between a fluid and a rigid body. This second-order method involves the use of one predictor and several correctors, thus requiring an iteration process during each time step. To reduce the number of subiterations, Aitken's dynamic underrelaxation scheme is used.
In this work, we investigate an adaptive, partitioned method for FSI problems. To discretize the problem in time, we use a refactorization of Cauchy's adaptive, one-legged θ -like method. Using this method, the problem is solved as a sequence of a backward Euler (BE) step and a forward Euler (FE) step, where the FE step is written as a linear extrapolation. When θ = 1 2 , the method corresponds to the midpoint rule [8], in which case it features second-order accuracy in time. This approach is combined with a partitioned strategy, where the fluid and structure subproblems are decoupled and solved separately using Robin boundary conditions at the interface. Hence, in the BE step, the fluid and structure subproblems are solved iteratively until the convergence is reached. Then, the variables are extrapolated in the FE step and the time step is adapted. This method, in case of fixed time-stepping, has been analyzed in [8], where it was shown that the subiterations used in BE step linearly converge, and that the method is unconditionally stable. Numerical examples show second-order convergence in time when θ = 1 2 . While the theoretical work done in [8] can easily be extended to adaptive time-stepping, the computational aspects of time-adaptivity, which have not been investigated, bring many questions and are the focus of this work.
To adapt the time step, we use an LTE estimator with two variations for calculating the error [52]. In the first variation, we use Milne's device [42,27,52,36,37,39], which uses two different schemes with the same oder of accuracy but different error constants. Namely, we subtract our second-order midpoint solution from a modified version of Adams-Bashforth twostep (AB2) method (modified as it is calculated at the midpoint) and multiply the result by an appropriate factor. In the second variation, the LTE is evaluated using Taylor expansions at t n+ 1 2 , t n-1 2 , and t n-3 2 . For each of the cases, if the LTE is above a tolerance, δ, the time-step is modified but the trial is rejected and solution is computed again using the modified time step. This process repeats until the LTE is smaller than δ, in which case the time step is modified and the solver moves onto the next time level. In our numerical section, both variations of LTEs are numerically investigated.
In addition to investigating different local truncation errors for the adaptive time-stepping, we numerically explore the effects of the combination parameter used in Robin boundary conditions at the interface. By testing different values of the combination parameter, we are able to investigate which cases yield the least amount of subiterations in the BE step of our algorithm. Furthermore, we compare the results obtained using the adaptive and fixed time-stepping. The optimal choices of adaptivity and other problem parameters obtained in the first numerical example are then tested on two common problems in FSI simulations: the Hron-Turek benchmark describing the flow in a two-dimensional channel which contains an elastic bar attached to a rigid cylinder, and a three-dimensional study of fluid flow in an elastic channel representing a simplified model of blood flow.
The remainder of this paper is as follows: Section 2 introduces the problem setting, with the numerical solver detailed out in Section 3. The computation of the LTE is discussed in Section 4. The numerical examples are presented in Section 5, and the conclusions are drawn in Section 6.
Suppose that a bounded, deformable domain (t) ⊂ R d (where d = 2, 3 is the spatial dimension and t ∈ [0, T ] denotes time) consists of two regions, F (t) and S (t), separated by a common interface (t). Region F (t) is occupied by an incompressible, viscous fluid and region S (t) is occupied by an elastic solid. The reference fluid and structure domains are denoted by ˆ F and ˆ S , respectively. An example of such fluid-structure interaction problem is shown in Fig. 1. In the following, we will define the solid equations in the Lagrangian framework and the fluid equations in the ALE framework. Denote the structure displacement by η and let the solid domain deformation be a smooth, injective mapping φ S : ˆ S × [0, T ] → S (t) from the reference to the deformed domain, given by
We denote the deformation gradient by F (x, t) = ∇ φ S (x, t) = I + ∇ η(x, t) and its determinant by ˆJ .
To track the deformation of the fluid domain in time, we introduce a smooth, invertible, ALE mapping φ
where ηF denotes the fluid displacement. We assume that ηF equals the structure displacement at the fluid-structure interface, and that it is arbitrarily extended onto the fluid domain ˆ F [17]. To simplify the notation, we will write
vn instead of
whenever we need to integrate a function vn on a domain F (t m ), for m = n.
To model the fluid flow, we use the Navier-Stokes equations in the ALE formulation [38,7,6], given by
where u is the fluid velocity, ρ F is fluid density, σ F is the Cauchy stress tensor and f F is the forcing term. For a Newtonian fluid, the Cauchy stress tensor is given by σ F (u, p) = -p I + 2μ F D(u), where p is the fluid pressure, μ F is the fluid viscosity and D(u) = (∇u + (∇u) T )/2 is the strain rate tensor. Notation ∂ t u| ˆ F denotes the Eulerian description of the ALE field ∂ t u • φF [22], i.e.,
and the domain velocity is denoted by
To model the elastic structure, we use the elastodynamics equations written in the first order form as
(2.3a)
where ξ and η are the solid velocity and displacement, respectively, ρS is the solid density, and Ŝ is the second Piola-Kirchhoff stress tensor. We note that Ŝ is a function of F . The particular choices of Ŝ will be specified in Section 5. The Cauchy stress tensor for the elastic structure is given by
The solid velocity and displacement in the reference configuration are related to their Eulerian counterparts via ξ = ξ • φ S and η = η • φ S , respectively.
To couple the fluid and structure, we prescribe the kinematic and dynamic coupling conditions. The kinematic (no-slip) coupling condition describes the continuity of velocity at the fluid-structure interface, given by u = ξ on (t) × (0, T ).
(2.4)
The dynamic coupling condition describes the continuity of stresses at the fluid-structure interface, prescribed as σ F n F + σ S n S = 0 on (t) × (0, T ), (2.5) where n F and n S are the outward unit normals to the boundaries of the deformed fluid and structure domains, respectively. The fluid and structure problems are complemented with initial and boundary conditions.
In this section, we propose a time-adaptive, partitioned numerical scheme for the fluid-structure interaction based on the strongly-coupled partitioned method developed in [8]. The time integration used in this work is based on the refactorized Cauchy's one-legged 'θ -like' method, which is combined with the partitioned approach used in [50]. The resulting method is adaptive, strongly-coupled, unconditionally stable when θ ∈ [ 1 2 , 1] and second order accurate when θ = 1 2 + O(τ n ), where τ n is the time step.
The main steps in the derivation of the refactorized Cauchy's one-legged 'θ -like' method are summarized as follows.
Let {t n } 0≤n≤N denote mesh points based on a variable time step τ n such that t n+1 = t n + τ n . For an initial value problem y = f (t, y(t)), where we assume that f is smooth, the adaptive Cauchy's one-legged 'θ -like' method is given as
for θ n ∈ [0, 1], where t n+θ n = t n + θ n τ n and y n+θ n = θ n y n+1 + (1θ n )y n . We note that this method is equal to the "classical" θ -method when f is linear. However, for nonlinear problems with variable time-stepping, the Cauchy's one-legged 'θ -like' method is stable, unlike the classical θ -method [12].
It was shown in [12] that problem (3.1) can be written as a sequence of a Backward Euler (BE) step followed by a Forward Euler (FE) step, where the latter can be expressed as a linear extrapolation, resulting in the following system equivalent to (3.1): BE:
Using this approach, the main computational load of the algorithm is related to the computation of the BE step, while computationally inexpensive linear extrapolations increase the accuracy of the scheme. The case when θ n = 1 2 , ∀n, corresponds to the midpoint rule, for which the method is second-order accurate, unconditionally A-stable and conservative.
To decouple the fluid and structure problems, we use similar approach as in [50]. The method studied in [50], based on the Robin-Robin coupling conditions, is non-iterative and unconditionally stable, but sub-optimally, O(τ
), accurate in time. While the refactorized Cauchy's one-legged 'θ -like' method can be applied to easily improve the accuracy of the first-order legacy codes, its direct application to the partitioned approach used in [50] would result in instabilities. Hence, a novel, strongly-coupled, partitioned method was introduced in [8], which was shown to be unconditionally stable for θ ∈ [ 1 2 , 1] and second-order accurate in time when θ = 1 2 . As it was shown in [12] that the adaptive refactorized Cauchy's one-legged 'θ -like' method retains the stability properties of the fixed time-stepping version, we propose to extend the method developed in [8] to include the adaptive time-stepping.
In the following, we denote
The resulting adaptive, partitioned fluid-structure interaction solver is presented as follows.
Algorithm 1. Given θ 0 , u 0 in ˆ F , and η 0 , ξ 0 in ˆ S , we first need to compute
with a second-order method. Then, for all n ≥ 2, compute the following steps:
Step 1 (BE Step). Set the initial guesses as the linearly extrapolated values:
and similarly for ξn+θ n (0) , u n+θ n (0) . The pressure initial guess is defined as
For κ ≥ 0, compute until convergence the backward-Euler partitioned problem: S:
Step 2 (FE Step). Now evaluate the following: S:
Step 3. Compute the local truncation error T n+1 , and given a tolerance δ, and parameters r min , r max and s, adapt the time step:
, and evolve the time interval t n+2 = t n+1 + θ n+1 τ n+1 . Otherwise, set τ n = τ ne w and go back to Step 1.
Step 1 are to be completed with boundary conditions on ˆ S \ ˆ and n+θ n F \ n+θ n . The analysis in [8] can easily be extended to variable time stepping to show that problem defined in Step 1 is linearly convergent and that Algorithm 1 is unconditionally stable for θ ∈ [
1 2 , 1]. The method is second-order accurate for θ = 1 2 + O(τ n ). Hence, the focus of this paper is on the adaptivity aspects of the proposed method. A few remarks are in order. Remark 1. The converged solutions of the structure and fluid BE problems outlined in Step 1, ηn+θ n (κ) , ξn+θ n (κ) , u n+θ n (κ) , p n+θ n (κ) κ→∞ ---→ ηn+θ n , ξn+θ n , u n+θ n , p n+θ n satisfy the following equations: S:
We note that the first problem converges to the structure equations solved with the Neumann interface condition, while the second problem converges to the fluid equations solved with the Dirichlet interface condition. In practice, the BE step of Algorithm 1 will be solved until the relative error between two consecutive approximations of the structure velocity and displacement, and the fluid velocity, is less than a given tolerance .
Remark 2. We note that the nonlinearities in the structure problem related to F and ˆJ are linearized using a fixed-point approach in the BE step of Algorithm 1. However, depending on the definition of the second Piola-Kirchhoff stress tensor, additional nonlinearities may appear, in which case they can be linearized using the same approach. This allows us to treat all nonlinearities in the problem in the same iterative procedure that is used to strongly couple the fluid and solid subproblems. Alternatively, a few iterations of Newton's method can be used to resolve structure nonlinearities.
Remark 3. In the description of Algorithm 1, a simple harmonic extension is used to compute the deformation of the fluid domain. However, this problem can be easily replaced with a different one, such as the linear/nonlinear elasticity equation or other methods mentioned in [53,59,51]. In fact, we use a nonlinear elastic model to realize the ALE mapping in Example 2, presented in Section 5.2.
The formula used to adapt the time step is based on the elementary stepsize selection algorithm commonly used in locally adaptive time-stepping [24]:
.
Numbers r min and r max in (3.2) are added so that the ratio of τ ne w and τ n stays between these values. This type of restriction guarantees the zero-stability of general one-leg variable step size methods (see e.g., [34,35,14,23,28]), and helps to keep the time step from changing too rapidly, which is especially important in stiff problems. The coefficient s ∈ [ 1 2 , 1) is a 'safety' parameter, routinely used to reduce the number of rejected time steps in the adaptive algorithm.
To compute the LTE, we consider two approaches. The first approach, is based on Milne's device [42,39], which is used in predictor-corrector methods [37,Chapter 4], evaluates the error by multiplying with an appropriate factor the difference between the second-order midpoint solution, denoted by y n+1 midpoint , and a modified, explicit, second-order Adams-Bashforth two-step (AB2) method, denoted by y n+1 A B2 . The difference between the modified and the classical AB2 formula is that in the modified AB2 method, the function values are evaluated at half-times, as described in [12], resulting in the following approximation:
.
The LTE for the AB2 method, under localization assumption [27,37], can be written as
where y is the exact solution to the generic problem y = f (t, y), which is approximated in (3.1), and
.
The LTE for the midpoint method, as well as the 'θ -like' method for θ n = 1 2 + 1 2 (τ n ) 2 , is given by
Using (4.1), the LTE of the midpoint method can be written as
Another approach to compute the LTE is by using the Taylor expansions. Namely, the Taylor expansions at t n+1/2 , t n-1/2 and t n-3/2 are used to evaluate y (t n ), as described in [12], resulting in the following formula:
Both approaches detailed above are used and compared in Section 5.
We present three numerical examples used to investigate the performance of the adaptive, partitioned method described in Algorithm 1. Example 1 is based on the method of manufactured solutions, where a simplified, linear model is considered. This example is used to examine which variables should be used to define the LTE, and to compare the two methods for computing the LTE described in Section 4. The performance of the method is then investigated on two more complex problems. In example 2, a moving domain FSI problem is considered, describing the flow in a two-dimensional channel containing an elastic bar attached to a rigid cylinder. Finally, in Example 3, the adaptivity properties of the proposed method are studied on a three-dimensional example representing blood flow in a coronary artery. In all the examples presented here, the problem is discretized in space using the finite element method, and software packages FreeFem++ [30] and NGSolve [49] are used to preform the simulations.
In the first example, we use the method of manufactured solutions to investigate the performance of Algorithm 1 based on the different parameters in the problem. We define the fluid and structure domains as F = [0, 1] × [0, 0.5] and S = [0, 1] × [0.5, 1], respectively. In this example, we assume that the fluid-structure interaction is linear, that the fluid domain is fixed and that the structure is described using a linearly elastic model, given by σ S = 2μ S D(η) + λ S (∇ • η)I, (5.1) where μ S and λ S are Lamé parameters. The linear model is derived with the assumption that F ≈ I , so in this case Ŝ = σ S .
We define the exact solution as
,
We note that the manufactured solution has no gradient discontinuity at the fluid-structure interface. The Neumann boundary conditions are imposed at the external boundaries of the structure domain. At the left and right boundaries of the fluid domain, we also impose the Neumann conditions while the Dirichlet conditions are imposed at the bottom boundary. The boundary conditions and the forcing terms are computed using the exact solution.
The coefficients related to the material properties of the fluid and solid, μ S , λ S , μ F , ρ S , ρ F are all set to 1, and parameter θ is set to 1 2 , in which case we expect the second-order accuracy in time. To discretize the problem in space, we use a 40 × 20 mesh and chose P 1 elements for the fluid pressure and P 2 elements for the fluid velocity, as well as the displacement and velocity of the solid. In this example, we use r min = 0.2, r max = 1.5 and s = 0.95. The final time is set to T = 5 s.
In this section, we investigate the optimal choice of defining the LTE, which is used for time-adaptivity in Algorithm 1.
The initial time step is set to τ 0 = 0.01 s and we use the adaptivity tolerance of δ = 10 -6 . The adaptivity process begins after two iterations have passed. The tolerance which controls the convergence of the subiterations in Step 1, , is set to ε = 10 -4 . The combination parameter is determined using
where H S and H F are the height of the solid and fluid domain, respectively, E is the Young's modulus and ν is the Poisson's ratio. To compute the LTE, we consider both the formula based on the AB2 method (4.1), and the formula based on the Taylor's method (4.3), as described in Section 4. To determine the size of the LTE, we first set
where T η is the LTE for the structure displacement, and T u is the LTE for the fluid velocity. Fig. 2 shows the evolution of the LTEs for the velocity (top-left) and the displacement (bottom-left), the number of trials over time (top-right), and the evolution of the time step (bottom-right) obtained for this choice of T , where LTEs are computed using the AB2 method (4.1). We can observe that T is determined by the LTE for the velocity because of the difference in the magnitudes of T η and T u .
In this case, the simulation ended before t = 1.5 was reached. The LTE for the fluid velocity exceeded the given tolerance and then started decreasing very slowly while the time step was adapting accordingly. Finally, the time step decreased to a machine zero, which caused the simulations to break before the LTE decreased below the given tolerance. As the time step decreased to zero, the LTE for the displacement decreased rapidly as well. The similar behavior is also obtained when the LTE is computed using the Taylor's method (4.3), however the figures are not included for clearer presentation and in order to avoid redundancy. We note that a similar behavior was observed with this definition of T when time-adaptive FSI problems with thin structures were studied in [9]. Hence, we consider redefining T by considering the LTE for the structure displacement only, T = T η . The results obtained with this selection are shown in Fig. 3. In this case, the simulation is completed successfully using both AB2 and Taylor's method. While the adaptivity is controlled by the LTE for the displacement, we note that the LTE for the velocity remains bounded in both cases since the problem is coupled. Using AB2 method, at most 4 trials are needed to advance in time, with the average of 1.18 trials. When Taylor's method is used, no more than 2 trials are needed at any point in time, with the average of 1.4815 trials. The Taylor's method results in larger values of the time step but also larger number of average trials compared to the AB2 method. In both cases, the time step does not approach zero. Based on our results in this section, we use T = T η in the remainder of the paper.
The combination parameter, α, used in the partitioning between the fluid and solid problems, is often determined heuristically, for example using (5.2). Since there are a variety of ways in which it can be computed, in this section we focus on the effects of the different values of α on the performance of Algorithm 1. We set τ 0 = 0.01 as in Example 5.1.1, and use δ = 10 -6 . It was shown in [8] that needs to be reduced at the same rate as the time step in order to obtain the optimal rate of convergence, so here we use the following formula to compute :
starting from 0 = 10 -4 . While the formula is based on (3.2), an outer min{} is added to prevent the tolerance from becoming too large, which could have a significant effect on the accuracy.
To test the effects of α, we use values of 100, 500, 1000, and α n F , which is adapted together with the time step using
(5.4)
The LTE is defined as T = T η and computed using formulas based on both AB2 and Taylor's methods, which are detailed in Section 4.
Table 1 shows the average number of subiterations required in Step 1 and the average number of trials used in the time adaptivity. A wide range of the number of subiterations is observed for both AB2 and Taylor's methods, where in both cases, the smallest value is obtained when α = α n . The number of trials used in the time adaptivity shows very small variability when different values of α are used. In fact, the only slight difference in the results is obtained when AB2 method is used with α = α n .
Our results indicate that formula (5.4) provides a good choice for the combination parameter. Hence, this definition will be used in the rest of this work.
In this section, we compare the results obtained using the variable time-stepping described in Algorithm 1 to the results obtained using the same algorithm, but with a fixed time step. To obtain a fair comparison, the simulations with the adaptive method are performed first. Then, the time step for the fixed time-stepping simulations is computed by dividing the final time reached using the adaptive time-stepping by the total number of trials (including the rejected trials).
In the adaptive algorithm, we use the initial time-step of τ 0 = 0.01, and the initial tolerance of ε 0 = 10 -4 , which is then adapted using (5.3). The combination parameter, α, is computed according to (5.4). The simulations are run until T = 10 is reached using adaptivity tolerances δ = 10 -6 and 10 -10 , and both AB2 and Taylor's methods to compute the LTE. In each case, we record the L 2errors for the fluid and structure velocities, e u and e ξ , respectively, the L 2 -error for the structure displacement, e η,L 2 , and the error for the structure displacement in the energy norm, e η,S , defined as follows:
,
where
The comparison of the relative errors obtained with the fixed and variable time-stepping algorithms is shown in Table 2. For all the adaptivity tolerances and methods for computing the LTE considered in this example, we obtain smaller errors when the adaptive time-stepping is used, compared to the fixed time-stepping. The smallest differences are obtained for the structure displacement in the energy norm. However, significant differences are observed for the other errors, indicating the efficiency of the adaptive time-stepping. When Taylor's method is used, a significant difference in the errors is obtained for δ = 10 -6 . However, the errors become much closer in size when the tolerance is reduced to δ = 10 -10 . While a similar trend is also observed using the AB2 method, the differences obtained using the smaller tolerance are more significant than when Taylor's method is used. Hence, in the remainder of this work, AB2 method based on T = T η will be used to compute the LTE.
This example can also serve to justify the computational cost of the variable time stepping. Indeed, our results indicate that more steps would have to be used when the time step is fixed in order to obtain the same error as with the variable time-stepping algorithm, even when rejected trials are taken into account.
To numerically compute the rates of convergence, we first use Algorithm 1 with τ 0 = 0.01, 0 = 10 -4 and α = 1 on a 10 × 5 mesh, and we use the adaptivity tolerance of δ = 10 -4 . The temporal mesh produced in this simulation using accepted trials consists of 21 points. Starting from this mesh, we create three more temporal meshes, each obtained by splitting each interval of the previous mesh in half. In the same way, we refine the spatial mesh. The resulting spatial and temporal meshes, paired corresponding to the level of refinement, are used to compute the rates of convergence. The relative errors for the fluid velocity and structure displacement, obtained using a range of values for θ , are shown in Fig. 4. As expected, we obtain first-order convergence when θ = 1, and second-order convergence for θ = 0.5. Furthermore, as θ decreases from 1 to 0.5, the errors decrease as well.
Here we consider the classical benchmark problem proposed by Turek and Hron [54] which models a two-dimensional incompressible channel flow around a rigid cylinder with an attached nonlinearly elastic bar. This moving-domain problem features large deformations, and is modeled using the Navier-Stokes equations for the fluid, and nonlinear structure equations for the solid with St. Venant-Kirchhoff materials:
The geometric setup is given in Fig. 5. For the boundary conditions, a parabolic velocity profile is prescribed at the left channel inflow
where u(0, y) = 1.5 Ū y(Hy) (H/2) 2 = 61.5 Ū 4.0 0.1681 y(0.41y), the stress-free boundary condition is prescribed at the right outflow boundary, and the no-slip condition is prescribed on all the other boundary parts, i.e. the circle and the top and bottom walls. We consider the FSI2 problem described in [54], where the physical parameters are given as follows:
A fully developed periodic flow appears around time t = 10. Here we take the final simulation time as T = 15. The quantities of interest are the displacement of the control point A at the end of the beam structure (see Fig. 5), and the lift and drag forces acting on the cylinder and the beam structure:
where S denotes the boundary between the fluid domain and the cylinder together with the elastic structure. We apply Algorithm 1 on a triangular mesh with 2017 fluid elements and 126 structure elements; see Fig. τ 0 = τ 1 = 2 × 10 -3 , with time adaptivity starting from the third iteration. The initial stopping tolerance of the subiterations in Step 1 is set to be 0 = 10 -4 , and is dynamically updated using (5.3). We take combination parameter to be α n , given by (5.2), and set θ n = 0.5 + τ n . Furthermore, the Aitken's 2 method [55] is used to accelerate the fixed point iteration in
Step 1, where the acceleration is performed separately on the fluid velocity update and the pressure update. A total of 2275 trials are used in the simulation, with 50 rejected trials. The time evolution of the quantities of interest is shown in Fig. 7, with zoom-in results in the fully developed flow regime on the time interval t ∈ [14,15] shown in Fig. 8, where we also present the reference data from [54]. Our results agree well with the reference data. The evolution of the time step size is shown in Fig. 9. From this figure, we observe that the time step size goes through three different stages: in the initial stage, it gradually increases from τ = 2 × 10 -3 to about τ ≈ 0.1 (at around time t = 1); in the transition stage, the step size gradually decreases to about τ = 0.04 (at around time t = 3) with a couple of rejected trials, and then further decreases in an oscillatory fashion to τ ≈ 4 × 10 -3 (at around time t = 10); in the final stage where the flow is fully developed, the time step size oscillates periodically around τ = 3.5 × 10 -3 , where the periodic time step size pattern can be seen more clearly on the right panel of Fig. 9. In each of these periods, two trials are rejected. The average number of subiterations required for convergence in Step 1 (corresponding to κ) is about 13.8, and the average stopping tolerance is ≈ 7.4 × 10 -6 . We present in Fig. 10 the evolution of the number of subiterations and the stopping tolerance . The evolution of the displacement error indicator T η and the fluid velocity error indicator T u are shown in Fig. 11. We observe that T η < δ = 5 × 10 -7 as expected, and that the fluid velocity error indicator, whose magnitude is about 4 × 10 -4
in the fully developed regime, is also under control.
We note that in [54], a constant time step of 2 • 10 -3 was used for the same problem, resulting in 7500 steps. Even with the rejected trials taken into account, only 2275 trials are needed in our simulation, giving the average time step of 6.59 • 10 -3 . Hence, the adaptive time-stepping in this example requires only about 30% of steps needed when the fixed stepping is used. Fluid density ρ f (g/cm 3 )
Wall density ρ s (g/cm 3 )
1 .1
In this example, we focus on a three-dimensional, simplified model of blood flow in a common carotid artery. The computational domain is a straight cylinder of length 4 cm. The fluid domain has a radius of 0.3 cm, and the thickness of the structure domain is 0.06 cm, see Fig. 12. The arterial wall elastodynamics are described using a linear model (5.1), and the blood flow is modeled using the Navier-Stokes equations (2.1)-(2.2). Even though a linear model is used for the structure, a moving domain problem is considered. Values of the parameters used in this example are given in Table 3.
At the fluid inlet section, we prescribe a fully developed time-dependent axial velocity, and zero normal stress is imposed at the outlet. The boundary conditions [58,46] are given as
at inlet and σ n = 0 at outlet, (5.5) where r(t) is defined as
and u D (t) is shown in Fig. 13 for t ∈ [0, 1]. Afterwards, it is repeated periodically until the final time (T = 6 s) is reached.
We assume that the structure is pinned at the edges, and that zero stress is imposed at the solid external boundary, η = 0 at inlet and outlet, σ S n = 0 at structure external boundary.
All initial conditions are set to zero. The fluid mesh used in this example consists of 10320 tetrahedral elements, while the structure mesh consists of 14400 tetrahedral elements. We initially set τ 0 = 10 -3 , and keep this value for 20 iterations, when the time adaptivity process starts. We use 0 = 10 -4 , θ = 0.5, δ = 4 • 10 -8 , r min = 0.2, r max = 1.5, and s = 0.95. The combination parameter, α, is computed in each step using (5.4). A total of 2578 trials are used in the simulation, with 689 rejected trials, giving a mean value of 1.36 trials per step. The time evolution of the local truncation errors for the structure displacement and the fluid velocity is shown in Fig. 14. We note that the errors are divided by the characteristic size of the corresponding variables, E = 1.7 • 10 -4 and U = 50. In both cases, a periodic-like behavior is observed. The local truncation error for the structure displacement oscillates throughout the simulation in a similar manner in each cycle. While the error for the fluid velocity shows similar shape in each period, its magnitude is slowly growing in time.
The evolution of the time step size and the number of subiterations in Step 1 are shown in Fig. 15. The time step size also develops a periodic-like behavior, where the time step drops to roughly 6 • 10 -5 when the inflow boundary condition begins to increase at the beginning of systole, and then increases to 7 • 10 -4 or more during the diastolic part of the cycle.
On average 20.9 subiterations are needed in Step 1 of the proposed algorithm. We note that the Aitken's acceleration was not used in this example. Finally, Fig. 16 shows a comparison of the displacement magnitude, velocity magnitude and pressure obtained using the adaptive and fixed time stepping. In each method, the variables are obtained at a point on the fluid-structure interface located in the center of the domain, shown in Fig. 12. To obtain the resulting using fixed time-stepping, Algorithm 1 was used with the fixed time step of τ = 3 • 10 -3 , which was obtained as an average value of the time step sizes used in the adaptive algorithm. In all three cases, an excellent agreement between the fixed and adaptive time-stepping is observed.
An adaptive, partitioned method for FSI problems with thick structures is presented and investigated in this work. The method is based on the Cauchy's one-legged 'θ -like' method, which is A-stable for θ ∈ [ 1 2 , 1] when used with adaptive time steps. When θ = 1 2 , the method corresponds to the midpoint rule, in which case the method is conservative and second-order accurate. The second-order accuracy is maintained also for θ = 1 2 + O(τ n ), in which case a small amount of dissipation is introduced. The Cauchy's method is combined with a partitioning approach, resulting in a strongly-coupled scheme. The main focus of this work is on the adaptivity aspects of the proposed algorithm. In particular, we investigated different choices of variables used in the calculation of the LTE, as well as two different formulas for the computation of LTE. Our results indicate that the best results are obtained when the LTE is calculated using the structure displacement only, as the LTE for the displacement is most sensitive to the changes in the time step size. The formulas for computing the LTE based on the AB2 and Taylor's methods were compared. Both methods gave comparably good results, however, the AB2 method slightly outperformed the Taylor's method. Using the displacement-based LTE, calculated using the AB2 formula, the performance of the method was tested on a two-dimensional benchmark problem describing the channel flow around a cylinder with an elastic flap, and on a three-dimensional problem describing blood flow in a coronary artery.
In both problems, the results obtained using the adaptive time stepping algorithm showed excelled agreement with the reference data. Even though the LTE used to adapt the time step is based on the displacement, the LTE for the fluid velocity remains bounded in the first two examples. While in the third example the velocity LTE shows linear growth in time, the errors between the values obtained using the adaptive and fixed time-stepping do not exhibit the same behavior. In Example 2, the time step used in the first part of the simulation is more than one order of magnitude larger than the time step used in the last part, when periodic regime was developed. Similarly, in Example 3, the time step used during systole is at least an order of magnitude smaller than the time step used during diastole. This demonstrates the potential for computational savings that can be achieved using adaptive time-stepping in problems with varying dynamics, while controlling the accuracy of the solution.
Journal of Computational Physics 473 (2023) 111708