1 Introduction↩︎

A number of studies have been reported in the literature on loosely coupled schemes for incompressible fluid-structure interaction (FSI) with thick-walled solids (see, e.g., [1]–[10]). In most of these works, either no rigorous error analysis is carried out or error estimates are provided which are not uniform with respect to the spatial discretization. Robin-Robin loosely coupled schemes, first proposed in [3] for FSI problems, were recently analyzed in [11]–[16]. Both the parabolic-parabolic and the hyperbolic-parabolic couplings have been considered in [14]. The FSI case corresponds to the hyperbolic-parabolic coupling with the additional difficulty coming from the fluid incompressibility constraint [11], [13], [15], [16]. The method is shown to be unconditionally stable and, for the first time, with a splitting error that is spatially uniform. Sub-optimal convergence rates in time of order \(\mathcal{O}(\sqrt{\Delta t})\) were first proved in [11]–[13]. However, this was later improved to nearly first-order accuracy in [14]–[16].

By leveraging the good stability properties and spatially uniform splitting error of the method in [13], [15], our aim is to construct a higher-order scheme in time by building on the prediction-correction methodology [17], with the Robin-Robin loosely coupled method as the prediction step. Although our ultimate goal is to develop a higher-order method for FSI problems, in this paper we start by constructing and analyzing a method for a coupled parabolic-parabolic problem. Indeed, the parabolic-parabolic problem is a classical model problem for fluid-fluid interaction, but it also involves some of the difficulties present in linear FSI problems, making it a good first step. It should be mentioned that correction methods have been used widely for FSI problems; see for example, [1]–[3], [13], [15], [18]. Burman and Fernández [1] already argued that one should expect a \(\Delta t\) improvement for every correction step. However, there does not seem to be an analysis of this for FSI problems (or for parabolic-parabolic problems). A method based on subcycling (several correction steps) was proposed in [19], [20], with linear convergence in the number of steps. In a domain decomposition framework, focusing on a multi-timestep approach, a Robin-Robin coupling for time-dependent advection–diffusion was introduced in [21] with numerical investigation of the stability.

In this paper, we aim to provide rigorous evidence of this by proving second-order convergence of a prediction-correction method with a single correction step, under an assumption of the convergence of the prediction step. The guiding principle of our correction method is to have discrete time differences of the prediction step error appearing in the right-hand side of the correction step error equations. Then, one uses the fact that the discrete time differences converge with second-order accuracy. We should mention that this strategy has also been used in [22] for a non-linear fluid problem, and the idea goes back to the work of Stetter and Böhmer et al. [17], [23], [24]. Typically, in the previous work on defect-correction for interface problems [22], the physical coupling conditions in fluid-fluid interaction introduced some dissipation, which made it possible to control the splitting error. This is often the case for fluid-fluid interactions modelling atmosphere-ocean interaction, where the two systems are coupled through the drag force on the ocean surface [25]–[32]. In our case, no such dissipative mechanism is present, and therefore the arguments of [22] fail. The conservative coupling presents an important challenge for the analysis, both for the low-order splitting scheme and the defect-correction scheme. Indeed, to the best of our knowledge, this is the first time second-order convergence is proved for a splitting scheme where the coupling conditions are conservative. A key milestone towards achieving second-order convergence of the correction step was reached in an accompanying paper [33], where we prove that many of the time difference errors of the prediction step are indeed second-order accurate in the \(L^2\)-norm. When the interface is perpendicular to the two sides of the domain, we also prove that the \(H^2\)-norm of the time difference errors of the prediction step also converges as \(\mathcal{O}((\Delta t)^2)\). This \(H^2\) estimate turns out to be a crucial ingredient in the analysis of the correction method. Although it is nontrivial to prove the \(H^2\) estimate in a more general setting, numerical results indicate that the \(H^2\) estimate also holds for more general cases.

Let us summarize the steps of the prediction-correction method at time \(t_{n+1}\):

• Solve the prediction step to obtain the defect solutions \(u_0^{n+1}\), \(w_0^{n+1}\) and \(\lambda_0^{n+1}\).

• Modify the right-hand side of the prediction step using the solutions \(u_0^{n+1}\), \(w_0^{n+1}\) and \(\lambda_0^{n+1}\), while keeping the left-hand side of the prediction step unchanged.

• The modification of the right-hand sides is designed in such a way that the right-hand sides of the error equations of the correction step consists of the term \(U_0^{n+1}-U_0^{n}\), \(W_0^{n+1}-W_0^{n}\), \(\Lambda_0^{n+1}-\Lambda_0^{n}\) and their gradients.

• Solve the new system with the modified right-hand side to obtain the correction solutions \(u_1^{n+1}\), \(w_1^{n+1}\) and \(\lambda_1^{n+1}\).

Here, \(\lambda_0\) denotes the Lagrange multiplier for the prediction step. The terms \(U_0^n\), \(W_0^n\) and \(\Lambda_0^n\) are the corresponding error terms of the variables \(u_0^n\), \(w_0^n\) and \(\lambda_0^n\) at time \(t_n\). As mentioned above, it has been observed numerically that \(U_0^{n+1}-U_0^{n}\), \(W_0^{n+1}-W_0^{n}\), and \(\Lambda_0^{n+1}-\Lambda_0^{n}\) converge with rate \(\mathcal{O}((\Delta t)^2)\) in the \(L^2\)-norm in general. Moreover, we prove that the first two errors are second-order accurate in the accompanying paper [33]. We also prove that the \(H^2\)-norm of \(U_0^{n+1}-U_0^{n}\) is second-order accurate for a special case. Note that this implies the second-order convergence of the Lagrange multiplier in \(L^2\)-norm. In Assumption 1 below, we assume the second-order accuracy of this error in order to carry out the analysis of the correction step. With this assumption we are able to prove second-order accuracy of the correction step. It should be noted that we perform the analysis for the time semi-discrete method. Extension to the fully discrete case are non-trivial and will be considered in the future.

