2
$\begingroup$

I want to solve the Poisson PDE for heat flow in a 3-D solid cube with given dimensions $x$, $y$, and $z$:

$$\rho C\frac{\partial T}{\partial t} = k \Delta T$$

The cube is irradiated with a constant heat flux $I(x,y)$ at the $z=0$ surface. The boundaries should meet the Neumann condition for radiative cooling given by $\frac{\partial T}{\partial x_i} = \pm\epsilon\sigma(T^4-300^4)$, where the sign depends on wether you look at $x_i=$max or $x_i=0$. Due to the irradiation the $z=0$ surface is an exception, where the condition should rather be $\frac{\partial T}{\partial x_i} = \epsilon\sigma(T^4-300^4) + I(x,y)$. The last term represents the intensity distribution of the radiative flux.

I decided to go for a numerical approach I found in this thread using Python. As a result, I came up with the following piece of code.

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits import mplotn3d

dt = 3
di = 0.01

#Thermal conductivity
k = 1.38
#Density
rho = 2202
#Specific Heat capacity
cp = 745
#Thermal diffusivity
alpha = k/(rho*cp)
#Emissivity
emissivity = 0.79

sigma = 5.670367e-8

y_max = 0.3
x_max = 0.3
z_max = 0.3
t_max = 60

P = 1
omega = 0.02

#Function that generates a gaussian intensity pattern I(x,y) with power P and beam waist omega

def get_beam(P, omega, di, x_max, ymax, z_max):
  x = np.arange(0,x_max+di,di)
  y = np.arange(0,y_max+di,di)
  z = np.arange(0,z_max+di,di)
  t = np.arange(0,t_max+dt,dt)
  r = len(t)
  cx = len(x)
  cy = len(y)
  cz = len(z)
  I = np.zeros([r, cx, cy, cz])
  for jx in range(0,cx-1):
      for jy in range(0,cy-1):
          I[:,jx,jy,0] = 2*P/(np.pi*omega**2)*np.exp(-2*(abs(di*(jx-cx/2))**2 + abs(di*(jy-cy/2))**2)/(omega**2))
  return I

def FTCS(dt,dy,t_max,x_max,y_max,z_max,k,T0, I):
  s = alpha*dt/di**2
  x = np.arange(0,x_max+di,di)
  y = np.arange(0,y_max+di,di)
  z = np.arange(0,z_max+di,di)
  t = np.arange(0,t_max+dt,dt)
  r = len(t)
  cx = len(x)
  cy = len(y)
  cz = len(z)
  #Initialize mesh
  T = np.ones([r,cx, cy, cz])*T0
  #time loop
  for n in range(0,r-1):  
      print(n)
      #y loop
      for jy in range(0,cy-1):   
          #x loop 
          for jx in range(0,cx-1):   
              #z loop 
              for jz in range(0,cz-1):  
                  #Compute T at next time step: Tn+1 = Tn + d^2T/dx^2 + d^2T/dy^2 + d^2T/dz^2  
                  T[n+1,jx,jy,jz] = T[n,jx,jy,jz] + s*(T[n,jx-1,jy,jz] - 2*T[n,jx,jy,jz] + T[n,jx+1,jy,jz]) + s*(T[n,jx,jy-1,jz] - 2*T[n,jx,jy,jz] + T[n,jx,jy+1,jz]) + s*(T[n,jx,jy,jz-1] - 2*T[n,jx,jy,jz] + T[n,jx,jy,jz+1]) 


      #Boundary Condition
      #Set boundary condition for x on both sides   
      T[n+1,cx-1,:,:] = T[n,cx-1,:,:] + s*(T[n,cx-2,:,:] - 2*T[n,cx-1,:,:] + T[n,cx-2,:,:] - 2*di*emissivity*sigma/k*(T[n,cx-1, :,:]**4-300**4))#- 2*di*4*emissivity*sigma/k*(T[n,cx-1,:,:]*300**3))
      T[n+1,0,:,:] = T[n,0,:,:] + s*(T[n,1,:,:] - 2*T[n,0,:,:] + T[n,1,:,:] - 2*di*emissivity*sigma/k*(T[n,0,:,:]**4-300**4))#- 2*di*4*emissivity*sigma/k*(T[n,0,:,:]*300**3))

      #Set boundary condition for y on both sides 
      T[n+1,:,cy-1,:] = T[n,:,cy-1,:] + s*(T[n,:,cy-2,:] - 2*T[n,:,cy-1,:] + T[n,:,cy-2,:] - 2*di*emissivity*sigma/k*(T[n,:,cy-1,:]**4-300**4))#- 2*di*4*emissivity*sigma/k*(T[n,:,cy-1,:]*300**3))
      T[n+1,:,0,:] = T[n,:,0,:] + s*(T[n,:,1,:] - 2*T[n,:,0,:] + T[n,:,1,:] - 2*di*emissivity*sigma/k*(T[n,:,0,:]**4-300**4))#- 2*di*4*emissivity*sigma/k*(T[n,:,0,:]*300**3))

      #Set boundary condition for z on both sides 
      T[n+1,:,:,cz-1] = T[n,:,:,cz-1] + s*(T[n,:,:,cz-2] - 2*T[n,:,:,cz-1] + T[n,:,:,cz-2] - 2*di*emissivity*sigma/k*(T[n,:,:,cz-1]**4-300**4))#- 2*di*4*emissivity*sigma/k*(T[n,:,:,cz-1]*300**3))
      T[n+1,:,:,0] = T[n,:,:,0] + s*(T[n,:,:,1] - 2*T[n,:,:,0] + T[n,:,:,1] - 2*di*emissivity*sigma/k*(T[n,:,:,0]**4-300**4) + 2*di/k*I[n,:,:,0])#- 2*di*4*emissivity*sigma/k*(T[n,:,:,0]*300**3) + 2*di*I[n,:,:,0])

  return x,y,z,t,T,r,s 

