M1: The numerical solution of coupled multiphysics problems

Project completed September 30, 2010

Background

The need to model systems involving multiple physical processes is a pressing one across the sciences. The mathematical modelling of many physical phenomena is composed of sub-problems that contain several unique physical phenomena. 

These multiphysics problems are usually modelled by systems of differential equations that depend on each other, which are generally too large for a direct numerical solution. It is often the case that there exist mature numerical methods to solve each of the sub-problems separately. Therefore it is desirable to decompose the system into a collection of smaller sub-systems and use these established numerical methods.

Researchers at the Oxford Centre for Collaborative Applied Mathematics (OCCAM) examined techniques for decoupling such multiphysics problems. 

Techniques and Challenges

The project is divided into finding techniques for dealing with two types of coupling: one-way coupling and weak coupling.

In one-way coupling with two sub-systems, the first sub-problem depends directly on the second sub-problem, but the second sub-problem does not depend at all on the first sub-problem. One-way coupling was explored using a model of the heart. 

In the heart, an electrical current is initiated at a node where it causes the ion concentrations within the heart tissue cells to change. This leads to tension in the tissue, causing it to contract and deform – a heart beat. The coupling is strong from the electrophysiology to the tissue deformation, but there is assumed to be no dependence in the other direction under normal circumstances.

The sub-problems were modelled separately. The electrophysiology was described using the Fithugh-Nagumo equations, and tissue deformation using linear elasticity. These are well-known models for these phenomena.

Each of the sub-problems was solved using the finite element method (FEM). With the FEM, the domain is split into elements, often either small triangles or quadrilaterals. The accuracy of the FEM depends on the size of the elements. The smaller the elements, the less error there is in the solution, but the smaller the elements are, the more complicated the solution is to compute, thus, there is a trade-off between element size and computation time. It is desirable to use the coarsest grid possible for a given error tolerance.

To identify where smaller elements are needed, a posteriori error estimation is used. Firstly, the error cannot be calculated simply from solving the main problem numerically. The solution to the dual problem must be found. The dual problem is similar, but not identical to the main problem and involves the quantity of interest. The numerical solution for the main and the dual problem are compared in each element to determine the error. If the error contribution is high, the element is broken up into smaller elements until the error is sufficiently small.

The second part of the project focused on the effect of ignoring weaker coupling in multiphysics problems and modelling the problem as one-way coupled.

An approach that is often used to solve two-way coupled systems is Gauss-Seidel iteration. Here, the first system is solved using a guess to the solution of the second system. Then, the determined solution for the first system is used to solve the second system. This is repeated until the solution is found within a given error tolerance.

The aim of this part of the project was to automatically determine when the coupling between the systems is weak in one direction in order to avoid using later Gauss-Seidel iterations. Using a combination of the computed solution, the quantity of interest and matrices representing the discretisations of the sub-systems, the error bound was determined. If the error is large, further Gauss-Seidel iterations can be performed until the error is small enough to assume one-way coupling.

Results

For one-way coupling, it was determined that the numerical error from solving the electrophysiology problem is much smaller than from solving the tissue deformation problem, meaning the electrophysiology solution can be computed on a coarser mesh.

In the second part of the project, a technique was developed that automatically determines when the coupling between sub-problems is weak, and importantly, the effect of ignoring the weak coupling can be quantified.

Overall, the project determined that it is possible to use a coarse mesh for the electrophysiology of the heart with little effect on the accuracy of the deformation solution assuming one-way coupling and using a posteriori analysis. This improves the efficiency of computing solutions and is promising for similar multiphysics problems.

The Future

A long-term aim of this project is to apply this methodology to more complex models of the heart. Since the electrophysiology and tissue deformation are only one-way coupled under normal circumstances, it is desirable to investigate the extreme circumstances and model the weak coupling from the tissue deformation to the electrophysiology, which would facilitate the study of heart attacks caused by impact or when a thump to the chest is used to restart a stopped heart.