A fast front-tracking approach and its analysis for a temporal multiscale flow problem with a fractional-order boundary growth

This paper is concerned with a blood flow problem coupled with a slow plaque growth at the artery wall. In the model, the micro (fast) system is the Navier-Stokes equation with a periodically applied force and the macro (slow) system is a fractional reaction equation, which is used to describe the plaque growth with memory effect. We construct an auxiliary temporal periodic problem and an effective time-average equation to approximate the original problem and analyze the approximation error of the corresponding linearized PDE (Stokes) system, where the simple front-tracking technique is used to update the slow moving boundary. An effective multiscale method is then designed based on the approximate problem and the front tracking framework. We also present a temporal finite difference scheme with a spatial continuous finite element method and analyze its temporal discrete error. Furthermore, a fast iterative procedure is designed to find the initial value of the temporal periodic problem and its convergence is analyzed as well. Our designed front-tracking framework and the iterative procedure for solving the temporal periodic problem make it easy to implement the multiscale method on existing PDE solving software. The numerical method is implemented by a combination of the finite element platform COMSOL Multiphysics and the mainstream software MATLAB, which significantly reduce the programming effort and easily handle the fluid-structure interaction, especially moving boundaries with more complex geometries. We present some numerical examples of ODEs and 2-D Navier-Stokes system to demonstrate the effectiveness of the multiscale method. Finally, we have a numerical experiment on the plaque growth problem and discuss the physical implication of the fractional order parameter.


Introduction
Multiscale problems have been extensively studied in the past two decades.For spatial multiscale problems, people have developed computable models such as quasi-continuum or atomistic-tocontinuum coupling (AtC) models and QM (quantum mechanics) and MM (molecular mechanics) coupling to simulate material behaviors [40,20,35,22,25,41].For temporal multiscale problems, an example is chemical reactions with concentrations of the species varying from seconds to hours while the time scale of the oscillations of the chemical bonds is in the order of femtoseconds [2].For a general introduction to multiscale methods, we refer to [5,6,39].
A common challenge in simulating these multiscale problems is the enormous computational cost when the microscale feature needs to be resolved and corresponding microscopic discretization is performed, or alternatively, the loss of microscopic information if the macroscopic discretization is performed.Macroscopic and microscopic processes should be properly coupled in order to solve such problems effectively and accurately.
The heterogeneous multiscale method (HMM) is one of the most prominent techniques to deal with multiscale problems, which relies on an efficient coupling between the macroscopic and microscopic models [4,1,7].For temporal multiscale problems with time scale separation, the macroscale quantities can be computed from the microscale subproblem, and large time steps can be employed to solve the macro-scale model in order to save the computational cost.HMM for temporal multiscale problems only considers local solutions of the microscopic subproblem, thus the initial condition on the local interval needs to be carefully designed and depend on some prior knowledge of the microscale behaviour of the system.
We shall consider in this paper a temporal multiscale problem of the atherosclerosis with a commonly slow plaque growth along the artery boundaries.Frei and Richter [8] pioneered the study in this direction and presented a basic model of two-way coupled blood and plaque growth in blood arteries, where its numerical simulation is carried out in the Arbitrary Lagrangian-Eulerian (ALE) framework.The ALE may complicate not only governing equations and the analysis of the fluid structure interaction but possibly also the treatment of more complex boundary growth of the plaque.The analysis of the multiscale method is thus done in [8] for a largely simplified coupled ODE system from the ALE transformed blood-flow-plaque-growth model.In this paper we propose to track the changing domain directly using the front-tracking approach instead of ALE at each time step, thus governing equations are not changed and the method may be more handily applied to general dynamic growth of the plaque.The front-tracking framework not only simplifies the design and analysis of the multiscale method in its original PDE form but also makes it easy to use existing PDE solving software to implement the developed algorithm.
Furthermore, we adopt a more general plaque growth process containing fractional derivatives where memory effects of the plaque accumulation or evolution may be included.Fractional calculus has been extensively studied in the last two decades, especially in the fields of fluid mechanics [37,3,18] and anomalous diffusion [16,21,33,42].Compared with integer order operator, the fractional order operators have a non-local structure, and are suitable for describing the memory and hereditary properties of many physical processes.For our applications, macrophages in the artery wall take up low density lipoproteins (LDL), which carry cholesterol and triglycerides to the tissues, and are finally transformed into foam cells, which are engorged with lipids [13].In the long term, macrophages and foam cells in the artery wall are influenced by a variety of other cells, and thus the diffusion is most likely to be anomalous [43,19].Therefore, the fractional operator with memory effect may be more suitable than the local integer operator to describe the anomalous diffusion process in the artery wall [23].More realistic plaque growth equations can be found in [43] and [44].
In this paper, we consider the blood flow problem with the atherosclerosis, where the incompressible Navier-Stokes equation with a time-periodic force is coupled with a fractional plaque growth model.Due to the slow plaque growth a significant long-term computation is necessary to observe the change in the domain and the flow properties.The nonlocal property of the fractional order operator makes such a long-term computation impossible.For an efficient long-term computation, it is necessary to develop a temporal multiscale method.The 2-D fluid structure interaction problem in this work faces the challenge that the computational domain changes with time (plaque grows slowly with time) in a fractional order and that the multiscale error analysis will be significantly more difficult in comparison to a simplified integer order ODE system in [8].We shall first simplify the procedure by directly tracking the changing domain at each time step using the front-tracking approach and then formulate an auxilary time-periodic system and an effective time-average equation.Based on the auxilary system, a multiscale method is then developed to deal with the macro (slow) and micro (fast) equations separately and two scale variables interact through the growing boundary so as to reduce the computational cost.A simple finite difference scheme in time and a finite element method with an adaptive mesh near the time-dependent boundary will be used to solve both the original and the multiscale method.We shall also introduce a fast iterative procedure to find the initial value of the time-periodic flow subproblem and analyze its convergence rate.The front-tracking framework, designed discrete schemes and the iterative procedure for solving the temporal periodic problem make it easy to implement the multiscale method with existing finite element software.The multiscale method is then implemented through a combination of COMSOL Multiphysics [29,26] and MATALB.The numerical framework may be applied to a wide range of problems with a periodic applied force and a slow boundary growth.