I = get_beam(P, omega, di, x_max, y_max, z_max)
x,y,z,t,T,r,s = FTCS(dt,di,t_max,x_max,y_max,z_max,k,300,I)


fig = plt.figure()
X, Y = np.meshgrid(x, y)
ax = plt.axes(projection='3d')
ax.plot_surface(X, Y, T[-1,:,:,0])
ax.set_xlabel('x-axis')
ax.set_ylabel('y-axis')

The result looks promising but is not equal to references I found. Furthermore, the result scales rapidly with the mesh size $d_i$ and crashes for huge point densities. I couldn't figure out whether and if so where I did a mistake. Can you help me?

//-----------------------------------//

I updated the code to include the correct boundaries. The loop looks like this now. I know its a bit clunky, but I thought I take care of it later. Nevertheless, the dependence from the mesh size is still there. I am afraid that I just forgot a $di$ somewhere, but I couldn't find out where.

    for n in range(0,r-1):
    print(n)
    for jy in range(0,cy-1):
        for jx in range(0,cx-1):
            for jz in range(0,cz-1):
                T[n+1,jx,jy,jz] = T[n,jx,jy,jz] + s*(T[n,jx-1,jy,jz] - 2*T[n,jx,jy,jz] + T[n,jx+1,jy,jz]) + s*(T[n,jx,jy-1,jz] - 2*T[n,jx,jy,jz] + T[n,jx,jy+1,jz]) + s*(T[n,jx,jy,jz-1] - 2*T[n,jx,jy,jz] + T[n,jx,jy,jz+1]) 

                if jx == cx-1:
                    T[n+1,jx,jy,jz] = T[n,jx,jy,jz] + s*(T[n,jx-1,jy,jz] - 2*T[n,jx,jy,jz] + T[n,jx-1,jy,jz] + 2*di*emissivity*sigma/k*(T[n,jx, jy, jz]**4-300**4)) + s*(T[n,jx,jy-1,jz] - 2*T[n,jx,jy,jz] + T[n,jx,jy+1,jz]) + s*(T[n,jx,jy,jz-1] - 2*T[n,jx,jy,jz] + T[n,jx,jy,jz+1])
                elif jx == 0:
                    T[n+1,jx,jy,jz] = T[n,jx,jy,jz] + s*(T[n,jx-1,jy,jz] - 2*T[n,jx,jy,jz] + T[n,jx-1,jy,jz] - 2*di*emissivity*sigma/k*(T[n,jx, jy, jz]**4-300**4)) + s*(T[n,jx,jy-1,jz] - 2*T[n,jx,jy,jz] + T[n,jx,jy+1,jz]) + s*(T[n,jx,jy,jz-1] - 2*T[n,jx,jy,jz] + T[n,jx,jy,jz+1]) 
                elif jy == cy-1:
                    T[n+1,jx,jy,jz] = T[n,jx,jy,jz] + s*(T[n,jx-1,jy,jz] - 2*T[n,jx,jy,jz] + T[n,jx+1,jy,jz]) + s*(T[n,jx,jy-1,jz] - 2*T[n,jx,jy,jz] + T[n,jx,jy-1,jz] + 2*di*emissivity*sigma/k*(T[n,jx,jy,jz]**4-300**4)) + s*(T[n,jx,jy,jz-1] - 2*T[n,jx,jy,jz] + T[n,jx,jy,jz+1]) 
                elif jy == 0:
                    T[n+1,jx,jy,jz] = T[n,jx,jy,jz] + s*(T[n,jx-1,jy,jz] - 2*T[n,jx,jy,jz] + T[n,jx+1,jy,jz]) + s*(T[n,jx,jy-1,jz] - 2*T[n,jx,jy,jz] + T[n,jx,jy-1,jz] - 2*di*emissivity*sigma/k*(T[n,jx,jy,jz]**4-300**4)) + s*(T[n,jx,jy,jz-1] - 2*T[n,jx,jy,jz] + T[n,jx,jy,jz+1]) 
                elif jz == cz-1:
                    T[n+1,jx,jy,jz] = T[n,jx,jy,jz] + s*(T[n,jx-1,jy,jz] - 2*T[n,jx,jy,jz] + T[n,jx+1,jy,jz]) + s*(T[n,jx,jy-1,jz] - 2*T[n,jx,jy,jz] + T[n,jx,jy+1,jz]) + s*(T[n,jx,jy,jz-1] - 2*T[n,jx,jy,jz] + T[n,jx,jy,jz-1] + 2*di*emissivity*sigma/k*(T[n,jx,jy,jz]**4-300**4)) 
                elif jz == 0:
                    T[n+1,jx,jy,jz] = T[n,jx,jy,jz] + s*(T[n,jx-1,jy,jz] - 2*T[n,jx,jy,jz] + T[n,jx+1,jy,jz]) + s*(T[n,jx,jy-1,jz] - 2*T[n,jx,jy,jz] + T[n,jx,jy+1,jz]) + s*(T[n,jx,jy,jz-1] - 2*T[n,jx,jy,jz] + T[n,jx,jy,jz-1] - 2*di*emissivity*sigma/k*(T[n,jx,jy,jz]**4-300**4) + 2*di*I[n,jx,jy,0])
