How to solve two coupled equations?

[CalebDmArcher]

Hi,
I need to solve two coupled equations at steady state:
eq1 is a thermal equation with solution variable 'temperature'; eq2 is an electrical equation with solution variable 'voltage'. The equations are:

eq1 = DiffusionTerm(coeff=Conductivity)
eq2 = (DiffusionTerm(coeff=Gamma) +
      source +
      (RobinCoeff * g).divergence -
      ImplicitSourceTerm(coeff=(RobinCoeff * a.dot(n)).divergence))

The temperature affects the Conductivity of eq1, the electrical equation (I also use Cond to test different cases of this Conductivity variable); the voltage gradient and the resistivity affect the heat source of eq2, the thermal equation.

I was trying to test if the coupled equation can be solved by having a for loop (50 iterations) with a while loop inside. The inner while loop is only for sweeping eq1 and eq2 and the outer for loop is mainly for updating the resistivity, conductivity, and heat source. (And I plan to use a while loop to replace the outer for loop with the stop condition that resistivity and heat source values converge.)

The problem I am having right now is that if I have Cond = 1/Resistivity to be the coefficient in eq1, the values of 'Cond' will change with the values of Resistivity, and eq1 cannot converge after 1 or 2 for loop iterations. Similar problems exist no matter if I use Resistivity as CellVariable or FaceVariable.

But if I have this conductivity be a constant (by setting Conductivity as a CellVariable or FaceVariable the first time I define Conductivity and its values will not update unless setValue() being used), eq1 will have solutions converge.

It is kind of weird that eq1 cannot have solutions converged if its coefficient changes as eq1 is a very simple equation, steady state with only an implicit diffusion term and the coefficient is not updated during sweeping.

Would really appreciate it if you could let me know which part I did wrong or any suggestions for how to solve this type of coupled equations.

The full code is here:
example12_08_24.zip

[CalebDmArcher]

Hi, I need to solve two coupled equations at steady state: eq1 is a thermal equation with solution variable 'temperature'; eq2 is an electrical equation with solution variable 'voltage'. The equations are:

eq1 = DiffusionTerm(coeff=Conductivity)
eq2 = (DiffusionTerm(coeff=Gamma) +
      source +
      (RobinCoeff * g).divergence -
      ImplicitSourceTerm(coeff=(RobinCoeff * a.dot(n)).divergence))

The temperature affects the Conductivity of eq1, the electrical equation (I also use Cond to test different cases of this Conductivity variable); the voltage gradient and the resistivity affect the heat source of eq2, the thermal equation.

I was trying to test if the coupled equation can be solved by having a for loop (50 iterations) with a while loop inside. The inner while loop is only for sweeping eq1 and eq2 and the outer for loop is mainly for updating the resistivity, conductivity, and heat source. (And I plan to use a while loop to replace the outer for loop with the stop condition that resistivity and heat source values converge.)

The problem I am having right now is that if I have Cond = 1/Resistivity to be the coefficient in eq1, the values of 'Cond' will change with the values of Resistivity, and eq1 cannot converge after 1 or 2 for loop iterations. Similar problems exist no matter if I use Resistivity as CellVariable or FaceVariable.

But if I have this conductivity be a constant (by setting Conductivity as a CellVariable or FaceVariable the first time I define Conductivity and its values will not update unless setValue() being used), eq1 will have solutions converge.

It is kind of weird that eq1 cannot have solutions converged if its coefficient changes as eq1 is a very simple equation, steady state with only an implicit diffusion term and the coefficient is not updated during sweeping.

Would really appreciate it if you could let me know which part I did wrong or any suggestions for how to solve this type of coupled equations.

The full code is here: example12_08_24.zip

I have rearranged the equation-solving part of the code and seems like it works this way (the unable to converge problem did not show up):

res1 = 1e100
res2 = 1e100
while max(res1, res2) > 0.001:
    print("Before eq1")
    res1 = eq1.sweep(var=voltage, dt=1e-5)
    print(">>> eq1 is solved")

    source = (voltage.grad.dot(voltage.grad)) / Resistivity

    res2 = eq2.sweep(var=temperature, dt=1e-5)
    print(">>> eq2 is solved")

    Resistivity.setValue(
        Resistivity_0 * (1 + alpha_cu * (temperature - 20)) * ~mask_fr4_cell + Resistivity_0 * mask_fr4_cell)
    Cond = 1 / Resistivity
    gradient_x = -J_x * Resistivity_at_faces
    voltage.faceGrad.constrain(gradient_x * mesh.faceNormals, where=mask_top_cu_right_face)