Outline
The rest of this manuscript is organized as follows.In Section 2, we describe the mathematical model and make necessary assumptions.In Section 3, we derive an time-periodic subproblem of the flow equations to approximate the original problem, and analyze the error of the temporal multiscale system at the continuous level.In Section 4, we present a time discretization scheme and implementation details of the temporal multiscale system, and analyze the error of its timediscrete scheme.An iterative method to find the initial value of the time-periodic flow equations is also shown in this Section.In Section 5, we demonstrate and validate the accuracy and efficiency of our multiscale method through several numerical examples.The effect of the fractional order parameter on plaque growth is also investigated.Finally, we conclude the paper in Section 6.

Notation
For domain Ω and m ≥ 0, we use the standard notation for the Sobolev space H m (Ω) and the Banach space L m (Ω).We use (•, •) to denote the inner product in L 2 (Ω).Throughout this paper, the letter C will denote a positive constant, with or without subscript, its value may change in different occasions.

Model problem
We consider the model which describes the biochemical processes leading to the growth of plaque in blood arteries, as shown in Figure .1.We assume that the plaque growth occurs only at the upper boundary of the blood artery, which is controlled by the concentration variable u(t).The blood flow is modeled as an incompressible Newtonian fluid, which is suitable for the description of large arteries [30].The two-way coupled model of blood flow and plaque growth with fractional derivatives is given as follows For simplicity of analysis we consider Dirichlet boundary condition, though the method developed here may be applied to other common boundary conditions.In (1), the velocity v and concentration u represent the micro (fast) variable and the macro (slow) variable, respectively.ρ is the density of blood.A periodic force f(t) = f(t + 1) is applied to the flow due to the periodic nature of heart pulse.The Cauchy stress tensor σ is defined as where ν is the kinematic viscosity.The function γ(u, x) characterizes the shape change of Ω with respect to time through the concentration u(t).
The reaction term R ≥ 0 describes the influence of wall shear stress on the boundary growth, which can be simplified as follows [9,8] where n denotes the outward facing unit normal vector at the deformation boundary ∂Ω.ε 1 is a small parameter that controls the change of u.D α 0 + is the Caputo fractional derivative of order 0 < α < 1 denoted by [27] The Riemann-Liouville (R-L) fractional integral for α ∈ (0, 1) on finite interval [0, T ] is defined as From the definition of the R-L integral and the Caputo derivative [27], for α ∈ (0, 1), we can obtain

Assumptions
Next, we present the essential assumptions which ensure the existence of solutions.
We assume that the incompressible Navier-Stokes equations on the moving domain Ω(u(t)) have a solution v(t) ∈ H 2 (Ω(t)) and p(t) ∈ H 1 (Ω(t)).The reaction term is bounded and has the following Lipschitz condition with respect to slow and fast variables Remark 1.The reaction term R given in (3) satisfies the above assumptions, and its proof can be seen in [8].
Remark 2. We would like to point out that u(t) is bounded based on the properties of fractional operators.Applying the operator I α 0 + on both sides of the third equation of (1), we obtain So to see significant (or O(1)) boundary growth we need to compute up to T = O(ε − 1 α ).
Assumption 2. Let u ∈ C 1 [0, T ] be fixed.We assume that the following incompressible Navier-Stokes equation have a unique periodic solution (v u , p u ) and the solutions are uniformly bounded Remark 3.For a fixed flow domain, the uniqueness of the periodic solution is guaranteed for moderate Reynolds numbers [8,10].The periodic Navier-Stokes system can serve as an auxiliary problem, which allows us to quickly solve the temporal multiscale problems.
We have the following assumption for the time changing shape function γ(u, x).
Assumption 3. Plaque growth is due to the increased concentration of foam cells [36], and the plaque growth process is irreversible.We present the following assumptions ∂γ ∂u ≥ 0, and ∂γ ∂t ≥ 0.
We remark that in the integer order case ∂γ ∂t ≥ 0 can be derived from (8) and ∂γ ∂u covers the case considered in [8].

Derivation and analysis of the fractional multiscale problem
In this section, we derive effective time-average equations for the temporal multiscale system with fractional plaque growth (1) based on the assumptions in the previous section.We note that the error analysis between the effective equation and the original equation ( 1) is performed for the Stokes problem.This is substantially different from the highly simplified system of ODEs [8], and can be extended to the full Navier-Stokes system in (1).

Derivation of the effective equation
According to the properties of fractional derivatives, we introduce a new variable where I α 0 + is the R-L fractional integral operator.By inserting R(v(s), U (t)) in ( 14), we have Proof.By direct calculation, we have Applying the operator I α 0 + on both sides of the equation ( 17), we have which completes the proof.
Applying equation ( 6), we have Let where C L33 depends on α, and the constants in Assumptions 1.
Proof.By using the Lipschitz condition of the reaction term R, Lemma 1 and Lemma 2, we have We thus have the following estimate for the equation ( 14) of U (t), The discretization of U (t) in a macro-time step T n → T n+1 = T n + T is not accurate because it involves the dynamic evolution of v(T n ) to v(T n+1 ) on the fast scale.The local periodicity can be helpful to effectively capture the velocity feature on the fast scale in [T n , T n+1 ].We thus introduce an auxiliary periodic problem at the end of this subsection whose solution v U (t) will be used to approximate the fast scale velocity feature.Here U (t) is an approximation of U (t) to be defined at the end of the Section.
Next, we approximate the fractional differential equation ( 25) by inserting R(v U (t) (s), U (t)) with a fixed U (t), that is, writing (25) as By using the Lipschitz condition of R, we have In the Theorem 1, we will prove that and therefore we have Here we remark that for integer order case (α = 1), the estimates can be much better and can be done much easier, that is, the estimate of ( 20), ( 27) and ( 28) is O(ε) and that of ( 29) is O(ε 2 ).See a relevant Remark 5 later.
It has been shown that blood flow would be approximately periodic with period 1 after a period of time [34,15].In other words, even if, for a given initial velocity, the blood flow is not immediately time-periodic, it would become periodic after a certain time.We assume that the blood flow has reached a periodic state initially at time zero since the blood flow has already run for a long time in the body.Therefore, this paper considers the initial value flow problem (1) that has already been time periodic of period 1 from the initial time when the boundary is fixed with initial u(0) = u 0 , in other word, we may assume v(0) = v u(0) (0) (the definition of v u (0) is given in (11)).This assumption will be used later in the analysis of Lemma 5.
For the temporal multiscale problem of (1), we give the following auxiliary time periodic flow problem based on (29).
From this, we give the framework diagram (Figure .2) of the multiscale algorithm.
) Stepping the macro variable U (t) with macro step size ∆T 1, and then computing the periodic solution v U (t) using the micro step size ∆t over a unit cycle for a fixed flow domain Ω(U ).

Error analysis of the fractional multiscale problem
In this section, we will carry out a temporal multiscale error analysis for the Stokes problem.We will show that problem (30) can be used as an effective approximation of problem (1).Lemma 4. Let the concentration u be fixed and 0 ≤ u ≤ u max , v u is the solution of the corresponding time periodic Stokes problem It holds that ≤ C L34a , and Here C L34a and C L34b depend on f, domain Ω, the constants in Assumptions 2 and 3.
Proof.Based on the Stokes equation and Assumption 2, we have For domain Ω(u), we can derive from a fixed reference domain Ω(0) using the ALE method [31] T : with Jacobian matrix and determinant given by The Stokes equations can be mapped to a fixed reference domain Two different representations of the Stokes (Navier-Stokes) equation in the ALE and in the Eulerian coordinates are equivalent (see [31], Chapter 2), it follows that Differentiating the equation (36) with respect to u.Let w u = ∂vu ∂u , we obtain Let 38) can be rewritten as with time periodic condition z u (0) = z u (1) and homogeneous Dirichlet boundary condition.Here, It is easy to check that both the determinant J and the elements associated with the matrix F are non-negative and bounded from below and above, Multiplying the second equation of ( 39) by z u and integrating with respect to the space variables on the domain Ω(0), we get Combining Young's inequality, Poincaré inequality, ( 37) and ( 40), we can derive the following inequalities, where γ is the Poincaré constant.Multiplying e C1t on both sides and integrating in t and using the periodic condition of z u , we obtain Integrating (42) in t on [0,1], we obtain Multiplying the second equation of (39) by −P ∇ • J∇z u F −1 F −T , where P is the L 2orthogonal projection from L 2 (Ω(0)) 2 to H, and Integrating on the domain Ω(0), using the properties of the projection (see Section. 2 in [32]) and Young inequality, we have with time periodic condition . Integrating (46) from 0 to 1 and following the proof of in [32,Theorem 2.22], we obtain Based on the relationship between w u and z u and with the help of equation ( 37), we have and The proof is completed.The following lemma gives the error estimate between the solution v(t) of the original problem (1) and the time periodic solution v u(t) (t).For a fixed u(t), we have the family of periodic solutions We note here that although v u(t) (s) is defined on [0, 1], it can be periodically extended to [0, T ].
Lemma 5. Let u ∈ C 1 [0, T ], v(t) be the solution of the original Stokes problem and the initial values of v and v u(t) are identical, i.e., v u(0 where C L36a and C L36b depend on C L34a and the constants in Assumptions 1. ξ = min{1, νe αγ } and γ is the the Poincaré constant. Proof.For v u(t) (t), by the chain rule, it holds that Combining (50), v u(t) satisfies the following PDE Multiplying (54) by w and then integrating on the domain Ω(u(t)) yields For the first term, with the help of Assumption 3, and integrating over the domain Ω(u(t)), we have Combining Poincaré inequality, Young's inequality, Lemma 4. Integrating from 0 to t and using Lemma 2, we obtain where λ = ν γ , γ is the the Poincaré constant and ξ = min{1, νe αγ }.Multiplying (54) by w t and integrating on the domain Ω(u(t)), we have It thus follows that Multiplying (54) by −P ∆w, P is the orthogonal projection from L 2 (Ω(u(t))) 2 to G G = w; div w = 0 in Ω(u(t)) and w Integrating on the domain Ω(u(t)), we obtain The proof concludes with integrating (61) in t on [t, t + 1] and using (59), Next, we prove a theorem to show the error estimation between the auxiliary time periodic problem (30) and the original problem (1).
u) be defined as the Stokes problem corresponding to (1) and (v U (t) , U (t)) be defined as the Stokes problem corresponding to (30).With the initial value C T 37a and C T 37b depend on C L34b , C L36a , C L36b and the constants in Assumption 1. ξ follows the definition in Lemma 5.
Using Assumption 1, we obtain For the first term, using (24), we have For the second term, combining Lemma 4, Lemma 5 and Cauchy-Schwarz inequality, we have Taking I α 0 + on both sides of (65), for T = O(ε − 1 α ), we obtain Applying Gronwall's inequality, it is easy to obtain Using (24), Lemma 4 and Lemma 5, we obtain (64) is estimated by inserting v u(t) (t) and using Lemma 4 and Lemma 5, we obtain The proof is completed.
Remark 5.In the integer order case (α = 1), it is obvious that u (t) = O(ε), the order of ( 57) and (62) are O(ε 2 ) by using Young's inequality.Following the proof of Theorem 1, we can easily obtain the order of ( 63) and ( 64) is O(ε).