The rest of the paper is organized as follows. In Section 2, we introduce a parabolic-parabolic interface problem and the corresponding prediction-correction method. In Section 4, we show that the correction method gives higher-order convergence in time under a crucial lemma when the interface is perpendicular to the two sides of \(\Omega\). In Section 5, we present two numerical examples that illustrate the theoretical findings. In Section 6, we consider the interface problem with Dirichlet boundary conditions and discuss an inconsistency that arises and briefly describe how to remedy it. We conclude the paper by drawing some future directions in Section 7.

2 Prediction-correction method for a parabolic\parabolic interface problem↩︎

Let \(\Omega\) be a bounded domain that can be decomposed as \(\Omega=\Omega_f \cup \Omega_s\) where the common interface \(\Sigma=\partial \Omega_f \cap \partial \Omega_s\). The subset \(\Gamma_D^f\) (resp., \(\Gamma_D^s\)) of the boundary \(\partial\Omega_f\) (resp., \(\partial\Omega_s\)) is the bottom side (resp., top side) of the domain \(\Omega_f\) (resp., \(\Omega_s\)). We then denote \(\Gamma_N^f=\partial\Omega_f\setminus(\Gamma_D^f\cap\Sigma)\) and \(\Gamma_N^s=\partial\Omega_s\setminus(\Gamma_D^s\cap\Sigma)\). See Figure 1 for an illustration.

We consider the following parabolic\parabolic interface problem. \[\tag{1} \begin{alignat}{2} \partial_{t}\mathsf{u}-\nu_f \Delta \mathsf{u}=&0, \quad && \text{ in } [0, T] \times \Omega_f, \nonumber\\ \mathsf{u}(0, x)=& \mathsf{u}_0(x), \quad && \text{ on } \Omega_f, \tag{2}\\ \mathsf{u}=&0, \quad && \text{ on } [0, T] \times \Gamma^f_D,\nonumber\\ \frac{\partial\mathsf{u}}{\partial n}=&0, \quad && \text{ on } [0, T] \times \Gamma^f_N\nonumber\\ \nonumber\\ \partial_{t}\mathsf{w}-\nu_s \Delta \mathsf{w}=&0, \quad && \text{ in } [0, T] \times \Omega_s, \nonumber\\ \mathsf{w}(0, x)=& \mathsf{w}_0(x), \quad && \text{ on } \Omega_s, \tag{3}\\ \mathsf{w}=&0, \quad && \text{ on } [0, T] \times \Gamma^s_D,\nonumber\\ \frac{\partial\mathsf{w}}{\partial n}=&0, \quad && \text{ on } [0, T] \times \Gamma^s_N\nonumber\\ \nonumber\\ \mathsf{w}- \mathsf{u}=& 0, \quad && \text{ in } [0,T] \times \Sigma, \\ \nu_s \nabla \mathsf{w}\cdot \boldsymbol{n}_s + \nu_f \nabla \mathsf{u}\cdot \boldsymbol{n}_f =& 0, \quad && \text{ in } [0,T] \times \Sigma, \end{alignat}\] where \(\boldsymbol{n}=\boldsymbol{n}_f\) and \(\boldsymbol{n}_s\) are the outward facing normal vectors of \(\Sigma\) for \(\Omega_f\) and \(\Omega_s\), respectively. Here \(\nu_f\) and \(\nu_s\) are piece-wise constant functions. We also assume that the initial data is smooth and that \(\mathsf{u}\) and \(\mathsf{w}\) is smooth on \(\Omega_f\) and \(\Omega_s\) respectively.

Note that we consider Neumann boundary conditions on two sides of the domain instead of Dirichlet boundary conditions, this is different from our previous work [14]. There are two main reasons we consider the problem 1 . First, the parabolic\parabolic interface problem 1 is more applicable in real-world scenarios, as numerous fluid-structure interaction problems involve a specific section within a broader domain where the fluid flows from left to right. Consequently, imposing Dirichlet boundary conditions to two vertical sides of the domain is unsuitable. Secondly, the numerical method for the Dirichlet problem has some inconsistencies which presents some difficulties in the analysis. We discuss this in Section 6.

2.1 Variational form↩︎

Let \((\cdot, \cdot)_i\) be the \(L^2\)-inner product on \(\Omega_i\) for \(i=f, s\). Moreover, let \(\bigl\langle\cdot, \cdot \bigr\rangle\) be the \(L^2\)-inner product on \(\Sigma\). Consider the spaces \[\label{eq:disspaces} \begin{alignat}{1} V_f:=&\{ v \in H^1(\Omega_f): v=0 \text{ on } \Gamma^f_D \}, \\ V_s:= &\{ v \in H^1(\Omega_s): v=0 \text{ on } \Gamma^s_D \}, \\ V_g:= & L^2(\Sigma). \end{alignat}\tag{4}\] We define the discrete time derivative as \[\begin{align} {1} \partial_{\Delta t}v^{n+1}:=\frac{v^{n+1}-v^n}{\Delta t}, \end{align}\] and the discrete average in time \[\begin{align} {1} v^{n+\frac{1}{2}}:=\frac{v^{n+1}+v^n}{2}. \end{align}\]

Furthermore, define the time step \(\Delta t\) and an integer \(N\), such that \(T = N\Delta t\), and let \(\mathsf{u}^{n} = \mathsf{u}(t_n, \cdot)\), where \(t_n = n\Delta t\) for \(n \in \{0, 1, 2, \dots, N\}\). Define \(\mathsf{l}=\nu_f \nabla \mathsf{u}\cdot \boldsymbol{n}_f\) and assuming that \(\mathsf{l}\in L^2(\Sigma)\) we get that the solution to 1 solves, for all \(t\in [0,T]\):

\[\begin{align} {2} (\partial_t \mathsf{u}, v)_f+\nu_f(\nabla \mathsf{u}, \nabla v)_f- \bigl\langle\mathsf{l}, v\bigr\rangle=&0, \quad &&\forall v \in V_f,\tag{5}\\ (\partial_t \mathsf{w},q )_s+\nu_s (\nabla \mathsf{w}, \nabla q)_s+ \bigl\langle\mathsf{l}, q\bigr\rangle=&0 , \quad && \forall q \in V_s, \tag{6}\\ \bigl\langle\mathsf{u}-\mathsf{w}, \mu \bigr\rangle=&0, \quad && \forall\mu \in V_g. \end{align}\]

