Convective cooling of a bar using double newman boundary condition.

[sbtristan98]

Hi, I would like to model the convective cooling of a 1D 400K heated bar in air (300K).

This is my current implementation:

from fipy import CellVariable, Grid1D, DiffusionTerm, TransientTerm, GeneralSolver
import matplotlib.pyplot as plt
import matplotlib
matplotlib.use('TkAgg')

# 1. Define the domain and discretization
L = 2.0e-6  # Length of the bar in meters (2 µm)
nx = 100  # Number of cells/grid points
dx = L / nx  # Size of each cell

# 2. Initialize the mesh
mesh = Grid1D(nx=nx, dx=dx)

k = CellVariable(name="thermal conductivity", mesh=mesh, value=1.0, hasOld=True)

# 4. Create and initialize the temperature variable
Temperature = CellVariable(name="Temperature", mesh=mesh, value=400.0, hasOld=True)  # Initial value is 400K

TInfinite = 300.0
hTopSurface = 10.0
kAir = 0.02

eq = DiffusionTerm(coeff=k, var=Temperature) == 1000*TransientTerm(var=Temperature)

dt =1.00e-8

for sweep in range(10000):  # Sweep through simulation

    #Every 10 sweeps plot
    if sweep % 10 == 0:
        #Clear the plot
        plt.clf()
        plt.plot(mesh.cellCenters[0], Temperature.value)
        plt.legend(["Fipy", "Analytical"])
        plt.xlabel('Position (m)')
        plt.ylabel('Temperature (K)')
        plt.title('Temperature Distribution in the Bar')
        plt.grid(True)
        plt.show(block=False)
        plt.pause(0.01)

    Temperature.faceGrad.constrain(-hTopSurface * (TInfinite - Temperature.faceValue) / kAir, where=mesh.facesLeft)
    Temperature.faceGrad.constrain(-hTopSurface * (Temperature.faceValue - TInfinite) / kAir, where=mesh.facesRight)
    residual = eq.sweep(dt = dt, solver=GeneralSolver(iterations=10))

    if sweep % 20 == 0:
        dt = dt + 5.00e-8
    print("Residual: ", residual)

However, while it initially looks like it is functioning as expected (Cooling down from 400K to 360K is fine), as we approach 350K, the solver starts to diverge rapidly for some reason. I am already using sweeps to help convergence, but I can't get it to reach 300K stably... This final equilibrium with the air temperature is essential since I later intend on adding heat generation terms to my equation, but I can't seem to figure out why this simple model is non-stable...

Cheers,

Tristan

[guyer]

A few issues here:

• You have a Robin condition, not a gradient constraint.
• You need to .updateOld() in order to advance the time step. Sweeping just solves the same time step over and over again in an attempt to converge any nonlinearities.
  • Sweeping is something you do a small number of times at each time step. I didn't change this, but sweeping 10000 times is confusing.
• Constraints should only be applied once. FiPy uses lazy evaluation, so a single expression before you start solving will automatically update each time your solution changes.

With the following changes, your code seems to exhibit the behavior you're looking for:

diff --git a/discussion1075.py b/discussion1075.py
index 20ce589..21a8691 100644
--- a/discussion1075.py
+++ b/discussion1075.py
@@ -1,4 +1,4 @@
-from fipy import CellVariable, Grid1D, DiffusionTerm, TransientTerm, GeneralSolver
+from fipy import CellVariable, FaceVariable, Grid1D, DiffusionTerm, TransientTerm, GeneralSolver, ImplicitSourceTerm
 import matplotlib.pyplot as plt
 import matplotlib
 matplotlib.use('TkAgg')
@@ -11,7 +11,11 @@ dx = L / nx  # Size of each cell
 # 2. Initialize the mesh
 mesh = Grid1D(nx=nx, dx=dx)
 
-k = CellVariable(name="thermal conductivity", mesh=mesh, value=1.0, hasOld=True)
+mask = mesh.facesLeft | mesh.facesRight
+
+k0 = 1.0
+k = FaceVariable(name="thermal conductivity", mesh=mesh, value=k0)
+k.setValue(0., where=mask)
 
 # 4. Create and initialize the temperature variable
 Temperature = CellVariable(name="Temperature", mesh=mesh, value=400.0, hasOld=True)  # Initial value is 400K
@@ -20,7 +24,30 @@ TInfinite = 300.0
 hTopSurface = 10.0
 kAir = 0.02
 
-eq = DiffusionTerm(coeff=k, var=Temperature) == 1000*TransientTerm(var=Temperature)
+# 5. Define Robin boundary conditions
+
+# mesh.cellDistanceVectors isn't defined for Grid1D
+# https://stackoverflow.com/a/60080485/2019542
+# https://github.com/usnistgov/fipy/issues/706
+from fipy.tools import numerix
+MA = numerix.MA
+
+tmp = MA.repeat(mesh._faceCenters[..., numerix.NewAxis,:], 2, 1)
+cellToFaceDistanceVectors = tmp - numerix.take(mesh._cellCenters, mesh.faceCellIDs, axis=1)
+
+tmp = numerix.take(mesh._cellCenters, mesh.faceCellIDs, axis=1)
+tmp = tmp[..., 1,:] - tmp[..., 0,:]
+cellDistanceVectors = MA.filled(MA.where(MA.getmaskarray(tmp), cellToFaceDistanceVectors[:, 0], tmp))
+
+dPf = FaceVariable(mesh=mesh,
+                   value=mesh._faceToCellDistanceRatio * cellDistanceVectors)
+n = mesh.faceNormals
+a = FaceVariable(mesh=mesh, value=n, rank=1)
+b = FaceVariable(mesh=mesh, value=kAir / hTopSurface, rank=0)
+g = FaceVariable(mesh=mesh, value=TInfinite, rank=0)
+RobinCoeff = mask * k0 * n / (dPf.dot(a) + b)
+
+eq = DiffusionTerm(coeff=k, var=Temperature) + (RobinCoeff * g).divergence - ImplicitSourceTerm(coeff=(RobinCoeff * a.dot(n)).divergence, var=Temperature) == 1000*TransientTerm(var=Temperature)
 
 dt =1.00e-8
 
@@ -39,10 +66,9 @@ for sweep in range(10000):  # Sweep through simulation
         plt.show(block=False)
         plt.pause(0.01)
 
-    Temperature.faceGrad.constrain(-hTopSurface * (TInfinite - Temperature.faceValue) / kAir, where=mesh.facesLeft)
-    Temperature.faceGrad.constrain(-hTopSurface * (Temperature.faceValue - TInfinite) / kAir, where=mesh.facesRight)
     residual = eq.sweep(dt = dt, solver=GeneralSolver(iterations=10))
 
     if sweep % 20 == 0:
         dt = dt + 5.00e-8
+        Temperature.updateOld()
     print("Residual: ", residual)

I encourage you to work though examples.diffusion.mesh1D in detail before attempting to set up your own problems.

[sbtristan98]

Thank you very much! This solves the issue I had.