Numerical methods
In this section we introduce simple discrete schemes to approximate this two-way coupled multiscale problem.For the time discretization of fast scale, we use the first-order linear backward Euler scheme.The finite element method is used for spatial discretization, the velocity and pressure are approximated by P2−P1 elements.For the slow scale, the discretization is based on the L1 scheme (See [11]).We can easily choose a higher-order scheme for discretization, but the computational cost of long-term simulation are prohibitive.

Direct method
For fractional differential equations, writes the L1 scheme [11,38,24] to discretize the Caputo derivative of order 0 < α < 1 where a j = (j + 1) 1−α − j 1−α (j ≥ 0) and |K(u The time discretization scheme of ( 1) is then as follows The L1 scheme is an implicit method.For simplicity, we change it to an explicit method and the order is O(∆t).
The direct method is used for the forward simulation, and the algorithm is as follows.
Algorithm 1 Algorithm for the direct method.Let v(0) = v 0 and u(0) = u 0 .for i = 1 : M do Step 1. Solve the third equation of (75) to obtain u i .
Step 2. Based on u i obtained in Step 1, a fixed domain is generated and meshed.
Step 3. Solve the coupled Navier-Stokes equation (the first and the second equation of (75) ) to obtain v i .end for

Multiscale method
For the auxiliary time periodic problem (30), we propose the semi-discretization scheme of U The semi-discretization scheme of the time periodic Navier-Stokes equations is as follows To solve the multiscale problem (30), it is necessary to identify the time-periodic solution.We first prove a lemma which provides an iterative method to find the time-periodic solution.Although the lemma is for the continuous Navier-Stokes equations, it can be proved in a very similar way for the temporal discrete scheme of the Navier-Stokes equations.Our numerical experiments demonstrate that the iterative method based on this lemma works very well in finding the initial value of the time-periodic problem of the Navier-Stokes equations.Lemma 6.For a fixed U (flow domain is fixed), let v U (t) be the solution of the initial value problem and v U (t) be the time-periodic solution.Let v U 0 = v 0 be the initial trial value.We have the following result In other words the solution of the initial value problem fastly approaches the initial value of the periodic problem under the uniqueness assumption (Assumption 2).
Proof.For a fixed U , let w(t) = v U (t) − v U (t).We have the following governing equations Multiplying the second equation of ( 79) by w and integrating with respect to the space variables on the domain Ω(U ), using Poincaré inequality, we obtain Multiplying ( 80) by e Ct , then integrating from t − 1 to t, we have For a constant 0 < e −C < 1, we obtain Based on Lemma 6, we can give an initial trial value (usually the inflow velocity) and the tolerance τ > 0 to reach the initial value (and the end value) of the periodic solution, and then perform the calculation.In this process, we calculate the error every 1 second until the given tolerance is reached.The algorithm is as follows Algorithm 2 Algorithm for the identification of the initial value of time-periodic solutions.
Given an inflow velocity v 0 , let τ > 0 be a given tolerance and let n = 0.
Besides the fast convergence we would like to point out that Algorithm 4.2 also makes it easy or effective to implement the solver of Navier-Stokes with the software COMSOL Multiphysics.Now, the fast and slow scales of (77) have been separated.Having Algorithm 2 to find the initial value for the time-periodic problem (77), we introduce the following multiscale algorithm.
Step 1.For a fixed U j−1 , solve the time periodic auxiliary problem (77) to obtain (v Uj−1 ) i based on Algorithm 2.
Step 2. Calculate the integral reaction term in equation ( 76) Step 3.
Step forward U j−1 → U j and go to Step 1 with the L1 approximation

Implementations
In the ALE formulation, the moving boundary may be transfered into a fixed boundary and then u(t) which defines the moving boundary enters the flow equations, which are usually much more complex than the original flow equations (See, for example, (36) after the ALE transformation).
In [8] their numerical algorithm is based on the ALE and these ALE transferred flow equations are further simplified into quasilinear ODEs in their theoretical analysis.As seen earlier, our algorithms are based on simple front-tracking and the flow equations remain as their original simple form.This plus the finite different scheme and the iterative algorithm to identify the initial value of the timeperiodic problem (See Section 4.1) makes it particularly easy to use existing software to implement our algorithms, for example, a combination of the finite element software COMSOL Multiphysics 5.6 and MATLAB.In our numerical experiments later, we will test examples of 2-D Navier-Stokes problem, we connect the finite element software COMSOL Multiphysics 5.6 and MATLAB 2016b to realize the coupled computation.The time discretization is performed in MATLAB, and the spatial discretization is handled in COMSOL by using adaptive mesh.Specifically, in the multiscale computation, we use MATLAB to solve the concentration U , and fix the flow domain with the U .
Then we pass the domain information to COMSOL and use it to find the time-periodic solution v U , and the obtained v U is returned to MATLAB to complete an implementation cycle.This is very effective for our numerical experiments of direct and multiscale algorithms, and saves time in programming.This procedure may also be directly applied to engineering applications.
To verify the performance of the multiscale method, we compare its solution with that of the original initial value problem (1) (called direct solution).In doing this direct computation, we set the time step t = 1/20 and use 210 spatial elements in the domain.For the multiscale method, the number of spatial elements remains the same, and we change the macro-scale time step ∆T and the micro-scale time step ∆t to test the efficiency and accuracy of the proposed scheme.We would like to point out that the mesh is locally refined at the deformed boundary.Compared with the model with integer derivatives, using the direct method to achieve long-term calculation for the model with fractional derivatives is almost impossible due to the nonlocality (or the integral in the time interval [0, t]), therefore the total computation time T is usually not too large in our error comparison.

Error Analysis of the semi-discrete scheme
There are already a lot of error analysis available for numerical schemes of Navier-Stokes equation (See e.g.[12,14] and references therein), although very few are for the Navier-Stokes equations coupled with a fractional differential equation.In this section we analyse the error of the temporal semi-discrete scheme for the fractional part of the system, assuming that the temporal discrete error estimate has been done for the integer order Navier-Stokes part.In view of this, we make the following assumption.Assumption 4. We assume that τ > 0 is the tolerance in Algorithm 2 to solve the time periodic auxiliary problem.The error of the fast scale first-order linear backward Euler scheme can be expressed as follows, for all time step t i , where C A4 is a constant.Remark 6.We can adapt the temporal error estimates in, for example, [17, equations (3.32)-(3.40)]for the initial value problem of Navier-Stokes equation, by assuming sufficient regularity and compatibility conditions and constant density.For the time-periodic problem the estimates may be obtained as well by additional applying the technique of initial value estimates ( 43)-( 44) to obtain the initial value error first.
Theorem 2. Let u ∈ C 2 [0, T ] and U ∈ C 2 [0, T ] be the solutions to the original problem (1) and the time averaged problem (30), respectively, U j be the solution of the discrete effective equations (77).For T j = O(ε − 1 α ), we have the following error estimate where C T 44 is a positive constant.
Proof.We first consider the regularity estimate of U .Based on the definition of the Caputo fractional derivative, by direct computation, we have Moreover, it is easy to see U (0) = 0, and (d t 1 0 R(v U (t) (s), U (t))ds) t=0 = 0. Therefore, we have the following estimate of U (t) by integrating by parts (87) Applying Gronwall inequality to the above inequality, we obtain Applying Taylor expansions of U (T j ) around T j−1 , Lemma 4, we obtain the error equation With the help of the Lipschitz condition of R, we have Combining Lemma 4, Assumption 4 and Cauchy-Schwarz inequality, we obtain ∆t (89), (90), and the inequality a j−1 > (1 − α)j −α lead to the following inequality, Now, we will prove the following estimate by induction, For j = 1, we can easily obtain |E 1 | ≤ C |E 0 | + ∆t + ∆T ε 1 α + τ .Suppose now (92) holds for 1, 2, ..., j − 1, we need then to prove that it holds also for j.From (91), we obtain which completes the proof of the estimate (92).With E 0 = 0, we obtain The result follows by Theorem 1, we have Remark 7.Here we assume 1 2 ≤ α < 1 in the theorem is only to avoid the weak singularity at t = 0 in derivtives of U (t) (see ( 86)).If 0 < α < 1/2, we expect that there will generally be a weak singularity in derivatives of U (t) near t = 0 and the convergence order of the L1 finite difference scheme will be affected there too.A further discussion on this is beyond the scope of this paper.We refer to [38] (Section 4) for a discussion on using a so-called graded temporal steps near t = 0 to improve the convergence order near the weak singularity.

Numerical experiments
In this section, we carry out two numerical experiments to test the accuracy and effectiveness of the proposed multiscale method.We first perform a simple numerical test with a system of ODEs as a simplified model, such that the exact solution can be obtained to verify the accuracy of the algorithm.Then we carry out numerical tests on the atherosclerotic problem with a plaque growth modeled by Navier-Stokes equations to show the advantages of the multiscale method.Finally, we evaluate the effect of the fractional parameter α, which describes the growth rate of plaque.

Test of ODEs system
To obtain the exact solutions of v(t) and u(t), the construction of reaction term R is very simple without using the form in (3).
This fractional problem has the exact solution: is the generalized Mittag-Leffler (M-L) function [28].
We set the macro step size ∆T = 1000, micro step size ∆T = 1/100 and trial value (v Uj ) 0 = 0.5 in the simulation.It can be seen from Figure 3 that the numerical results of the multiscale method agree well with the exact solution.This also verifies the accuracy of Algorithm 1-3.The micro equation is discretized using the standard first order backward Euler scheme, and we further choose different macro step size to test the convergence rate of the discrete scheme (76) of the macro equation.Quantitative error measurements and convergence rates are shown in Table 1.It is observed that our numerical method is first-order convergence for macro variable, which is consistent with the results of numerical analysis.

Test of 2-D flow problem
In this subsection, we carry out numerical tests for Navier-Stokes flow.The results obtained by the direct method are used as a reference to verify the effect of the multiscale method and give the convergence rates.
We test the full Navier-Stokes system in (1).The flow domain is shown in Figure 1 with the shape function γ(u, x) = ue −x 2 .
The time periodic Dirichlet condition is set on the left inflow boundary, which appears to make the problem periodic.We set v in = 30 1 − y 2 4 sin 2 (πt) cm/s.On the right outflow boundary we set a frequently used pressure condition: −p n + ρν ∂v ∂ n = 0.The parameters in the fluid and reaction term are the same as in [8] which are claimed to mimic the real fluid-structure interaction problem: ρ = 1 g/cm 3 , ν = 0.04 cm 2 /s and σ 0 = 30 g/(cm s 2 ).Moreover, we set fractional parameter α = 0.6, ε = 8 • 10 −4 s −1 , computation time T = 8000 s and u(0) = u 0 = 0.2.
For the multiscale method, we set the tolerance for the time periodic solution as v U (t + 1) − v U (t) 2 L 2 (Ω(U )) ≤ τ = 10 −6 .We choose different macro step sizes for comparison with the direct method and define Error= |u j − U j |.The errors and CPU time are shown in Figure 6.Even for not so large computing time T = 8000 s, the CPU time to solve such a fluid-structure interaction problem using the direct method is 77 hours, which demonstrates the excellent performance of our multiscale method.
We test the convergence rate of Scheme ( 76)-( 77).The result of numerical analysis (85) shows that the effect of the macro step size ∆T is dominant when ∆T > ε −1 ∆t, and vice versa.We fix the micro step size ∆t = 1/40 and macro step size ∆T = 80 to test the convergence order of ∆T and ∆t, respectively.We then refine the macro time step or micro time step by 2 and calculate the error at T = 8000 s.It is observed from Figure 5 that our numerical method achieves the expected first-order convergence for ∆T and ∆t.
It should be noted that the direct long-term fractional calculation poses challenge to the computer capacity (in terms of CPU time and RAM) since the solution has nonlocal dependence on previous steps.The multiscale method constructed and theoretically justified in this paper significantly reduces the computational cost associated with the micro-scale steps needed in the direct computation and can thus solve the problem very effectively.

The effect of fractional order parameter on plaque growth
Finally, we consider the effect of the fractional parameter α on the atherosclerotic plaque growth and hemodynamics.Here, the vessel width is proportionally enlarged to 3 cm, and the variable domain is given by Ω(u) = (x, y) : |x| < 5cm, −1.5 cm < y < (1.5 − ue −x 2 ) cm .
The velocity magnitude and plaque growth when the fractional order parameter α = 0.85 and α = 0.95 are shown in Figure 6, with a strong narrowing of the flow domain.We expect that the computational cost would be huge for such a large T by using the direct method since the fractional derivative is involved with an integral from 0 to T .For a long-term simulation, the computational cost of the multiscale method is significantly reduced.Figure 7 shows the effect of different fractional order parameter on plaque growth.The concentration of macrophages in the vessel wall increases with the increase of α, which leads the accumulation of foam cells and plaque formation.Therefore, α can be regarded as a parameter to describe the plaque growth rate, a key index in the study of atherosclerosis.

Conclusions and remarks
In this paper, we study a fluid-structure interaction problem with periodic forcing and temporal multiscale features (fast flow/slow plaque growth) and the slow variable equation contains fractional derivatives where a memory effect of the plaque evolution is included.We use the simple front tracking method to deal with the boundary growth and then formulate an auxiliary time periodic flow equation.Under the front tracking framework and for the linear Stokes flow, we analyze the error between solutions of the auxiliary periodic problem and the original problem.The error estimate is also expected to hold for the full Navier-Stokes equation.Based on this auxiliary time periodic problem, we then design an efficient multiscale algorithm and implement it in a combination of COMSOL Multiphysics and MATLAB.An iterative procedure to solve the timeperiodic problem is designed and its exponential convergence is analyzed.An error analysis for a fractional order time-discrete scheme of the multiscale algorithm is also provided.Several numerical experiments are conducted to illustrate the accuracy and efficiency of the proposed multiscale algorithm.The test results for simplified ODE systems and coupled Navier-Stokes systems show the great performance of the algorithm and that refining the macro-scale time step size can reduce the error more significantly than refining the micro-scale time step size.The final numerical example of the plaque growth problem shows the effect of different fractional order parameter on the plaque growth, and illustrates that the fractional order parameter may be used as an alternative index to reflect the plaque growth rate.
For future work, we shall consider to include the reaction-diffusion equation into the plaque growth model (PDE/PDE system) and explore the possibility of designing and analysing an effective multiscale method.

Figure 1 :
Figure 1: Schematic diagram of atherosclerotic with a plaque growth.

Figure 2 :
Figure 2: Schematic diagram of the multiscale method.Stepping the macro variable U (t) with macro step size ∆T 1, and then computing the periodic solution v U (t) using the micro step size ∆t over a unit cycle for a fixed flow domain Ω(U ).

Remark 4 .
It is possible to obtain the same results for a nonhomogeneous Dirichlet boundary condition since it can be transformed to the homogeneous Dirichlet boundary condition with appropriate assumptions on the smoothness of the domain and the regularity of the boundary value function.

Figure 3 :
Figure 3: Comparison of the numerical solution and exact solution with t = 1/100 and T = 1000.

Figure 4 :Figure 5 :
Figure 4: Error result (left) and CPU time (right)of different macro step size T with fixed t = 1/20 s.

Figure 7 :
Figure 7: Effects of different fractional parameter α on plaque growth.

Table 1 :
Error of macro variable U at t = 14000 s, convergence rates and CPU time with fixed ∆t = 1/100 and ∆T ε