Boundary condition for heat equation that allows "outflow" of temperature

5 months ago by
I just started using FEniCS (and FEM calculations too). I would like to solve a 2d heat equation problem where an external wind velocity $\overline{v}$v may transport the heat into a specific direction, e.g.
$\partial_tT_1=\frac{\alpha}{\rho_0c_p}I-\left(\overline{v}\cdot\nabla-\frac{\kappa}{\rho_0c_p}\nabla^2\right)T_1$tT1=αρ0cp I(v·κρ0cp 2)T1,
where  $T_1$T1  is the temperature change with respect to some ambient temperature,  $\alpha,\rho_0,\kappa$α,ρ0,κ are material constants, and  $I$I is an intensity distribution which I assume to be Gaussian for all times. My attempt to solve this problem resulted in the following code:
from fenics import *
import matplotlib.pyplot as plt
import numpy as np

# Define system parameter
rho0 = 1.176685    # mass density (kg.m^-3)
c_p = 1.0049e3     # specific heat capacity at constant pressure (^-1.K^-1)
kappa = 25.87e-3   # thermal conductivity (W.m^-1.K^-1)
T0 = 293.15        # ambient temperature
vel = [0.1, 0.05]  # wind speed (m.s^-1)
w0 = 10e-2         # beam radius
alpha = 1e-6       # absorption coefficient due to aerosols
P0 = 10e3          # beam power (W)

# Define spatial grid
w_extend = 10 # weight of the spatial extend
nx = ny = 100 # Number of grid points

# Define time steps
T = 10.0         # final time (s)
num_steps = 50   # number of time steps
dt = T/num_steps # time step size

# Derived parameter
K = kappa/rho0/c_p    # thermal diffusivity
vel = Constant((vel[0], vel[1]))   # convert to numpy array
I0 = 2*P0/np.pi/w0**2 # peak intensity
beta = alpha/rho0/c_p # prefactor in source term

# Create mesh and define function space
mesh = RectangleMesh(Point(-w_extend*w0/2, -w_extend*w0/2), Point(w_extend*w0/2, w_extend*w0/2), nx, ny, 'crossed')
V = FunctionSpace(mesh, 'P', 1)

# Define source term
I = Expression('I0*exp(-(pow(x[0], 2) + pow(x[1], 2))/pow(w0, 2))',
                 degree=2, I0=I0, w0=w0)
I_ini = interpolate(I, V)

# Define boundary condition
def boundary(x, on_boundary):
    return on_boundary

bc = DirichletBC(V, Constant(0), boundary)

# Define initial value
u_0 = Constant(0)         # ambient temperature with respect to T_0
u_n = interpolate(u_0, V) # save u_0 to u_n since this is kept during time stepping

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)

f = beta*I
F = u*v*dx + dt*dot(vel, grad(u))*v*dx  + dt*K*dot(grad(u), grad(v))*dx - (u_n + dt*f)*v*dx
a, L = lhs(F), rhs(F)

# Create VTK file for saving solution
vtkfile = File('solution/solution.pvd')

# Time-stepping
u = Function(V)
t = 0.0
vtkfile << (u, t)

for n in range(num_steps):
    # Update current time
    t += dt
    # Compute solution
    solve(a == L, u, bc)
    # Save to file and plot solution
    vtkfile << (u, t)
    print('t = {0:0.2f}'.format(t))
    # Update previous solution
The result looks OK instead for late times, when the wind results in a temperature change begins to rise at the edges of my geometry. I am sure this is due to my boundary condition which is  $T_1=0$T1=0 on the edges. However, this is only true for t=0 (where  $T_1=0$T1=0 everywhere) and I am searching for some smarter boundary condition. Is there something like an "open boundary" condition to avoid this problem? To illustrate my problem, I attached a figure: solution.pdf

Community: FEniCS Project

How about Newton cooling, i.e. a Robin boundary condition where the outward normal heat flux is proportional to the difference between domain and ambient temperature:
$\mathbf{q}\cdot\mathbf{n}=h\left(T-T_{\mathrm{amb}}\right)\quad\mathrm{on}\quad\partial\Omega$q·n=h(TTamb) on ∂Ω

written 5 months ago by klunkean  
Thank you for your reply. I am not sure how to apply this. h seems to be some constant? How big do I have to chose it? And I don't know how to implement this in FEniCS. From your equations I get the connections to the solution $T_1$T1 via $\overline{q}\cdot\overline{n}=-\kappa\nabla T_1\cdot\overline{n}$q·n=κT1·n. How can I implement this in FEniCS?
written 5 months ago by Phillip Springer  

Along the derivation of your weak form, when carrying out the integration by parts on the  $\nabla^2$2  term, you should encounter a boundary integral term which includes the outward normal heat flux (through reinserting Fourier's law). There you simply insert the cooling expression.

The constant is the heat transfer coefficient

written 5 months ago by klunkean  
Thanks again. I added the term 
- h/kappa*u*v*ds​

(where  $h=50\frac{W}{m^2K}$h=50Wm2K ) to F but this does not change much (still this pixel looking area where the temperature rises near the edges). I think this is because I still have the boundary condition that it should be 0 at the edges.

written 5 months ago by Phillip Springer  
Did you remove the DirichletBC? You can't prescribe essential *and* natural boundary conditions.
written 5 months ago by klunkean  
Ok, I think I got it. The term I added is not correct, a dt*K was missing. With that, an additional sign flip in that term and a much smaller h=1W/(m^2 K) it seems to work. Any idea why the sign flip is needed and why I cannot use an arbitrary h?
written 5 months ago by Phillip Springer  
I'll just consider the following simplified semi-discrete version:

$T=T_0+\Delta t\frac{\kappa}{\rho_0c_p}\nabla^2T$T=T0+Δtκρ0cp 2T

multiplying by a test function  $\delta T$δT  and integrating over the domain yields

$\int\left(T-T_0\right)\delta T\mathrm{d}x-\frac{\kappa\Delta t}{\rho_0c_p}\int\nabla^2T\delta T\mathrm{d}x=0$(TT0)δTdxκΔtρ0cp 2TδTdx=0
Integration by parts on the second integral gives

  $\int\left(T-T_0\right)\delta T\mathrm{d}x+\frac{\kappa\Delta t}{\rho_0c_p}\int\nabla T\cdot\nabla\delta T\mathrm{d}x-\frac{\kappa\Delta t}{\rho_0c_p}\int\nabla T\cdot n\delta T\mathrm{d}s=0$(TT0)δTdx+κΔtρ0cp T·δTdxκΔtρ0cp T·nδTds=0  

With Fourier's law  $q=-\kappa\nabla T$q=κT  we have

$\int\left(T-T_0\right)\delta T\mathrm{d}x+\frac{\kappa\Delta t}{\rho_0c_p}\int\nabla T\cdot\nabla\delta T\mathrm{d}x+\frac{\Delta t}{\rho_0c_p}\int q\cdot n\delta T\mathrm{d}s=0$(TT0)δTdx+κΔtρ0cp T·δTdx+Δtρ0cp q·nδTds=0

Now in the last surface integral there's the term where you put in your cooling law -- with a positive sign.
The choice of  $h$h  is everything but unambiguous and I don't have much experience with the choice. But roughly spoken, the larger the value, there more heat will be flowing outward.
written 5 months ago by klunkean  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »