Solving Multiphysics Problems
Walter Frei December 16, 2013
Here we introduce the two classes of algorithms used to solve multiphysics finite element problems in COMSOL Multiphysics. So far, we’ve learned how to mesh and solve linear and nonlinear single physics finite element problems, but have not yet considered what happens when there are multiple different interdependent physics being solved within the same domain.
A Simple SteadyState Multiphysics Problem
Let’s start by considering a very simple steadystate multiphysics problem: A coupling of steadystate electric current flow through a metal busbar, heat transfer in the bar, and structural deformations. Resistive heating arises due to the current flow, which raises the temperature of the bar and causes it to expand. In addition, the temperature rise will be significant enough that the electrical, thermal, and structural material property variations with temperature must be considered. We want to find the current flow, temperature fields, and deformations and stresses under steadystate conditions. The figure below shows a schematic of the problem being solved.
The multiphysics problem at hand.
The Associated Equations
There are here three governing partial differential equations being solved. First off, the equation that describes the voltage distribution within the domain is:
After discretizing via the finite element method, we can write a system of equations as:
where the subscript, _{V}, denotes the voltage unknowns, and the system matrix, \mathbf{K}_V, is dependent upon the temperature unknowns, \mathbf{u}_T. Assuming that the voltage distribution is known, then the volumetric resistive heating can be computed from:
where \bf{E}, the electric field, is: \nabla V. This heat source shows up in the governing equation for temperature:
And this equation gives us the system of equations:
Once we have the temperature distribution within the domain, we can solve for the structural displacements:
where the elasticity matrix, \bf{C}, is computed based on the temperaturedependent Young’s Modulus, E(T). The imposed thermal strain is \epsilon_{\Delta T} = \alpha(TT_0) and the strain is \epsilon = 1/2 [{\nabla \mathbf{u}^\mathbf{T}_D + \nabla \mathbf{u}_D}]. The system of equations that solve for the displacements is written as:
where the subscript, {_D}, indicates the displacement unknowns.
We can combine these systems of equations together:
\left\{
\begin{array}{c} \mathbf{f}_V \\ \mathbf{f}_T \\ \mathbf{f}_d \end{array}
\right\}
=
\left[
\begin{array}{ccc}
\mathbf{K}_V(\mathbf{u}_T) %26 \mathbf{0} %26 \mathbf{0} \\
\mathbf{0} %26 \mathbf{K}_T(\mathbf{u}_T) %26 \mathbf{0} \\
\mathbf{0} %26 \mathbf{0} %26 \mathbf{K}_D(\mathbf{u}_T)
\end{array}
\right]
\left\{
\begin{array}{c} \mathbf{u}_V \\ \mathbf{u}_T \\ \mathbf{u}_D \end{array}
\right\}