2.2 Defect-correction method↩︎

We propose the following defect-correction method. The method consists of two steps, namely, the prediction step and the correction step. The prediction step is a Robin-Robin coupling method [14]. The correction step is designed by modifying the right-hand side of the prediction step while keeping the left-hand side of the prediction step unchanged as we discussed in the introduction.

2.2.0.1 Prediction step

Find \(w_0^{n+1} \in V_s\), \(u_0^{n+1} \in V_f\), and \(\lambda_0^{n+1} \in V_g\) such that for \(n\geq 0\), \[\tag{7} \begin{alignat}{2} (\partial_{\Delta t}w_0^{n+1}, z)_s+\nu_s(\nabla w_0^{n+1}, \nabla z)_s + \alpha \bigl\langle w_0^{n+1}-u_{0}^{n}, z \bigr\rangle+ \bigl\langle\lambda_{0}^{n}, z\bigr\rangle=&0 , \quad && \forall z \in V_s, \tag{8}\\ (\partial_{\Delta t}u_0^{n+1}, v)_f+\nu_f(\nabla u_0^{n+1}, \nabla v)_f- \bigl\langle\lambda_0^{n+1}, v\bigr\rangle=&0, \quad &&\forall v \in V_f, \tag{9}\\ \bigl\langle(\lambda_0^{n+1}- \lambda_{0}^{n})+ \alpha (u_0^{n+1}-w_0^{n+1}) , \mu \bigr\rangle=& 0, \quad &&\forall \mu \in V_g. \tag{10} \end{alignat}\]

2.2.0.2 Correction step

Find \(w_1^{n+1} \in V_s\), \(u_1^{n+1} \in V_f\), and \(\lambda_1^{n+1} \in V_g\) such that for \(n\geq 0\), \[\tag{11} \begin{alignat}{2} &(\partial_{\Delta t}w_1^{n+1}, z)_s+\nu_s(\nabla w_1^{n+1}, \nabla z)_s + \alpha \bigl\langle w_1^{n+1}-u_{1}^{n}, z \bigr\rangle+ \bigl\langle\lambda_{1}^{n}, z\bigr\rangle\quad &&\tag{12}\\ &=\nu_s(\nabla w_0^{n+1}, \nabla z)_s + \alpha \bigl\langle w_0^{n+1}-w_{0}^{n}, z \bigr\rangle+ \bigl\langle\lambda_{0}^{n}, z\bigr\rangle-\nu_s(\nabla w_0^{n+\frac{1}{2}}, \nabla z)_s - \bigl\langle\lambda_{0}^{n+\frac{1}{2}}, z\bigr\rangle, \quad && \forall z \in V_s,\nonumber \\ &(\partial_{\Delta t}u_1^{n+1}, v)_f+\nu_f(\nabla u_1^{n+1}, \nabla v)_f- \bigl\langle\lambda_1^{n+1}, v\bigr\rangle\quad && \tag{13}\\ &=\nu_f(\nabla u_0^{n+1}, \nabla v)_f- \bigl\langle\lambda_0^{n+1},v\bigr\rangle-\nu_f(\nabla u_0^{n+\frac{1}{2}}, \nabla v)_f+ \bigl\langle\lambda_0^{n+\frac{1}{2}}, v \bigr\rangle, \quad && \forall v \in V_f, \nonumber\\ & \bigl\langle(\lambda_1^{n+1}-\lambda_{1}^{n})+ \alpha (u_1^{n+1}-w_1^{n+1}) , \mu \bigr\rangle= \bigl\langle\lambda_0^{n+1}-\lambda_{0}^{n}, \mu \bigr\rangle, \quad && \forall \mu \in V_g. \end{alignat}\] In practice, the splitting method should be implemented sequentially. Assuming the information from the previous time-step is known, we solve the prediction step 7 for \(w_0^{n+1}\), \(u_0^{n+1}\) and \(\lambda_0^{n+1}\). We then compute the right-hand side of the correction step 11 using \(w_0^{n+1}\), \(u_0^{n+1}\) and \(\lambda_0^{n+1}\) and solve for \(w_1^{n+1}\), \(u_1^{n+1}\) and \(\lambda_1^{n+1}\) using 11 . Moreover, by design, each problem 7 , 11 can be split into two sub-problems.

3 Estimates for the prediction step↩︎

The key to the analysis of the correction method is the error estimates of the time difference quantities of the prediction step 7 . The proof of these results are established in [33]. In this section, we summarize these estimates. We define the errors \[\label{eq:errorsdef0} \begin{align} W_0^n&=\mathsf{w}^n-w_0^n,\\ U_0^n&=\mathsf{u}^n-u_0^n,\\ \Lambda_0^n&=\mathsf{l}^n-\lambda_0^n. \end{align}\tag{14}\] In [33], we show the following estimates for the prediction step 7 . \[\label{eq:pdstepesti} \begin{align} &\max_{0\le n\le N-1}\|W_0^{n+1}-W_0^n\|^2_{L^2(\Omega_s)}+\max_{0\le n\le N-1}\|U_0^{n+1}-U_0^n\|^2_{L^2(\Omega_f)}\\ &+\Delta t\sum_{n=0}^{N-1}\nu_f\|\nabla(U_0^{n+1}-U_0^n)\|^2_{L^2(\Omega_f)}+\Delta t\sum_{n=0}^{N-1}\nu_s\|\nabla(W_0^{n+1}-W_0^n)\|^2_{L^2(\Omega_s)}\le C (\Delta t)^4 \mathcal{Y}, \end{align}\tag{15}\] Additionally, we obtain the following estimate in [33], \[\label{eq:sdesti} \Delta t\sum_{n=0}^{M-1}\nu_f \|\frac{\nabla (U_0^{n+1}-2U_0^{n}+U_0^{n-1})}{(\Delta t)^2}\|_{L^2(\Omega_f)}^2\le C(\Delta t)^2\mathfrak{Y}.\tag{16}\] Here \(\mathcal{Y}\) and \(\mathfrak{Y}\) are defined in Appendix 8. At last, the following relation is established in [33] as well, \[\label{eq:L0rela} \Lambda_0^{n+1}=\nu_f\nabla U_0^{n+1}\cdot\boldsymbol{n}_f \quad\text{on}\quad\Sigma.\tag{17}\]