In this code, source will update with voltage.grad and Resistivity; Resistivity will update itself by setValue(); Cond will update right after Resistivity updates; and other parameters that need to update with Resistivity will also get the update in the while loop (hopefully). This way I will not need to run the 50 iterations and the condition for when to accept the converged solution variables purely depends on the two residuals from the two equations sweeping.

But I still very much appreciate if you could give me any suggestions.

[guyer]

Can you please post the math for the equations and boundary conditions you are trying to solve?

[CalebDmArcher]

Hi Guyer, it is so great to hear from you.

This is eq1 which is the governing equation for voltage distribution in steady-state. In this equation, 𝜙 is the voltage and 𝜌 is the resistivity (So, 1/𝜌 is the conductivity).

$$ \left( \frac{\partial^2 \phi}{\partial x^2} + \frac{\partial^2 \phi}{\partial y^2} + \frac{\partial^2 \phi}{\partial z^2} \right) \frac{1}{\rho} = 0 $$

Below is eq2 which is the governing heat equation in steady-state. T is the temperature, k is the thermal conductivity, q_dot is the self-heat-generation term.

$$ \left( \frac{\partial^2 T}{\partial x^2} + \frac{\partial^2 T}{\partial y^2} + \frac{\partial^2 T}{\partial z^2} \right) k + \dot{q} = 0 $$

The coupling parameters between eq1 and rq2 are:

1. The electrical resistivity, 𝜌, needs to change with temperature by the expression below. rho_0, alpha, T_0 are constants.

$$ \rho = \rho_0 \left[ 1 + \alpha (T - T_0) \right] $$

I have two materials in this simulation, Cu and FR4. Only the resistivity of Cu needs to change with the temperature. (FR4 is an insulator so I set a very high constant resistivity for it rho_fr4 = 1.e+14)

2. The self-heat generation term source in eq2 needs to update with voltage.grad, the gradient of 𝜙 in eq1, and Resistivity, 𝜌 in eq1, by the Joule heating expression (This P_Joule is indeed source in my code):

$$ P_{\text{Joule}}(x, y, z) = \frac{\lvert \nabla \phi(x, y, z) \rvert^2}{\rho(x, y, z, T)} $$

Boundary conditions:

1. The boundary condition for eq1, the electrical equation, are one side of the top Cu layer has a contact 25 V voltage and a constant current goes through the top Cu layer from the 25 V side to the opposite side.
   
   The gradient of the voltage on the surface (opposite to the 25 V surface) can be expressed as the formula:

$$ \lvert \nabla \phi(x, y, z) \rvert = -\left( \frac{I}{S} \right) \rho(x, y, z, T) $$

In this expression, I is the current and S is the area of that surface. Since the electro resistivity, 𝜌, varies with temperature, this Neumann boundary condition also needs to be updated through iterations.
There is no other specific boundary conditions for the other surfaces.

2. The boundary condition for eq2, the thermal equation, is the top and the bottom surfaces in the image have free air convection (the surface index is not the same as shown in the graph below. eg. the top surface in the image below is actually mesh.facesBack). So the Robin boundary condition is used for the top and the bottom surface. Other surfaces are adiabatic so no specific boundary conditions are being coded.
   

Therefore, in the while loop of equation solving, other than the sweeping of two equations, there are also Resistivity, source, Cond (This is the conductivity in eq1), voltage.faceGrad.constain() (for the current go through the top Cu layer) need to be updated.

In my previous post, I forgot updateOld(), but I actually cannot spot any differences in the result with it. The updated while loop with updateOld() would be

res1 = 1e100
res2 = 1e100
while max(res1, res2) > 0.001:
    voltage.updateOld()
    temperature.updateOld()

    # print("Before eq1")
    res1 = eq1.sweep(var=voltage, dt=1e-5)
    print(">>> eq1 is solved in the current iteration")

    source = (voltage.grad.dot(voltage.grad)) / Resistivity

    res2 = eq2.sweep(var=temperature, dt=1e-5)
    print(">>> eq2 is solved in the current iteration")

    Resistivity.setValue(
        Resistivity_0 * (1 + alpha_cu * (temperature - 20)) * ~mask_fr4_cell + Resistivity_0 * mask_fr4_cell)
    Cond = 1 / Resistivity
    gradient_x = -J_x * Resistivity_at_faces
    voltage.faceGrad.constrain(gradient_x * mesh.faceNormals, where=mask_top_cu_right_face)

