Multi-domain, multi-modeling problems find important applications in science and engineering. They usually involve coupled processes or systems in more than one geometric domains or physical fields, and arise from classical parallel computing, domain decomposition, to contemporary multi-scale and multi-physics computation. Mathematically, they are modeled by different partial differential equations(PDEs) in different domains with certain interface coupling conditions, such as fluid-structure interactions (FSI)[9, 16, 23], coupling fluid flow with porous media flow [1, 7, 10], etc..
Numerical methods for multi-domain, multi-modeling problems are generally classified as monolithic and partitioned approaches. The monolithic approach usually leads to coupled schemes[3, 20, 21, 24]. Although unconditionally stable and convergent in general, those coupled schemes usually result in significant difficulties and inflexibility in the design and choice of mesh generation, PDE discretization, algebraic solvers, as well as software development.
On the other hand, certain partitioned approaches, often called loosely coupled, explicit coupling, or decoupled schemes, have also been investigated in the literature [2, 4, 5, 6, 11, 12, 13, 17, 18], where the sub-model in each domain is solved locally so that existing legacy solvers may be applied by easy software integration. There are many other important physical and numerical considerations that appeal to partitioned approaches for treating different sub-models in their own physical regions independently, particularly because of the very different physical properties and time scales. However, decoupled approaches are usually theoretically much more difficult. Decoupling could easily lead to numerical instability or non-convergence. For instance, many decoupled methods in FSI applications are so-called unconditionally instable due to the artificial added-mass effect [8, 14]. Even with classical non-overlapping domain decomposition, for a long time people have been still searching for and trying to understand how to decouple the computation effectively to ensure stability and convergence. The most typical and conventional approach is to match the given interface conditions, for instance the continuity of the Dirichlet and Neumann data, in one way or another by solving local sub-models with certain given boundary conditions and then exchanging the information across the interface in certain ad hoc or intuitive way, which often leads to decoupled methods of the Dirichlet-Dirichlet or Dirichlet-Neumann type. These methods are sometimes stable and convergent, sometime not, or not very efficient. A Robin type condition has been introduced by combining the original Dirichlet and Neumann interface conditions to obtain a set of mathematically equivalent interface conditions, which results in some other types of decoupled methods such as Dirichlet-Robin, Neumann-Robin, and Robin-Robin type of schemes. Sometimes, this may help to improve the effectiveness and efficiency. However, such an approach is artificial and lack of physical justification, which results in both theoretical and practical numerical difficulties for the Robin type decoupling approach. For instance, in applications like FSI the Dirichlet data correspond to velocity, while the Neumann data correspond to interfacial force, therefore, there is a mismatch between these two different types of physical quantities. In the classical Robin approach, this physical mismatch is somehow meredied by mathematically introducing a relaxation parameter, which is to be tuned for the purposes of numerical stability and convergence. The optimal relaxation parameter is then compensated for the mathematical effect of the physical mismatch. However, the optimal parameter analysis is often very difficult, even for certain simple model problems. The choice and tuning of the relaxation parameter with the existing numerical methods are practically very difficult and even impossible for most of the real applications, never mentioning the theoretical analysis and understanding of the numerical stability and convergence.
In an interesting FSI application of a thin-wall structure coupled with a fluid flow, a so-call intrinsic Robin condition was introduced with certain formal derivation, which leads to a stable decoupled scheme and overcomes the added-mass instability effect [12, 13]. Motivated by this, we believe that this could provide certain insight and help understand the interface mechanism, as well as provide guidance for developing decoupling methodology for general multi-domain and multi-physics applications.
We investigate a typical multi-modeling model problem. By cruising various stages of numerical approximation and decoupling, and track how the information is transmitted across the interface, we clearly demonstrate how the interface conditions evolve and affect the data and error transmission across the interface during spacial discretization, time discretization, as well as decoupling. An approximate intrinsic or inertial type Robin condition is rigorously derived for a semi-discrete model with finite element spacial discretization, which makes sense both mathematically and physically in contrast to the classical Robin condition. Based on this new interface condition, an equivalent semi-discrete model is introduced, which provides a general framework for devising effective decoupled numerical methods. Numerical experiments also confirm the effectiveness of this new decoupling approach.
The paper is organized as follows. The multi-domain model problem is described in Section 2. The derivation and analysis of the new intrinsic or inertial type Robin interface condition are carried out in Section 3. The decoupling approach based on the intrinsic or inertial type Robin interface condition is introduced in Section 4, followed by numerical experiments in Section 5. Concluding remarks are given in Section 6.
2 Model Problem
We consider a scalar coupled model consisting of two heat equations in two adjacent subdomains coupled through interface conditions. Let and be the two subdomains, with the interface as showed in Figure 1, where and
are the unit outward normal vectors toand , respectively.
The coupled PDE model is described as:
where and are temperature, and are density, and are material diffusivity, and are external source, and the Dirichlet and Neumann conditions are imposed on the interface by (3) and (4). For illustration and without loss of generality, we will assume .
It serves as a typical model for investigating decoupling techniques for various important applications, bridging classical domain decomposition and parallel computing with contemporary multi-physics computing. When and , this coupled model is equivalent to a underline single model of a heat equation defined over the global domain, being partitioned to subdomain problems with non-overlapping domain decomposition so that various domain decomposition time marching schemes or iterative methods may be devised for parabolic or elliptic problems with parallel computing. When , and/or , it may be used as a model problem to study decoupling techniques for jumping coefficients and multi-scale applications. The scalar model problem may also be used to investigate and understand the numerical and physical properties of decoupling techniques for being further extended to other important applications in multi-modeling, multi-physics, multi-scale computation such as FSI.
Define the function spaces
Model C: The Weak Form of the Continuous Coupled Model
Find , such that
3 Semi-Discretization and Approximate Intrinsic or Inertial Type Robin Conditions
The weak formulation is mathematically equivalent to the strong form of the original coupled PDE model in the conventional sense, and its solution still satisfies the original Dirichlet and Neumann conditions (2.3) and (2.4) on the interface. There could be various equivalent interface conditions to be derived for various purposes. To understand what should be the appropriate interface conditions for a numerical scheme to satisfy for the sake of stability and convergence, in particular when we want to devise effective decoupling techniques later, let us start with the original continuous model and investigate the mechanism on how the interface conditions would change from the original model due to discrete approximations.
Note that the issue of numerical stability during time marching in numerical computation is caused by the discretization in both space and time. To understand how the spacial discretizatin affects the interface conditions, yet without the stability factor, let us concentrate on the finite element semi-discrete approximation:
Model : The Semi-Discrete Model with Finite Element Approximation in Space
Find , such that
where is a standard finite element subspace of . Note now that Model is no longer equivalent to Model C
However, it is essentially important to point out that the approximate solution of Model no longer satisfies the interface conditions of the original coupled model although the Dirichlet condition remains enforced in the solution space , to which little attention has been paid in the past, mainly because the fully implicit time discretization will then lead to an unconditionally stable scheme, and the convergence and error analysis are ready in the classical finite element analysis, so it appears unnecessary to further investigate what happens on the interface. However, if decoupling is necessary as a further approximation, it would then become crucially important to find out what is indeed the corresponding NEW interface coupling conditions implied by this approximate semi-discrete model due to the spacial discretization, which would then provide us insight for decoupling the numerical computation properly in order to keep maintaining numerical stability and convergence for decoupled methods.
Unfortunately, it is difficult to derive the local interface conditions for such a weak formulation defined in the integration form over the entire geometric domain. For finding out such an interface condition approximately, motivated by , we propose to apply the lumped-mass technique , which will lead to a further approximation to the semi-discrete model.
Let us now briefly review the lumped-mass technique. Let be a matrix and be the corresponding lumped-mass matrix
Denote consists of continuous, piecewise linear functions on a quasiuniform family of triangulations of any domain , where or . The inner product in is defined by
Let be a triangle of the triangulation and be its vertices. Define the quadrature formula for any function
then define the lumped-mass approximation of inner product by
Let be the finite element basis, the following two properties hold :
From the properties, if is generated by , then is generated by .
Denote be the corresponding discrete spaces with continuous piecewise linear functions. Applying the lumped-mass technique to , we may build connections among the inner products in and . Define two discrete lifting operators
such that the nodal values of and vanish out of for all . The interface operator is defined by , where is the adjoint operator of w.r.t the lumped-mass inner product in . is also defined similarly. Based on the discrete operators introduced above, we have
The following Lemmas are useful in the analysis below. For simplicity, we denote by .
For all , we have
Similarly, we also have
For all , we have
We may now define the lumped-mass semi-discrete model as a further approximation, denoted by Model .
Model : The Lumped-Mass Semi-Discrete Model
Find , such that
Under certain regularity assumptions, there are also standard error estimates for the solutions of the two semi-discrete models, as well as the original continuous model . In particular, if has regularity, for each subdomain we have
There are many important applications in numerical computation where the Dirichlet condition enforced in the solution space should be relaxed at certain steps for decoupling purposes. Such cases arise from matching the Dirichlet data on the interface by the classical non-overlapping domain decomposition approach for constructing parallel time marching schemes for parabolic problems or for constructing parallel iterative methods for elliptic problems by simulating the time evolution to the equilibrium state, to today’s decoupling for multi-physics computation. For understanding the interface mechanism in order to construct effective decoupled numerical methods as well as to conduct the corresponding numerical analysis in the future, it will be fundamentally important to examine the following modified problem to the above Model by removing the Dirichlet interface condition from the solution space , where the modification of is defined as .
Problem : A Perturbation to Model
Consider a function satisfying
The following theorem is revealing.
If satisfies (20) in Problem , then satisfies the following condition:
Given , taking and in (20), we have
The lumped-mass property (16) in Lemma 3.1 allows us to condense the inner product in subdomain to the interface :
In order to pass the Dirichlet information from to , denote , we have
Reorganizing the terms by separating the information on each subdomain yields
Note that the left hand side in (26) consists of two part: the expression in the bracket corresponds to the terms in (20) from Problem for subdomain , and the time derivative term being condensed to the interface from subdomain , but then exchanged the Dirichlet interface data with subdomain yet leaving an error term in the right hand side, which will play the key role in the intrinsic Robin condition to be derived.
For symmetry and later on to compare with the semi-discrete model, let us add and subtract the same term to equation (26), which gives
The fundamental significance of 1 as well as its derivation is in that the error equation (21) reveals the interface mechanism of the semi-discrete model, in particular the impact of the error in the Dirichlet data on the interface conditions when the Dirichlet data is transmitted across the interface, which results instead in the error of the time derivative of the Dirichlet data in the form of on the interface. Mathematically, we note that the effect of the Dirichlet interface error at a given time here leads to an interface error of in the time evolution model, which implies how the Dirichlet type of information should be properly related on different time levels when time discretization is applied, in particular when decoupling is necessary along the interface. It explains why in most of the classical explicit type of decoupled schemes, the explicit use of Dirichlet data from previous time levels, such as the Dirichlet-Dirichlet, Dirichlet-Neumann type of decoupling methods, would usually lead to stability issues. Most importantly, we note that when this scalar model problem is extended to other applications such as FSI, the Dirichlet variable will physically correspond to velocity, while the related time derivative interface terms in (21) will correspond to certain interface inertial force quantities due to spacial discretization, which will make the new interface condition to be derived physically meaningful in contrast to the classical Robin interface condition.
Notice that and thus become zero in the error equation (21) if the Dirichlet condition is enforced on the interface. Therefore, as a special case of (21), for (18) in Model we have the following interface relationship.
If is a solution of Model , then satisfies the following interface condition:
Mathematically, (28) is equivalent to (18) in Model since the additional two terms and are canceled each other due to the Dirichlet interface condition being enforced in the solution space. However, they could make significant differences numerically when further approximation is applied for time discretization as well as decoupling.
To further understand the interface mechanism implied by (28), let us assume all the regularity requirements, well-posedness, as well as error estimates for the related continuous and semi-discrete models in order to focus on illustrating the derivation and understanding the mechanism of the intrinsic Robin condition to be derived. Let us also first concentrate on the terms inside the two brackets in (28), denoted by
which actually corresponds to all the terms in (18) of the lumped-mass model. Because (18) in Model is an approximation to (9) in Model due to the lumped-mass condensation of the inner product to the interface, we may relate the terms in back to (9) in Model by approximating the lumped-mass terms using and with an error of , where the error between the inner product and its corresponding lumped-mass approximation in is as showed in . So, we have the estimate for with the solution of Model , yet in terms of the expression in Model as follows.
Now that Model is also the standard finite element approximation to Model C, while Model is an approximation to Model C, too. Therefore, by using (19) we may approximate in (30) by the continuous solution with the error of to further relate the expression to the weak form of the coupled model, yet restricted to the test function space :
Under sufficient regularity assumptions for the continuous solution and following the standard procedure, integrating by parts in (31) and using the strong form of the coupled PDE model leads to
Therefore, we see that actually corresponds to the Newmann condition in the original interface matching conditions, or physically presenting the jump of the flux across the interface.