How to solve the equation only on a portion of the whole structure?

[CalebDmArcher]

Hi, hope someone could help. This might be a very basic question but I am having difficulty finding the solution.

I am trying to simulate the temperature distribution of the following structure

It has 3 layers and only the top layer has a constant voltage and a constant current applied on it which generates heat. And I want to simulate the temperature distribution over the whole structure in the steady state.

The electrical equation is:

$$\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$$

ϕ is the voltage;
ρ is the resistivity.

Since the middle layer is an insulator, the electrical behavior only happens within the top layer.

Is there a way that let the solver only solve this equation on the top layer part (which might help with the calculation load)?
And if there is a way to do so, will this cause any trouble if I want to couple this equation with a thermal equation which needs to be calculated over the whole structure?

I have tried other two approaches which both does not work:
To set the current, I used the Robin boundary method:

# Set conductance
Conductances = FaceVariable(mesh=mesh, value=1/rho_cu)
# FR4 conductance is approximately 0
mask_fr4 = (z > t_cu) & (z < (t_cu + t_fr4))
# conductance_fr4 = 1/(1.e+12)
conductance_fr4 = 0
Conductances.setValue(conductance_fr4, where=mask_fr4)

# 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))

Gamma_ele = Conductances * ~mask_top_cu_right_face

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

# Define current source
current = 10.
I_S = current / S

n = mesh.faceNormals
a_ele = FaceVariable(mesh=mesh, value=0.*n, rank=1)
b_ele = Conductances
g_ele = FaceVariable(mesh=mesh, value=-I_S, rank=0)

RobinCoeff = mask_top_cu_right_face * Conductances * n / (dPf_ele.dot(a_ele) + b_ele)

eq = (DiffusionTerm(coeff=Gamma_ele) +
      (RobinCoeff * g_ele).divergence -
      ImplicitSourceTerm(coeff=(RobinCoeff * a_ele.dot(n)).divergence))

Method 1: set the coefficient of the middle part to 0, which triggers error

RuntimeWarning: invalid value encountered in divide
  return self._BinaryOperatorVariable(lambda a, b: a / b, other)

Method 2: Set the coefficient to a very small value (I tried 1/(1.e+12)).
This way no error was reported; however, the result shows that the middle layer does not behave like an insulator at all.
This is why I am trying to only calculate for the top layer and not touch other parts of the structure at all.

[guyer]

Is there a way that let the solver only solve this equation on the top layer part (which might help with the calculation load)?

Not really. You can solve a completely separate problem on a mesh that only describes the top layer and then transcribe the values to the full mesh for the thermal equation, but there's no way to turn the solve off for part of a mesh.

And if there is a way to do so, will this cause any trouble if I want to couple this equation with a thermal equation which needs to be calculated over the whole structure?

Transcribing values is prone to stability issues, but they're not insurmountable.

n = mesh.faceNormals
a_ele = FaceVariable(mesh=mesh, value=0.*n, rank=1)
b_ele = Conductances
g_ele = FaceVariable(mesh=mesh, value=-I_S, rank=0)

RobinCoeff = mask_top_cu_right_face * Conductances * n / (dPf_ele.dot(a_ele) + b_ele)

If a_ele is zero everywhere, then you don't have a Robin condition, you have a Neumann condition.

Method 1: set the coefficient of the middle part to 0, which triggers error

Set which coefficient to 0? RobinCoeff or Gamma_ele. Making Gamma_ele zero or very small should be OK. Setting RobinCoeff to zero doesn't make any sense to me.

Method 2: Set the coefficient to a very small value (I tried 1/(1.e+12)). This way no error was reported; however, the result shows that the middle layer does not behave like an insulator at all. This is why I am trying to only calculate for the top layer and not touch other parts of the structure at all.

Without seeing the result or a functional script, I don't know what to tell you

[CalebDmArcher]

insurmountable

Hi, thank you for replying to my previous post.

Here is the code:
example09_19_24.zip