We then make the following crucial assumption.

Due to Assumption 1, the following Corollary is valid.

4 Error analysis for the correction step↩︎

In this section we provide the error analysis of 11 under Assumption 1. As argued in [14], the prediction step is first-order accurate. In this section, we will prove that the correction method produces a second-order accurate method provided Assumption 1 is satisfied.

4.1 Error equations↩︎

To obtain error equations, we write the exact solution in a way that can be compared to the numerical method.

\[\tag{20} \begin{alignat}{2} &(\partial_{\Delta t}\mathsf{w}^{n+1}, z)_s+\nu_s (\nabla \mathsf{w}^{n+1}, \nabla z)_s + \alpha \bigl\langle\mathsf{w}^{n+1}-\mathsf{u}^{n}, z \bigr\rangle+ \bigl\langle\mathsf{l}^{n}, z\bigr\rangle\quad &&\tag{21}\\ &\quad=\nu_s(\nabla \mathsf{w}^{n+1}, \nabla z)_s + \alpha \bigl\langle\mathsf{w}^{n+1}-\mathsf{u}^{n}, z \bigr\rangle+ \bigl\langle\mathsf{l}^{n}, z\bigr\rangle\nonumber\\ &\quad\quad-\nu_s(\nabla \mathsf{w}^{n+\frac{1}{2}}, \nabla z)_s - \bigl\langle\mathsf{l}^{n+\frac{1}{2}}, z\bigr\rangle+ T_1(z), \quad && \forall z \in V_s, \nonumber\\ &(\partial_{\Delta t}\mathsf{u}^{n+1}, v)_f+\nu_f(\nabla \mathsf{u}^{n+1}, \nabla v)_f- \bigl\langle\mathsf{l}^{n+1}, v\bigr\rangle\quad && \tag{22}\\ &\quad=\nu_f(\nabla \mathsf{u}^{n+1}, \nabla v)_f- \bigl\langle\mathsf{l}^{n+1},v\bigr\rangle-\nu_f(\nabla \mathsf{u}^{n+\frac{1}{2}}, \nabla v)_f+ \bigl\langle\mathsf{l}^{n+\frac{1}{2}}, v \bigr\rangle+ T_2(v), \quad && \forall v \in V_f, \nonumber\\ & \bigl\langle(\mathsf{l}^{n+1}-\mathsf{l}^{n}) +\alpha(\mathsf{u}^{n+1}-\mathsf{w}^{n+1}) , \mu \bigr\rangle= \bigl\langle\mathsf{l}^{n+1}-\mathsf{l}^{n}, \mu \bigr\rangle, \quad && \forall \mu \in V_g, \end{alignat}\] where \[\tag{23} \begin{alignat}{1} T_1(z):= &(\partial_{\Delta t}\mathsf{w}^{n+1}, z)_s+\nu_s(\nabla \mathsf{w}^{n+\frac{1}{2}}, \nabla z)_s + \bigl\langle\mathsf{l}^{n+\frac{1}{2}}, z\bigr\rangle, \\ T_2(v):=& (\partial_{\Delta t}\mathsf{u}^{n+1}, v)_f+\nu_f(\nabla \mathsf{u}^{n+\frac{1}{2}}, \nabla v)_f - \bigl\langle\mathsf{l}^{n+\frac{1}{2}}, v\bigr\rangle.\tag{24} \end{alignat}\] Now we define the errors \[\label{eq:errorsdef1} \begin{align} W_1^n&=\mathsf{w}^n-w_1^n,\\ U_1^n&=\mathsf{u}^n-u_1^n,\\ \Lambda_1^n&=\mathsf{l}^n-\lambda_1^n. \end{align}\tag{25}\] Subtract 11 and 20 , we obtain the error equations for \(n\ge0\), \[\tag{26} \begin{alignat}{2} &(\partial_{\Delta t}W_1^{n+1}, z)_s+\nu_s(\nabla W_1^{n+1}, \nabla z)_s + \alpha \bigl\langle W_1^{n+1}-U_1^{n}, z \bigr\rangle+ \bigl\langle\Lambda_1^{n}, z\bigr\rangle=R_1(z) +T_1(z) \quad && \forall z \in V_s, \tag{27}\\ &(\partial_{\Delta t}U_1^{n+1}, v)_f+\nu_f(\nabla U_1^{n+1}, \nabla v)_f- \bigl\langle\Lambda_1^{n+1}, v\bigr\rangle= R_2(v)+T_2(v) \quad && \forall v \in V_f, \tag{28}\\ & \bigl\langle\Lambda_1^{n+1}-\Lambda_1^{n}+ \alpha(U_1^{n+1}-W_1^{n+1}), \mu \bigr\rangle= \bigl\langle\Lambda_0^{n+1}-\Lambda_0^{n}, \mu \bigr\rangle\quad && \forall\mu \in V_g.\tag{29} \end{alignat}\] Here \[\begin{alignat}{1} R_1(z) :=&\frac{1}{2}\nu_s(\nabla (W_0^{n+1}-W_0^n), \nabla z)_s + \alpha \bigl\langle(W_0^{n+1}-W_0^{n}), z \bigr\rangle- \frac{1}{2}\bigl\langle\Lambda_0^{n+1}-\Lambda_0^{n}, z\bigr\rangle, \\ R_2(v): =& \frac{1}{2} \nu_f(\nabla (U_0^{n+1}-U_0^n), \nabla v)_f - \frac{1}{2}\bigl\langle\Lambda_0^{n+1}-\Lambda_0^{n}, v\bigr\rangle.\label{eq:R2} \end{alignat}\tag{30}\] It follows from 29 that, \[\label{errorstrong} \Lambda_1^{n+1}-\Lambda_1^{n}+ \alpha ( U_1^{n+1}-W_1^{n+1})= \Lambda_0^{n+1}-\Lambda_0^{n}\quad\text{on}\quad\Sigma.\tag{31}\] We use the convention \(u_1^i=\mathsf{u}^i=\mathsf{u}_0\) and \(w_1^i=\mathsf{w}^i=\mathsf{w}_0\) for \(i<0\).