[guyer]

This is eq1 which is the governing equation for voltage distribution in steady-state. In this equation, 𝜙 is the voltage and 𝜌 is the resistivity (So, 1/𝜌 is the conductivity).

( ∂ 2 ϕ ∂ x 2 + ∂ 2 ϕ ∂ y 2 + ∂ 2 ϕ ∂ z 2 ) 1 ρ = 0

Because $\rho = \rho(T(x, y, z))$, this is not a correct representation of the continuity equation. Rather,

$$ \nabla\cdot\left(\frac{1}{\rho}\nabla \phi\right) = 0$$

which is what DiffusionTerm(coeff=Cond) represents.

Below is eq2 which is the governing heat equation in steady-state. T is the temperature, k is the thermal conductivity, q_dot is the self-heat-generation term.

( ∂ 2 T ∂ x 2 + ∂ 2 T ∂ y 2 + ∂ 2 T ∂ z 2 ) k + q ˙ = 0

Because $k$ is constant, this is fine... actually, no it's not. You assign different values to $k$ / $D$ / $\Gamma$ (pick one!) for different points in space, so more correct to write

$$\nabla\cdot\left(k\nabla T\right) + q = 0$$

Issues I see:

• eq1 is defined in terms of Cond, which is defined by Conductivity, which is assigned the value 1/rho_cu_0 or 1/rho_fr4. Calling Resistivity.setValue() in the loop doesn't change Conductivity or Cond. Define Conductivity in terms of Resistivity and let lazy evaluation do the work:

diff --git a/fipy_example43.py b/fipy_example43.py
index 6393590..d3c9b36 100644
--- a/fipy_example43.py
+++ b/fipy_example43.py
@@ -126,18 +126,17 @@ mask_top_cu_left_face = left_face & (z >= (t_cu + t_fr4))
 voltage.constrain(25., where=mask_top_cu_left_face)
 
 # Set conductance
-Conductivity = CellVariable(mesh=mesh, value=1/rho_cu_0)
 mask_fr4_face = (z > t_cu) & (z < (t_cu + t_fr4))
 mask_fr4_cell = (zc > t_cu) & (zc < (t_cu + t_fr4))
-conductance_fr4 = 1/rho_fr4
-Conductivity.setValue(conductance_fr4, where=mask_fr4_cell)
 
 # Set resistivity
 Resistivity_0 = CellVariable(mesh=mesh, value=rho_cu_0)
 Resistivity_0.setValue(rho_fr4, where=mask_fr4_cell)
 
-Resistivity = CellVariable(mesh=mesh, value=Resistivity_0)
+# Update Resistivity based on the new temperature
+Resistivity = Resistivity_0 * (1 + alpha_cu * (temperature - 20)) * ~mask_fr4_cell + Resistivity_0 * mask_fr4_cell
 
+Conductivity = 1 / Resistivity
 
 # Set a mask for the current flow going out from the right face of the top Cu layer
 mask_top_cu_right_face = right_face & (z >= (t_cu + t_fr4))
@@ -208,15 +207,9 @@ for t in range(50):
     voltage.updateOld()
     temperature.updateOld()
 
-    # Update Resistivity based on the new temperature
-    Resistivity.setValue(
-        Resistivity_0 * (1 + alpha_cu * (temperature - 20)) * ~mask_fr4_cell + Resistivity_0 * mask_fr4_cell)
     # Cond = 1 / Resistivity
     # Cond = 1 / (Resistivity.arithmeticFaceValue)

• You solve eq1 twice in succession, once with dt=1 and then with dt=1e-5:

diff --git a/fipy_example43.py b/fipy_example43.py
index e6094a1..b4706ec 100644
--- a/fipy_example43.py
+++ b/fipy_example43.py
@@ -190,7 +190,6 @@ for t in range(50):
     res2 = 1e100
     while max(res1, res2) > 0.1:
         print("Before eq1")
