Implementing Convective and Stephan Boltzmann Radiative cooling of a 1D Rod using robin boundary condition.

[sbtristan98]

Hello,

I am a bit confused as to how to create the robin boundary condition for radiative cooling.

I had a read through https://www.ctcms.nist.gov/fipy/documentation/USAGE.html and #1007

Let's say we have this equation: -k∇Tn = h(T-Tinf) + ϵσ[T⁴-Tinf⁴]
According to the documentation, we need to rearrange it into the form: n(aϕ + b∇ϕ) = g

Therefore, based on my understanding, this would resolve to:
a = hn, b = k, g = hTinf + ϵσ[Tinf⁴ - T⁴]

However, I am unsure what the best way is to add the information from T⁴ into the boundary equation.

I was thinking of implementing it like this:

from fipy import CellVariable, FaceVariable, Grid1D, DiffusionTerm, TransientTerm, GeneralSolver, ImplicitSourceTerm
from fipy.tools import numerix
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)

mask = mesh.facesLeft | mesh.facesRight

k0 = 100.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

TInfinite = 300.0
hTopSurface = 10.0
gold_emissivity = 0.3
Stephan_Boltzmann_constant = 5.67e-8 #W/(m^2*K^4)

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=hTopSurface*n, rank=1)
b = FaceVariable(mesh=mesh, value=k0, rank=0)

Stephan_Boltzmann_Radiation = gold_emissivity*Stephan_Boltzmann_constant*(numerix.power(TInfinite,4) - numerix.power(Temperature[-1],4))

g = FaceVariable(mesh=mesh, value=hTopSurface*TInfinite + Stephan_Boltzmann_Radiation, 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

for sweep in range(10000):  # Sweep through simulation
    plt.show(block=False)
    plt.pause(0.01)
    print(g)

    #Every 10 sweeps plot
    if sweep % 10 == 0:
        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)

    residual = eq.sweep(dt = dt, solver=GeneralSolver(iterations=10))

    if sweep % 20 == 0:
        dt = dt + 5.00e-8
        Temperature.updateOld()

However, the FaceVariable g never updates, even though the temperature is constantly changing due to convection...

Cheers,

Tristan

[guyer]

• You need a variable expression, not a variable with its initial condition set to the expression
  • g = FaceVariable(mesh=mesh, value=hTopSurface*TInfinite + Stephan_Boltzmann_Radiation, rank=0) declares a FaceVariable with the instantaneous value (at time of definition) of hTopSurface*TInfinite + Stephan_Boltzmann_Radiation.
  • g = (hTopSurface*TInfinite + Stephan_Boltzmann_Radiation) * mask declares a variable expression with the value hTopSurface*TInfinite + Stephan_Boltzmann_Radiation that updates whenever its constituents update.
• Don't assume you know how data is arranged in FiPy. Temperature[-1] is the last value in the Temperature array, but you should not assume that corresponds to the right-most value. FiPy is designed around arbitrary, irregular meshes.

diff --git a/discussion1076.py b/discussion1076.py
index c0b5879..ba04e61 100644
--- a/discussion1076.py
+++ b/discussion1076.py
@@ -41,9 +41,9 @@ n = mesh.faceNormals
 a = FaceVariable(mesh=mesh, value=hTopSurface*n, rank=1)
 b = FaceVariable(mesh=mesh, value=k0, rank=0)
 
-Stephan_Boltzmann_Radiation = gold_emissivity*Stephan_Boltzmann_constant*(numerix.power(TInfinite,4) - numerix.power(Temperature[-1],4))
+Stephan_Boltzmann_Radiation = gold_emissivity*Stephan_Boltzmann_constant*(TInfinite**4 - Temperature.faceValue**4)
 
-g = FaceVariable(mesh=mesh, value=hTopSurface*TInfinite + Stephan_Boltzmann_Radiation, rank=0)
+g = (hTopSurface*TInfinite + Stephan_Boltzmann_Radiation) * mask
 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)

Note: there's nothing intrinsically wrong with numerix.power(TInfinite,4); it's just twice as much typing for no added benefit.

[Sandip433]

I want to implement this equation to check but the calculated result is different from analytical result

the code is
`

from fipy import *
import matplotlib.pyplot as plt
from tqdm import tqdm
L = 10
nx = 1000
dx = L / nx
mesh = Grid1D(nx=nx, dx=dx)
x = mesh.cellCenters[0]
u = CellVariable(mesh=mesh, value=100 * numerix.sin (numerix.pi * x/10), hasOld=True)
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
mask = mesh.facesLeft | mesh.facesRight
k0 = 0.01
k = FaceVariable(mesh=mesh, value=k0)
k.setValue(0., where=mask)
a1 = FaceVariable(mesh=mesh, value=2 * n, rank=1)
b1 = FaceVariable(mesh=mesh, value=3, rank=0)
g1 = 5
a2 = FaceVariable(mesh=mesh, value=4 * n, rank=1)
b2 = FaceVariable(mesh=mesh, value=-1, rank=0)
g2 = 10
RobinCoeff1 = mesh.facesLeft * k0 * n / (dPf.dot(a1) + b1)
RobinCoeff2 = mesh.facesRight * k0 * n / (dPf.dot(a2) + b2)
A1 = (RobinCoeff1 * g1).divergence - ImplicitSourceTerm(coeff=(RobinCoeff1 * a1.dot(n)).divergence, var=u)
A2 = (RobinCoeff2 * g2).divergence - ImplicitSourceTerm(coeff=(RobinCoeff2 * a2.dot(n)).divergence, var=u)
eq = TransientTerm(var=u) == DiffusionTerm(coeff=k, var=u) + A1 + A2
dt = 1.00e-3
T = 10
time = Variable(0.)
solver = LinearPCGSolver(precon=None, tolerance=1e-10, iterations=2000)
progress_bar = tqdm(total=int(T/dt), desc='Solving PDE')
while time()<T:
u.updateOld()
for i in range(1):
res=eq.sweep(dt=dt, solver=solver)
time.setValue(time() + dt)
progress_bar.update()
progress_bar.close()
plt.plot(x, u.value, 'red')
CAnalytical = CellVariable(mesh=mesh, value= 100 * numerix.sin (numerix.pi * x/10) * numerix.exp(-T * (numerix.pi * x/100) **2))
plt.plot(x, CAnalytical, 'blue')
plt.xlabel('x')
plt.ylabel('u')
plt.show`