4.2 Preliminary estimate for the correction step↩︎

In this section, we establish some preliminary estimates that will be used in the convergence analysis. The techniques we use here is similar but slightly different to those in [33]. First, we define the following quantities for \(n\ge0\), \[\begin{align} {1} \mathsf{Z}_1^{n+1}:= & \frac{1}{2}\|W_1^{n+1}\|_{L^2(\Omega_s)}^2+ \frac{1}{2} \|U_1^{n+1}\|_{L^2(\Omega_f)}^2 +\frac{\Delta t\alpha}{2} \|U_1^{n+1}\|_{L^2(\Sigma)}^2 + \frac{\Delta t}{2\alpha} \|\Lambda_1^{n+1}\|_{L^2(\Sigma)}^2, \\ \mathsf{S}_1^{n+1}:= & \Delta t(\nu_f \|\nabla U_1^{n+1}\|_{L^2(\Omega_f)}^2+ \nu_s \|\nabla W_1^{n+1}\|_{L^2(\Omega_s)}^2) + \frac{1}{2} (\|W_1^{n+1}-W_1^n\|_{L^2(\Omega_s)}^2 + \|U_1^{n+1}-U_1^n\|_{L^2(\Omega_f)}^2)\\ &+ \frac{\alpha \Delta t}{2}\|(U_1^n-U_1^{n+1})+ \frac{1}{\alpha}(\Lambda_1^n-\Lambda_1^{n+1}) \|_{L^2(\Sigma)}^2. \end{align}\]

We will use the following identity where \((\cdot,\cdot)\) is an inner product and \(\|\cdot\|\) is the corresponding norm: \[\label{eq:innerid} (\varphi-\psi,\vartheta)=\frac{1}{2}(\|\varphi\|^2-\|\psi\|^2+\|\vartheta-\psi\|^2-\|\vartheta-\varphi\|^2).\tag{32}\] We start with the following lemma.

In the following subsections, we give preliminary estimates of all the terms on the right-hand side of 33 . To this purpose, we will need the following trace inequality, which follows from the continuity of the trace operator and from Poincaré’s inequality.

4.3 Estimate for \(F_1^{n+1}\)↩︎

4.4 Estimate for \(R^{n+1}\)↩︎

4.5 Estimate for \(T^{n+1}\)↩︎

4.6 Estimate for \(\frac{\Delta t}{\alpha}\bigl\langle\Lambda_1^{n+1}, \Lambda_0^{n+1}-\Lambda_0^{n} \bigr\rangle\)↩︎

In this section we estimate the term \[\label{eq:L1term} \frac{\Delta t}{\alpha}\sum_{n=0}^{M-1}\bigl\langle\Lambda_1^{n+1}, \Lambda_0^{n+1}-\Lambda_0^{n} \bigr\rangle.\tag{43}\] Note that \(\nabla U_0^{n+1}\cdot\boldsymbol{n}_f\) is well-defined on \(\Omega_f\). Hence we may extend \(\nabla U_0^{n+1}\cdot\boldsymbol{n}_f\) into the domain \(\Omega_f\) in a similar fashion as in [33]. Let \(\phi:\Omega_f\rightarrow \mathbb{R}\) be an affine function that is \(1\) on \(\Sigma\) and vanishes on \(\Gamma_D^f\). Therefore, the following properties are satisfied: \[\begin{align} &\phi\in V_f,\tag{44}\\ &0\le\phi\le 1,\tag{45}\\ &\phi=1\quad\text{on}\quad\Sigma,\tag{46}\\ &\|\nabla\phi\|_{L^2(\Omega_f)}\le C.\tag{47} \end{align}\] Define the function \[\tilde{\Lambda}_0^{n+1}:=\phi(\nu_f\nabla U_0^{n+1}\cdot\boldsymbol{n}_f).\] Note that \(\tilde{\Lambda}_0^{n+1}=\Lambda_0^{n+1}\) on \(\Sigma\) and \(\tilde{\Lambda}_0^{n+1}\in V_f\), so it is indeed an extension. It then follows that \[\bigl\langle\Lambda_1^{n+1},\Lambda_0^{n+1}-\Lambda_0^{n}\bigr\rangle=\bigl\langle\Lambda_1^{n+1},\tilde{\Lambda}_0^{n+1}-\tilde{\Lambda}_0^{n}\bigr\rangle.\] Denote \(\mathsf{L}^{n+1}:=\Lambda_0^{n+1}-\Lambda_0^{n}\) and \(\tilde{\mathsf{L}}^{n+1}:=\tilde{\Lambda}_0^{n+1}-\tilde{\Lambda}_0^{n}\). Furthermore, owing to 28 with \(v=\tilde{\mathsf{L}}^{n+1}\), we have \[\begin{align} \bigl\langle\Lambda_1^{n+1},\tilde{\mathsf{L}}^{n+1}\bigr\rangle&=(\partial_{\Delta t}U_1^{n+1}, \tilde{\mathsf{L}}^{n+1})_f+\nu_f(\nabla U_1^{n+1}, \nabla \tilde{\mathsf{L}}^{n+1})_f-R_2(\tilde{\mathsf{L}}^{n+1})-T_2(\tilde{\mathsf{L}}^{n+1}). \end{align}\] Therefore, we obtain \[\label{eq:L1ext} \begin{align} \frac{\Delta t}{\alpha}\sum_{n=0}^{M-1}\bigl\langle\Lambda_1^{n+1},\Lambda_0^{n+1}-\Lambda_0^{n}\bigr\rangle&=\frac{1}{\alpha}\sum_{n=0}^{M-1}( U_1^{n+1}-U_1^n, \tilde{\mathsf{L}}^{n+1})_f+\frac{\Delta t}{\alpha}\sum_{n=0}^{M-1}\nu_f(\nabla U_1^{n+1}, \nabla \tilde{\mathsf{L}}^{n+1})_f\\ &\quad-\frac{\Delta t}{\alpha}\sum_{n=0}^{M-1} R_2(\tilde{\mathsf{L}}^{n+1})-\frac{\Delta t}{\alpha} \sum_{n=0}^{M-1}T_2(\tilde{\mathsf{L}}^{n+1}). \end{align}\tag{48}\] The following lemma states the estimate for \(\eqref{eq:L1term}\), based on the previous identity.