\left\{
\begin{array}{c} \mathbf{b}_V \\ \mathbf{b}_T(\mathbf{u}_V) \\ \mathbf{b}_D(\mathbf{u}_T) \end{array}
\right\}
\]
We can see by examination that this is a nonlinear problem, and as we learned earlier, this requires that we find the solution by taking NewtonRaphson iterations until we get convergence:
That is really all there is to it! There is no conceptual difference at all between solving a single physics nonlinear problem and solving a coupled physics problem. Everything that we have already learned about solving nonlinear single physics problems, including all of the discussions about damping, load and nonlinearity ramping, as well as meshing, is just as valid for solving a multiphysics problem.
Fully Coupled Approach
But it is also important to understand a (sometimes very serious) drawback to the above approach. During the NewtonRaphson iteration, we need to evaluate the derivative, \mathbf{f’(u}^{i}), so let’s write that out:
=
\left[
\begin{array}{ccc}
\mathbf{ K}_V%26\mathbf{K}_{V,T}%26\mathbf{0} \\
\mathbf{b}_{T,V}%26\mathbf{ K}_T%26\mathbf{0} \\
\mathbf{0}%26(\mathbf{K}_{D,T}\mathbf{b}_{D,T})%26\mathbf{K}_D
\end{array}
\right]
where the comma derivative notation is used, e.g.: \mathbf{K}_{V,T}=\partial \mathbf{K} _{V}(\mathbf{u}_T)/\partial \mathbf{u}_T.
Clearly, the above matrix is nonsymmetric, and this can lead to a problem: If the system matrix is not definite, then we may need to use the more memoryintensive direct solvers. (Although iterative solvers, with the right choice of preconditioner, can solve a wider class of problems they cannot be guaranteed to handle all cases.) Solving such a multiphysics problem with a direct solver will be both memory and timeintensive.
However, there is an alternative. The above method, called a Fully Coupled approach, assumes that all of the couplings between the physics must be considered at the same time. In fact, for the purposes of solving many types of multiphysics problems, we can neglect these offdiagonal terms during the solution, and solve using a more memory and timeefficient Segregated approach.
Segregated Approach
The Segregated approach treats each physics sequentially, using the results of the previously solved physics to evaluate the loads and material properties for the next physics being solved. So, using the above example, we first take a NewtonRaphson iteration for the voltage solution:
where, for the first iteration, we must have a starting guess for voltage and temperature ( \mathbf{u}_V^{i=0} , \mathbf{u}_T^{i=0} ). The material properties are evaluated using the initial conditions given for the temperature field. Next, the temperature solution is evaluated:
where, for the first iteration, i=0, the initial conditions given for the temperature field are used to evaluate the materials properties, \mathbf{K}_T(\mathbf{u}_T^{i=0}) , but the loads are evaluated based upon the voltage solution that was just computed: \mathbf{b}_T(\mathbf{u}_V^{i=1}) . Similarly, the displacement field is solved:
where the material properties and loads for the structural problem are computed using the temperature field computed above.
These iterations are then continued: voltage, temperature, and displacement are repeatedly computed in sequence. The algorithm is continued until convergence is achieved, as defined earlier.
The great advantage to the above approach is that the optimal iterative solver can be used in each linear substep. Not only are you now solving a smaller problem in each substep, but you can also use a solver that is more memoryefficient and generally solves faster. Although the segregated approach generally does require more iterations until convergence, each iteration takes significantly less time than one iteration of the fully coupled approach.
The algorithm used by the segregated solver for a model composed of n number of different physics is:
 Choose initial conditions for all physics in the model
 Initialize a counter for the number of iterations
 Solve for the first physics in the segregated sequence, using the previous step to evaluate material properties
 Solve for the second physics, using the part of the solution that has been computed to this point
 …
 Solve for the n^{th} physics, using the (n 1)^{th} previously computed parts of the solutions
 Repeat 26 until convergence, or until the desired peak number of iterations is exceeded
For general multiphysics problems, you will still have to choose the order in which the physics are solved, but the software has default suggestions as to an appropriate sequence for all builtin multiphysics interfaces. COMSOL Multiphysics will provide default linear solver settings for each physics in the segregated sequence.
When the segregated approach is applicable, it will converge to the same answer as the fully coupled approach. The segregated approach will usually take more iterations to converge; however, the memory and time requirements for each substep will be lower, so the total solution time and memory usage can be lower with the segregated approach.
Summary of Solving Multiphysics Problems
In this blog post, we have outlined the two classes of algorithms used to solve multiphysics problems — the Fully Coupled and the Segregated approach. The Fully Coupled approach is essentially identical to the NewtonRaphson method already developed for solving single physics nonlinear problems. It was shown to be very memoryintensive, but is useful, and generally needed, for multiphysics problems that have very strong interactions between the various physics being solved. On the other hand, the Segregated approach assumes that each physics can be solved independently, and will iterate through the various physics in the model until convergence.
Comments
Dear Frei,
Thanks for your instructions about the core functionality of multiphysics problems solving, and I have a question about the equation in “Fully Coupled Approach” about f(ui). I believe that the f(ui) is derived from “combined system equations”, but the Kv,t in the first row, and the Kd,t in the thrid row are not easy to be understood. If the Kv,t and Kd,t are remained in the matrix, then it may be not consistent with the “combined system equations”.
Wish to receive your reply.
Best wishes!
Dear Lingling, It is precisely because COMSOL does keep all of these terms in the system matrix that the fully coupled approach can be so powerful. Although it can be tricky to see how these affect the system, you can try to solve a strongly coupled problem (such as free convective flow) with a segregated approach and observe the differences. Of course, the memory requirements of the fully coupled approach can be quite a bit larger, so there is a cost to the benefit.
Dear Walter
I am going to debug my solver and find the best configuration for automatic selected solver. Would you please introduce me a basic source or assist me to find a start point.
Dar Maziyar, Keep in mind that COMSOL will always generate default solver suggestions which are usually appropriate. To understand when and why to override these defaults, keep reading the rest of this series of blogs on solvers!
Dear Walter,
I want to solve two coupled Governing Equations (Heat Transfer and Fluid Diffusion), and I am using COMSOL version 5.2. I want to choose a best solver. In the Heat Transfer Equation there is a spatial derivative of Pressure in the Convection Coefficient, and there is a time derivative of the Temperature in the Source Term of Fluid Diffusion Equation. Would you mind help me with this process?
Best Regards