By further checking on it, I think the code (currently I haven't coupled this equation with another equation yet) works now which has results close to the theoretical results.

But there are still some difficulties I am not sure what I can do:

The first one is to set part of the diffusion coefficient, conductance, to 0 for the insulator layer. (Currently, I am using a small value conductance_fr4 = 1/(1.e+12), which seems working fine) Although I guess I probably should not set this value to 0 because there is no absolute zero conductance material, I was trying to do so since I thought it might make the calculation load lighter.
The conductance variable is the 1/ρ term in the equation below.

$$\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$$

The error when set part of this variable to 0 is:

/usr/local/lib/python3.10/dist-packages/fipy/variables/variable.py:1141: RuntimeWarning: invalid value encountered in divide
  return self._BinaryOperatorVariable(lambda a, b: a / b, other)

I just realized that it might be better if I use the Neumann boundary condition instead of modifying the Robin boundary condition.

But if there is another case in which I want this current to flow through the entire top layer with a non-uniform current distribution (different from the example above that the current has the same value and the same direction through the entire layer), is there a way I can set this current flow?

[guyer]

The error when set part of this variable to 0 is:

/usr/local/lib/python3.10/dist-packages/fipy/variables/variable.py:1141: RuntimeWarning: invalid value encountered in divide
  return self._BinaryOperatorVariable(lambda a, b: a / b, other)

This is not an error; it's a warning. You should not ignore warnings, but you also don't need to panic about them.

I just realized that it might be better if I use the Neumann boundary condition instead of modifying the Robin boundary condition.

I agree. I don't think you have a Robin condition.

But if there is another case in which I want this current to flow through the entire top layer with a non-uniform current distribution (different from the example above that the current has the same value and the same direction through the entire layer), is there a way I can set this current flow?

I would rearrange Ohm's law to give a Neumann expression for the gradient in potential in terms of the desired current density: $\vec{i} = -\frac{1}{\rho}\nabla \phi \rightarrow \nabla\phi = -\rho\vec{i}$

[CalebDmArcher]

The error when set part of this variable to 0 is:

/usr/local/lib/python3.10/dist-packages/fipy/variables/variable.py:1141: RuntimeWarning: invalid value encountered in divide
  return self._BinaryOperatorVariable(lambda a, b: a / b, other)

This is not an error; it's a warning. You should not ignore warnings, but you also don't need to panic about them.

I just realized that it might be better if I use the Neumann boundary condition instead of modifying the Robin boundary condition.

I agree. I don't think you have a Robin condition.

But if there is another case in which I want this current to flow through the entire top layer with a non-uniform current distribution (different from the example above that the current has the same value and the same direction through the entire layer), is there a way I can set this current flow?

I would rearrange Ohm's law to give a Neumann expression for the gradient in potential in terms of the desired current density: i → = − 1 ρ ∇ ϕ → ∇ ϕ = − ρ i →

1. With the previous code, if set conductance_fr4 = 0., there will be the warning and the output will be 0 which is not correct.
   The comparison of the results between setting a very small value and 0 are shown below:
   

2. I am trying to set the gradient inside the shape and build the new equation with the following code (The change made starts from line 127 of the previous code).
   This is inspired by the example of Internal Fixed Gradient. And this method is trying to utilize the equation
   ∇ ϕ = − ρ i →

The problem with this code is that it shows a constant 25 V voltage across the entire top surface no matter how big the current I set (but overall it shows a high voltage on the top surface and low voltage on the bottom surface which seems current).

I suspect that the gradient is not correct especially the part that finding the corresponding faceNormals direction might be wrong. But I am having difficulty finding the definition of elements of faceNormals. Would really appreciate if you could point out where I could find it.

(I also attached the full code here just in case)
example09_25_24.zip

# Set conductance
Conductances = FaceVariable(mesh=mesh, value=1/rho_cu)
# FR4 conductance is approximately 0
mask_fr4 = (z > t_cu) & (z < (t_cu + t_fr4))
conductance_fr4 = 1/(1.e+12)
# conductance_fr4 = 0.
Conductances.setValue(conductance_fr4, where=mask_fr4)

# set the mask of where the current flows
top_layer_mask = FaceVariable(mesh=mesh, value=(z >= (t_cu + t_fr4)))
Conductances.setValue(0, where=top_layer_mask)

# Define largeValue to enforce the condition in masked areas
largeValue = 1.e7

# Define the desired current density vector
current = 70.
J_x = current / S
gradient_x = -J_x * rho_cu
gradient_vector = (gradient_x, 0., 0.) # The current is only in the x direction

# Extract faceNormals values to be used in cell equations, by averaging over the cells
faceNormals_x = mesh.faceNormals[0]  # Accessing the x-component of face normals
faceNormals_y = mesh.faceNormals[1]  # Accessing the y-component of face normals
faceNormals_z = mesh.faceNormals[2]  # Accessing the z-component of face normals

# Create the equation with the fixed current gradient
eq = (DiffusionTerm(coeff=Conductances)  # Original diffusion term
      + DiffusionTerm(coeff=largeValue * top_layer_mask)  # Large value diffusion term in masked regions
      - ImplicitSourceTerm(coeff=(top_layer_mask * largeValue *
                                  (gradient_vector[0] * faceNormals_x +
                                   gradient_vector[1] * faceNormals_y +
                                   gradient_vector[2] * faceNormals_z)).divergence))

[guyer]

1. With the previous code, if set conductance_fr4 = 0., there will be the warning and the output will be 0 which is not correct.

Understood. This is probably a consequence of the Robin condition, which I think we agree is not what you want anyway.

2. I am trying to set the gradient inside the shape and build the new equation with the following code (The change made starts from line 127 of the previous code).
   This is inspired by the example of Internal Fixed Gradient. And this method is trying to utilize the equation
   ∇ ϕ = − ρ i →

Why are you setting the gradient (or current) inside the shape? Where does it come from? Current needs to:

• come in from a boundary
• be induced by an electric field
• come in (out) at a source (sink)

The problem with this code is that it shows a constant 25 V voltage across the entire top surface no matter how big the current I set (but overall it shows a high voltage on the top surface and low voltage on the bottom surface which seems current).

I suspect that the gradient is not correct especially the part that finding the corresponding faceNormals direction might be wrong. But I am having difficulty finding the definition of elements of faceNormals. Would really appreciate if you could point out where I could find it.

Accepting that deus ex machina you somehow have a current within the domain, multiplying by the faceNormals isn't what you want. The faceNormals describe the orientation of the faces between the cells. This is a characteristic of the discretized mesh, not of the physical world. Your current presumably has an orientation in the physical world. By multiplying the constant current density by the faceNormals and then taking the divergence, you get zero contribution to the equation because the faceNormals describe a circuit around each cell, so the divergence of those orientations is zero.

[guyer]

I strongly recommend that, until you have things working satisfactorily, you work in 2D with, e.g.,

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

Both Gmsh and 3D just make troubleshooting and visualization unnecessarily harder.

[guyer]

multiplying by the faceNormals isn't what you want. The faceNormals describe the orientation of the faces between the cells. This is a characteristic of the discretized mesh, not of the physical world. Your current presumably has an orientation in the physical world.

I should clarify that the reason we use faceNormals for boundary conditions is because, physically, the flux normal to the boundary usually is what we want, e.g., $\hat{n}\cdot\vec{J} = J_0$, consistent with the natural boundary conditions from applying Gauss' theorom to the PDE.

But if there is another case in which I want this current to flow through the entire top layer with a non-uniform current distribution (different from the example above that the current has the same value and the same direction through the entire layer), is there a way I can set this current flow?

I would rearrange Ohm's law to give a Neumann expression for the gradient in potential in terms of the desired current density:
$\vec{i}=-\frac{1}{\rho}\vec{\phi} \rightarrow \vec{\phi} = -\rho\vec{i}$

On further thought, while you can do this to set a Neumann boundary condition, I think you're better off keeping the current and setting a fixed flux boundary condition.