4.7 Main result: error estimate↩︎

In this section, we state and prove our main result under Assumption 1 and Corollary 1. To this purpose, we first need the following lemma.

We now state and prove the main result of this paper.

5 Numerical experiments↩︎

In this section, we present two examples which numerically illustrate the analysis and observations of the previous sections. We use the following notations to denote errors of the correction step, \[\begin{align} {2} &e_{u1}:= \|u_1^N-\mathsf{u}^N\|_{L^2(\Omega_f)},\quad&& e_{w1}:= \|w_1^N-\mathsf{w}^N\|_{L^2(\Omega_s)},\\ &e_{du1}:= \|\nabla(u_1^N-\mathsf{u}^N)\|_{L^2(\Omega_f)}, \quad && e_{\lambda}:=\|\lambda_0^N-\mathsf{l}^N\|_{L^2(\Sigma)}, \\ & e_{1,\lambda}:=\|(\lambda_0^N-\mathsf{l}^N)-(\lambda_0^{N-1}-\mathsf{l}^{N-1})\|_{L^2(\Sigma)}. \end{align}\] All numerical experiments are performed using FEniCS and multiphenics [34], [35].

We can clearly see \(\mathcal{O}((\Delta t)^2)\) convergence rates for \(e_{u1}\), \(e_{w1}\) and \(e_{du}\), which agrees with our analysis. One also observes the optimal convergence rates for \(e_{1,\lambda}\) which match with Corollary 1.

6 The Dirichlet interface problem↩︎

In this section, we briefly discuss the defect-correction methods for the following parabolic\parabolic interface problem with Dirichlet boundary conditions on all sides (cf. [14]). As we shall demonstrate there is an inconsistency in the formulation at the the end points of \(\Sigma\) when \(\Sigma\) is not perpendicular to the two sides of \(\partial \Omega\). Then, this inconsistency makes the second-order rates of \(\Lambda_0^{n+1}- \Lambda_0^n\) deteriorate, which means that the Corollary 1 is not valid in this case. This in turn destroys the second-order convergence rate for the correction method measured in the \(H^1\) norm. Nonetheless, when measured in the \(L^2\) norm the approximations of the correction method seem to be second-order accurate.

The Dirichlet problem is as follows. \[\label{PPsplitF} \begin{alignat}{2} \partial_{t}\mathsf{u}- \nu_f\Delta \mathsf{u}=&0, \quad && \text{ in } (0, T) \times \Omega_f, \\ \mathsf{u}(0, x)=& \mathsf{u}_0(x), \quad && \text{ in } \Omega_f, \\ \mathsf{u}=&0 , \quad && \text{ on } (0, T) \times \partial \Omega_f \setminus \Sigma, \end{alignat}\tag{57}\]

\[\label{PPsplitS} \begin{alignat}{2} \partial_{t}\mathsf{w}-\nu_s \Delta \mathsf{w}=&0, \quad && \text{ in }(0, T) \times \Omega_s, \\ \mathsf{w}(0, x)=& \mathsf{w}_0(x), \quad && \text{ in } \Omega_s, \\ \mathsf{w}=&0, \quad && \text{ on } (0, T) \times \partial \Omega_s \setminus \Sigma, \end{alignat}\tag{58}\] and \[\label{PPsplitI} \begin{alignat}{2} \mathsf{w}- \mathsf{u}=& 0, \quad && \text{ in } (0, T) \times \Sigma, \\ \nu_s\nabla \mathsf{w}\cdot \boldsymbol{n}_s + \nu_f \nabla \mathsf{u}\cdot \boldsymbol{n}_f =& 0 \quad && \text{ in } (0, T) \times \Sigma, \end{alignat}\tag{59}\] where \(\boldsymbol{n}_f\) and \(\boldsymbol{n}_s\) are the outward facing normal vectors for \(\Omega_f\) and \(\Omega_s\), respectively.

We define the following spaces: \[\label{eq:disspacesdiri} \begin{alignat}{1} \tilde{V}_f=&\{ v \in H^1(\Omega_f): v=0 \text{ on } \partial\Omega_f\setminus\Sigma \}, \\ \tilde{V}_s= &\{ v \in H^1(\Omega_s): v=0 \text{ on } \partial\Omega_s\setminus\Sigma \}, \\ \tilde{V}_g= & L^2(\Sigma). \end{alignat}\tag{60}\]

Note that the variational form for 57 -59 is the same as that of 1 . The defect-correction method for 57 -59 is almost identical to that of 1 , the only difference is that we replace the \(V\) spaces with the \(\tilde{V}\) spaces.

There is an inconsistency of \(\lambda_0\) at the end points of \(\Sigma\). As in 10 we have \(\lambda_0^{n+1}-\lambda_0^n=\alpha(w_0^{n+1}-u_0^{n+1})\) on \(\Sigma\). Moreover, \(w_0^{n+1}\) and \(u_0^{n+1}\) have zero Dirichlet boundary conditions at two end points of \(\Sigma\) which force \(\lambda_0^{n+1}=\lambda_0^n\) at those endpoints. Therefore, \(\lambda_0^{n+1}\) does not change in time hence it does not approximate \(\mathsf{l}^{n+1}\) at those endpoints and consequently affect the convergence rates on \(\Sigma\).

We can clearly see that \(e_{1,\lambda}\) is not second-order accurate and hence Corollary 1 is not valid in this case. Moreover, the gradient error of \(\mathsf{u}\) of the correction method, \(e_{du1}\), is not second-order. Finally, we observe that \(e_{u1}\) and \(e_{w1}\) are second-order accurate. It is worthwhile mentioning that if the interface is perpendicular to two of the sides of \(\partial \Omega\) then there is no inconsistency for the Dirichlet problem and, numerically, one regains second-order convergence for all the variables above. We omit the numerical experiments here.