-        eq1.solve(var=voltage)
         res1 = eq1.sweep(var=voltage, dt=1e-5)
         print(">>> eq1 is solved")
         res2 = eq2.sweep(var=temperature, dt=1e-5)

• Defining source inside the loop doesn't do anything because this just creates a new variable called source that is not the same source used to define eq2. source is already a variable expression defined in terms of voltage.grad.

diff --git a/fipy_example43.py b/fipy_example43.py
index 9875941..a518a79 100644
--- a/fipy_example43.py
+++ b/fipy_example43.py
@@ -184,7 +184,6 @@ for t in range(50):
 
     # source = (voltage.grad ** 2) / Resistivity
     # source.setValue(voltage.grad.dot(voltage.grad) / Resistivity)
-    source = (voltage.grad.dot(voltage.grad)) / Resistivity
 
     res1 = 1e100
     res2 = 1e100

Therefore, in the while loop of equation solving, other than the sweeping of two equations, there are also Resistivity, source, Cond (This is the conductivity in eq1), voltage.faceGrad.constain() (for the current go through the top Cu layer) need to be updated.

Do not keep calling .constrain() in a loop. This does not replace the constraint; it just adds to a list. At best, the performance will get worse and worse, but the results may not work right, either.

In my previous post, I forgot updateOld(), but I actually cannot spot any differences in the result with it. The updated while loop with updateOld() would be

You have no TransientTerm in your equations, so updateOld() doesn't do anything.

Further recommendations:

• Particularly when troubleshooting, don't use such extreme differences in coefficients. E.g., rho_fr4 = rho_cu_0 * 1e10 seems to give stable solutions. Anything bigger looks questionable.
• Steady-state solutions often don't work well. Consider relaxing the solution to steady-state by adding a TransientTerm to each equation and taking large time steps.
• Consider actually coupling the equations, e.g. eq = eq1 & eq2. Sometimes convergence is better; sometimes it's worse.
• I strongly recommend troubleshooting in 2D. You have no pertinent 3D physics in this problem right now, it'll solve much faster, and it's much easier to visualize what's going on, e.g.,

mesh = Grid2D(dx=[dx] * nx,
              dy=[dz_cu] * nt_cu + [dz_fr4] * nt_fr4 + [dz_cu] * nt_cu)

[CalebDmArcher]

Thank you for the response.
After making the following modifications (no transient term yet):

#######################
# Electrical Parameters
# Define the electrical resistivity (constant for simplicity)
rho_cu_0 = 1.724e-8  # Resistivity of Copper (Ohm meter)\
alpha_cu = 0.00393
T_0 = 20.0
rho_fr4 = 1.e+12
# rho_fr4 = rho_cu_0*10

# Calculate the cross-sectional area in the yz-plane
S = w * t_cu

# Apply the voltage boundary condition
# Set 25 V on the left face of the top Cu layer
mask_top_cu_left_face = left_face & (z >= (t_cu + t_fr4))
voltage.constrain(25., where=mask_top_cu_left_face)

# Set resistivity
mask_fr4_face = (z > t_cu) & (z < (t_cu + t_fr4))
mask_fr4_cell = (zc > t_cu) & (zc < (t_cu + t_fr4))
# Set original resistivity
Resistivity_0 = CellVariable(mesh=mesh, value=rho_cu_0)
Resistivity_0.setValue(rho_fr4, where=mask_fr4_cell)
# Set updated resistivity
Resistivity = Resistivity_0 * (1 + alpha_cu * (temperature - 20)) * ~mask_fr4_cell + Resistivity_0 * mask_fr4_cell

# Set conductivity
Conductivity = 1 / (Resistivity)

# Set a mask for the current flow going out from the right face of the top Cu layer
mask_top_cu_right_face = right_face & (z >= (t_cu + t_fr4))

# Define the desired current density vector
current = 80.
J_x = current / S
Resistivity_at_faces = Resistivity.arithmeticFaceValue
gradient_x = -J_x * Resistivity_at_faces

# Define the boundary condition due to the current
voltage.faceGrad.constrain(gradient_x * mesh.faceNormals, where=mask_top_cu_right_face)

# Define equation 1
eq1 = DiffusionTerm(coeff=Conductivity)

#######################
# Thermal Parameters
# Define thermal conductivities
k_cu = 400.
k_fr4 = 0.8