Improve this question
$\endgroup$
13
  • $\begingroup$ At least one of your problems is that the forward-in-time (Euler) integration scheme is only conditionally stable. For stability, you must have $ \Delta t \le (\Delta x)^2/2 $. I suspect your large $\Delta t$ is why your solution "crashes" for fine meshes. For more information search for Courant, Friedrichs, Lewy condition. $\endgroup$
    – Bill Greene
    Commented Jul 25, 2018 at 19:57
  • $\begingroup$ Thank you for the comment. The crash-problems is pretty much solved when I adjust the time resolution. While writing this thread I found another mistake regarding the boundary conditions. I neglected the second derivative of the two remaining dimensions. So I included additional if statements in the inner loop: $\endgroup$
    – Phillip
    Commented Jul 25, 2018 at 20:03
  • $\begingroup$ I have updated my initial post with the new piece of code I added. Still, the other Problem occurs. $\endgroup$
    – Phillip
    Commented Jul 26, 2018 at 15:15
  • $\begingroup$ @BillGreene But the CFL condition is only for convection/advection and NOT for diffusion. For diffusion it should be $\Delta t \leq \frac{\left(\Delta x\right)^2 c_p \rho}{2 \lambda}$ (for one dimension). $\endgroup$
    – JE_Muc
    Commented Jul 26, 2018 at 16:15
  • $\begingroup$ @Phillip What is the ramining problem? The failing equality to the reference? $\endgroup$
    – JE_Muc
    Commented Jul 26, 2018 at 16:26

1 Answer 1

Reset to default
1
$\begingroup$

Thanks for asking your question. It really made me delve alot deeper into the topic of boundary conditions. I am still a rookie in this topic, but I'll try to answer your question anyways.

Your major problem seems to be that your units are not correct. Your equation for radiative heat flux has the unit $[\frac{\text{W}}{\text{m}^2}]$, while the Neumann boundary condition needs a unit of $[\frac{\text{K}}{\text{m}}]$. Your equation for the heat flux should say:

$$\frac{dq}{dt} = \epsilon \sigma \left(T^4 - 300^4 \right) + I(x,y)$$

You almost get your units correctly by multiplying this with $di/k$, resulting in $[\text{K}]$ for the radiative Stefan-Boltzmann heat transfer, but for the heat flux based on constant heat flux, you forgot the $1/k$ while still multiplying with $di$, resulting in $[\frac{\text{W}}{\text{m}}]$.

What I am not sure about is the factor of $2$. I guess this stems from a central difference scheme and thus shouldn't be correct here. But if anyone knows an explanation why this is correct/incorrect, I'd really like to know it.
Furthermore I am not sure if multiplying with $s$ is correct here. I guess you need to try all cases to and compare it with the reference solution, unless someone else with more knowledge in this topic answers.

Improve this answer
$\endgroup$
2
  • $\begingroup$ Thank you very much for the comprehensive answer. You are absolutely right, that die Neumann condition should have units of K/m. This should be corrected by the 'k' in '2*diemissivitysigma/k*(T[n,:,0,:]**4-300**4)), which has the units 'W/mK' and therefor cancels out the wrong units. But I forgot the 'k' right before 'I(x,y)'. $\endgroup$
    – Phillip
    Commented Jul 27, 2018 at 15:12
  • $\begingroup$ Did it help you finding your problem? $\endgroup$
    – JE_Muc
    Commented Jul 30, 2018 at 13:36

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.