We finish this section by stating that there is a remedy for the Dirichlet problem, when \(\nu_s=\nu_f=\nu\), in the case \(\Sigma\) is straight line and that is not parallel to the two sides of \(\Omega\). To be specific, we denote \(\boldsymbol{\tau}\) as the tangent vector of \(\Sigma\) and \(\mathbf{s}\) as the tangent vector of the side of \(\Omega\). It is easy to see that there exists constants \(a\) and \(b\) such that \(\mathbf{n}_f=a\boldsymbol{\tau}+b\mathbf{s}\), and hence \(\nabla \mathsf{u}\cdot\boldsymbol{n}_f=a\nabla \mathsf{u}\cdot\boldsymbol{\tau}+b\nabla \mathsf{u}\cdot\mathbf{s}\). Note that \(\mathbf{n}_s=-a\boldsymbol{\tau}-b\mathbf{s}\).

We then define the following modified defect-correction method.

6.0.0.1 Prediction step

Find \(w_0^{n+1} \in \tilde{V}_s\), \(u_0^{n+1} \in \tilde{V}_f\), and \(\lambda_0^{n+1} \in \tilde{V}_g\) such that for \(n\geq 0\), \[\tag{61} \begin{alignat}{2} &(\partial_{\Delta t}w_0^{n+1}, z)_s+\nu(\nabla w_0^{n+1}, \nabla z)_s + b\alpha \bigl\langle w_0^{n+1}-u_{0}^{n}, z \bigr\rangle\tag{62}\\ &\quad+ a\nu\bigl\langle\nabla w^{n+1}_0\cdot\boldsymbol{\tau},z\bigr\rangle+\bigl\langle\lambda^n,z\bigr\rangle=0 , \quad && \forall z \in \tilde{V}_s, \nonumber\\ &(\partial_{\Delta t}u_0^{n+1}, v)_f+\nu(\nabla u_0^{n+1}, \nabla v)_f- \bigl\langle\lambda_{0}^{n+1}, v\bigr\rangle-a\nu\bigl\langle\nabla u^{n+1}_0\cdot\boldsymbol{\tau},v\bigr\rangle=0, \quad &&\forall v \in \tilde{V}_f, \tag{63}\\ &\bigl\langle(\lambda_0^{n+1}- \lambda_{0}^{n})+ \alpha (u_0^{n+1}-w_0^{n+1}) , \mu \bigr\rangle= 0, \quad &&\forall \mu \in \tilde{V}_g. \tag{64} \end{alignat}\]

6.0.0.2 Correction step

Find \(w_1^{n+1} \in \tilde{V}_s\), \(u_1^{n+1} \in \tilde{V}_f\), and \(\lambda_1^{n+1} \in \tilde{V}_g\) such that for \(n\geq 0\), \[\label{correction1m} \begin{alignat}{2} &(\partial_{\Delta t}w_1^{n+1}, z)_s+\nu(\nabla w_1^{n+1}, \nabla z)_s + b\alpha \bigl\langle w_1^{n+1}-u_{1}^{n}, z \bigr\rangle\nonumber \\ &\quad\quad+ a\nu\bigl\langle\nabla w^{n+1}_1\cdot\boldsymbol{\tau},z\bigr\rangle+\bigl\langle\lambda_{1}^{n}, z\bigr\rangle\quad &&\nonumber\\ &=\nu(\nabla w_0^{n+1}, \nabla z)_s + b\alpha \bigl\langle w_0^{n+1}-w_{0}^{n}, z \bigr\rangle+ \bigl\langle\lambda_{0}^{n}, z\bigr\rangle+a\nu\bigl\langle\nabla w^{n+1}_0\cdot\boldsymbol{\tau},z\bigr\rangle\\ &\quad\quad-\nu(\nabla w_0^{n+\frac{1}{2}}, \nabla z)_s - \bigl\langle\lambda_{0}^{n+\frac{1}{2}}, z\bigr\rangle-a\nu\bigl\langle\nabla w^{n+\frac{1}{2}}_0\cdot\boldsymbol{\tau},z\bigr\rangle, \quad && \forall z \in \tilde{V}_s,\nonumber \\ &(\partial_{\Delta t}u_1^{n+1}, v)_f+\nu(\nabla u_1^{n+1}, \nabla v)_f- \bigl\langle\lambda_1^{n+1}, v\bigr\rangle-a\nu\bigl\langle\nabla u^{n+1}_1\cdot\boldsymbol{\tau},v\bigr\rangle\quad && \nonumber\\ &=\nu(\nabla u_0^{n+1}, \nabla v)_f- \bigl\langle\lambda_0^{n+1},v\bigr\rangle-a\nu\bigl\langle\nabla u^{n+1}_0\cdot\boldsymbol{\tau},v\bigr\rangle-\nu(\nabla u_0^{n+\frac{1}{2}}, \nabla v)_f\\ &\quad\quad+ \bigl\langle\lambda_0^{n+\frac{1}{2}}, v \bigr\rangle+a\nu\bigl\langle\nabla u^{n+\frac{1}{2}}_0\cdot\boldsymbol{\tau},v\bigr\rangle, \quad && \forall v \in \tilde{V}_f, \nonumber\\ & \bigl\langle(\lambda_1^{n+1}-\lambda_{1}^{n})+ \alpha (u_1^{n+1}-w_1^{n+1}) , \mu \bigr\rangle= \bigl\langle\lambda_0^{n+1}-\lambda_{0}^{n}, \mu \bigr\rangle\quad && \forall \mu \in \tilde{V}_g. \end{alignat}\tag{65}\]

We then show the numerical results of the modified method 65 for Example 3 in Table 4. We can clearly see the improvement of \(e_{\lambda}\), \(e_{1,\lambda}\) and \(e_{du1}\) that match with our observation.

7 Future directions↩︎

There are several future directions. One, of course, is to prove Assumption 1 for problems with more general interface. This seems to be a challenging problem and new ideas are needed. The second is to extend the methodology to the wave-parabolic interface problem. Preliminary results are very promising. Once we accomplish this, we will be ready to tackle the linear FSI problem. Finally, we will consider the dynamic FSI problem.