# Define free air convection coefficient
h = 5.

# Define heat source
source = (voltage.grad.dot(voltage.grad)) / Resistivity

# Side surfaces (adiabatic)
# Top & bottom surfaces are free air convection: Robin boundary condition
free_air_surfaces = mesh.facesFront | mesh.facesBack

# Define the thermal conductivities for the equation
D = FaceVariable(mesh=mesh, value=k_cu)
D.setValue(k_fr4, where=mask_fr4_face)

Gamma = D * ~free_air_surfaces

dPf = FaceVariable(mesh=mesh,
                   value=mesh._faceToCellDistanceRatio * mesh.cellDistanceVectors)

n = mesh.faceNormals
a = FaceVariable(mesh=mesh, value=h*n, rank=1)
b = D
g = FaceVariable(mesh=mesh, value= h*T_ambient, rank=0)

RobinCoeff = free_air_surfaces * D * n / (dPf.dot(a) + b)

eq2 = (DiffusionTerm(coeff=Gamma) +
      source +
      (RobinCoeff * g).divergence -
      ImplicitSourceTerm(coeff=(RobinCoeff * a.dot(n)).divergence))

res1 = 1e100
res2 = 1e100
# coupledEqn = eq1 & eq2
# coupledEqn.solve(dt=..., solver=...)

while max(res1, res2) > 0.1:
    print("Before eq1")
    res1 = eq1.sweep(var=voltage, dt=1e-4)
    print(">>> eq1 is solved in the current iteration")
    res2 = eq2.sweep(var=temperature, dt=1e-4)
    print(">>> eq2 is solved in the current iteration")
    print(">>> Residual is ", max(res1, res2))

It seems like after making rho_fr4 smaller, eq1 will be able to converge. However, there are two problems:

1. The current physical model is designed to test if the code would work with more complex 3D problems. Therefore, the ‘3D’ problem in this simulation can be considered a 1D problem (current only flows from the 25 V surface to the opposite surface and heat only flows vertically); hence, the solution variables can be calculated by iterating Newton’s method (let’s call it theoretical values). FR4 is an insulator with very high impedance, the results from a low rho_fr4 value are always much lower than the theoretical results (as I increase rho_fr4, the final results get closer to the theoretical results.)

2. If set rho_fr4 smaller than what I set previously but still large enough, like 1e+12, the max residual from sweepings will be very large like in the million level. Here comes the problem of how to decide on a proper residual value.

But since when rho_fr4 is too large (1e+14), it leads to eq1 cannot converge, and in the ideal case, the conductance of the insulator is 0. So I tried to have the FR4 layer to have 0 conductance and let the high resistance only used for heat generation source. To do so, only make change with Conductivity: Conductivity = (1 / Resistivity)* ~mask_fr4_cell + (0.) * mask_fr4_cell. Surprisingly, this worked pretty well and the results are very close to the theoretical results without any convergence issue or huge residual issue.

I also tried to replace the equation-sweeping while loop into

coupledEqn = eq1 & eq2
coupledEqn.solve()
# coupledEqn.solve(dt=..., solver=...)

But get the following error:

ERROR:fipy:Uncaught exception:
Traceback (most recent call last):
  File "/mnt/d/.../thermal_model_full_package_2d_01052022/fipy_example45.py", line 197, in <module>
    coupledEqn = eq1 & eq2
  File "/usr/local/lib/python3.10/dist-packages/fipy/terms/term.py", line 487, in __and__
    return _CoupledBinaryTerm(self, other)
  File "/usr/local/lib/python3.10/dist-packages/fipy/terms/coupledBinaryTerm.py", line 39, in __init__
    raise SolutionVariableNumberError
fipy.terms.SolutionVariableNumberError: Different number of solution variables and equations.

I've tried both Conductivity = 1 / (Resistivity) and Conductivity = 1 / (Resistivity.arithmeticFaceValue) but both showed the above error.

I am not very sure how to transfer this 3D problem into a 2D problem with Fipy, as the current and the heat flow directions are different (current flows from left to right, and the heat flows vertically). Were you saying that I should make some modifications to the model itself and try to design it in a way that the current and the heat flow in the same 2D surface?

I will try to add transient terms to both equations to see if I can get better results. It is just that I saw some posts previously, for adding transient terms, I will need to adjust the diffusion coefficients. I haven't figured out how to do the modification yet. Will let you know if I achieve anything.

Again, really appreciate your help. It means a lot to me.

I have put the code for the theoretical calculation by using Newton's method, which might not be helpful, but in case if you are interested.
theoretical_calculation.zip

[guyer]

It seems like after making rho_fr4 smaller, eq1 will be able to converge. However, there are two problems:

1. The current physical model is designed to test if the code would work with more complex 3D problems. Therefore, the ‘3D’ problem in this simulation can be considered a 1D problem (current only flows from the 25 V surface to the opposite surface and heat only flows vertically); hence, the solution variables can be calculated by iterating Newton’s method (let’s call it theoretical values).

If the problem can be considered a 1D problem[*], then you should be trying to get it working in 1D first.

[*] I don't see how, if heat flux is orthogonal to current, but regardless, it's not intrinsically a 3D problem.

FR4 is an insulator with very high impedance, the results from a low rho_fr4 value are always much lower than the theoretical results (as I increase rho_fr4, the final results get closer to the theoretical results.)

Many physical problems cannot be solved numerically or, if they can, require very specialized techniques. Get the numerical problem working and then worry about "physical" coefficients. The model may not be suitable for what you want to do at all, or you may be able to gain insight even though the coefficients aren't "correct". Remember, "The purpose of computing is insight, not numbers". If your problem really depends on the resistivity varying by 22 orders of magnitude, then everything about this is probably wrong; wrong mesh, wrong equations, wrong tool.

2. If set rho_fr4 smaller than what I set previously but still large enough, like 1e+12, the max residual from sweepings will be very large like in the million level. Here comes the problem of how to decide on a proper residual value.

The magnitude of the residual depends on magnitude of the coefficients in your equations and on which solver suite and solver you are using[*]. You need to examine the residuals for your problem, confirm that they are decreasing with sweeps, and then pick a value that is "small enough".

To do so, only make change with Conductivity: Conductivity = (1 / Resistivity)* ~mask_fr4_cell + (0.) * mask_fr4_cell. Surprisingly, this worked pretty well and the results are very close to the theoretical results without any convergence issue or huge residual issue.

OK? That doesn't sound like a problem.

I also tried to replace the equation-sweeping while loop into

coupledEqn = eq1 & eq2
coupledEqn.solve()
# coupledEqn.solve(dt=..., solver=...)

But get the following error:

ERROR:fipy:Uncaught exception:
Traceback (most recent call last):
  File "/mnt/d/.../thermal_model_full_package_2d_01052022/fipy_example45.py", line 197, in <module>
    coupledEqn = eq1 & eq2
  File "/usr/local/lib/python3.10/dist-packages/fipy/terms/term.py", line 487, in __and__
    return _CoupledBinaryTerm(self, other)
  File "/usr/local/lib/python3.10/dist-packages/fipy/terms/coupledBinaryTerm.py", line 39, in __init__
    raise SolutionVariableNumberError
fipy.terms.SolutionVariableNumberError: Different number of solution variables and equations.

As it says in the link I gave previously, "for coupled equations the variable that a given term operates on now needs to be specified when the equation is generated"

I am not very sure how to transfer this 3D problem into a 2D problem with Fipy, as the current and the heat flow directions are different (current flows from left to right, and the heat flows vertically).

Right. That's two directions. Map (x, y, z) -> (x, z).

Were you saying that I should make some modifications to the model itself and try to design it in a way that the current and the heat flow in the same 2D surface?

Replace the mesh with what I showed previously:

mesh = Grid2D(dx=[dx] * nx,
              dy=[dz_cu] * nt_cu + [dz_fr4] * nt_fr4 + [dz_cu] * nt_cu)

You will get errors if you try to unpack three dimensions from cell and face centers; change these to unpack two:

-xc, yc, zc = mesh.cellCenters
-x, y, z = mesh.faceCenters
+xc, zc = mesh.cellCenters
+x, z = mesh.faceCenters

You will get errors about mesh.facesFront and mesh.facesBack. Replace these with mesh.facesBottom and mesh.facesTop.

I will try to add transient terms to both equations to see if I can get better results. It is just that I saw some posts previously, for adding transient terms, I will need to adjust the diffusion coefficients.

What?