Acknowledgement↩︎

This material is based upon work supported by the National Science Foundation under Grant No. DMS-1929284 while the last author was in residence at the Institute for Computational and Experimental Research in Mathematics in Providence, RI, during the "Numerical PDEs: Analysis, Algorithms, and Data Challenges" program. EB was supported by EPSRC grants EP/J002313/2 and EP/W007460/1. EB and MF were supported by the IMFIBIO Inria-UCL associated team.

8 Notations↩︎

For readability, we collect various definitions in Section 3 here. We define \(\mathcal{Y}\) as, \[\begin{align} {1} \mathcal{Y}:=& \frac{1}{\nu_s}\|\partial_{t}^3 \mathsf{w}\|_{L^2(0,T;L^2(\Omega_s))}^2 + (\frac{1}{\nu_f}+1)\|\partial_{t}^3 \mathsf{u}\|_{L^2(0,T;L^2(\Omega_f))}^2) \\ & + (\frac{\nu_f}{\alpha^2}+\frac{\nu_f^2}{\alpha}) \|\partial_{t}^2 \mathsf{u}\|_{L^2(0,T;H^1(\Omega_f))}^2+ \frac{(\nu_f)^3}{\alpha^2} \|\partial_{t}^2 \mathsf{u}\|_{L^2(0,T;H^2(\Omega_f))}^2\\ & +\frac{\alpha^2}{\nu_s}\|\partial_{t}^2 \mathsf{u}\|_{L^2(0,T;L^2(\Sigma))}^2 + \frac{1}{\nu_f}\|\partial_{t}^2 \mathsf{l}\|_{L^2(0,T;L^2(\Sigma))}^2+\frac{\nu_f^2}{\alpha^2} \|\partial_{t}^2 \mathsf{u}\|_{L^\infty(0,T;H^1(\Omega_f))}^2, \end{align}\] define \(\mathfrak{Y}\) as, \[\begin{align} {1} \mathfrak{Y}:=& \frac{1}{\nu_s}\|\partial_{t}^4 \mathsf{w}\|_{L^2(0,T;L^2(\Omega_s))}^2 + (\frac{1}{\nu_f}+1)\|\partial_{t}^4 \mathsf{u}\|_{L^2(0,T;L^2(\Omega_f))}^2) \\ & + (\frac{\nu_f}{\alpha^2}+\frac{\nu_f^2}{\alpha}) \|\partial_{t}^3 \mathsf{u}\|_{L^2(0,T;H^1(\Omega_f))}^2+ \frac{(\nu_f)^3}{\alpha^2} \|\partial_{t}^3 \mathsf{u}\|_{L^2(0,T;H^2(\Omega_f))}^2\\ & +\frac{\alpha^2}{\nu_s}\|\partial_{t}^3 \mathsf{u}\|_{L^2(0,T;L^2(\Sigma))}^2 + \frac{1}{\nu_f}\|\partial_{t}^3 \mathsf{l}\|_{L^2(0,T;L^2(\Sigma))}^2+\frac{\nu_f^2}{\alpha^2} \|\partial_{t}^3 \mathsf{u}\|_{L^\infty(0,T;H^1(\Omega_f))}^2, \end{align}\] define \(Y\) as, \[\begin{align} Y:=& \frac{1}{\nu_s}\|{\partial}_x\partial_{t}^2 \mathsf{w}\|_{L^2(0,T;L^2(\Omega_s))}^2 + (\frac{1}{\nu_f}+1)\|{\partial}_x\partial_{t}^2 \mathsf{u}\|_{L^2(0,T;L^2(\Omega_f))}^2) \\ & + (\frac{\nu_f}{\alpha^2}+\frac{\nu_f^2}{\alpha}) \|{\partial}_x\partial_{t}\mathsf{u}\|_{L^2(0,T;H^1(\Omega_f))}^2+ \frac{(\nu_f)^3}{\alpha^2} \|{\partial}_x\partial_{t}\mathsf{u}\|_{L^2(0,T;H^2(\Omega_f))}^2\\ & +\frac{\alpha^2}{\nu_s}\|{\partial}_x\partial_{t}\mathsf{u}\|_{L^2(0,T;L^2(\Sigma))}^2 + \frac{1}{\nu_f}\|{\partial}_x\partial_{t}\mathsf{l}\|_{L^2(0,T;L^2(\Sigma))}^2+\frac{\nu_f^2}{\alpha^2} \|{\partial}_x\partial_{t}\mathsf{u}\|_{L^\infty(0,T;H^1(\Omega_f))}^2, \end{align}\] and define \(\mathbb{Y}\) as, \[\label{eq:bbY} \begin{align} \mathbb{Y}:=& \frac{1}{\nu_s}\|{\partial}_x\partial_{t}^3 \mathsf{w}\|_{L^2(0,T;L^2(\Omega_s))}^2 + (\frac{1}{\nu_f}+1)\|{\partial}_x\partial_{t}^3 \mathsf{u}\|_{L^2(0,T;L^2(\Omega_f))}^2) \\ & + (\frac{\nu_f}{\alpha^2}+\frac{\nu_f^2}{\alpha}) \|{\partial}_x\partial_{t}^2 \mathsf{u}\|_{L^2(0,T;H^1(\Omega_f))}^2+ \frac{(\nu_f)^3}{\alpha^2} \|{\partial}_x\partial_{t}^2 \mathsf{u}\|_{L^2(0,T;H^2(\Omega_f))}^2\\ & +\frac{\alpha^2}{\nu_s}\|{\partial}_x\partial_{t}^2 \mathsf{u}\|_{L^2(0,T;L^2(\Sigma))}^2 + \frac{1}{\nu_f}\|{\partial}_x\partial_{t}^2 \mathsf{l}\|_{L^2(0,T;L^2(\Sigma))}^2+\frac{\nu_f^2}{\alpha^2} \|{\partial}_x\partial_{t}^2 \mathsf{u}\|_{L^\infty(0,T;H^1(\Omega_f))}^2. \end{align}\tag{66}\]

References↩︎