# How to solve transient 3D heat equation with robin boundary conditions

Good afternoon!

I'm trying to solve the following heat equation:

with the following boundary conditions and initial value:

Nut I'm getting error while solving it with NDSolve:

s = NDSolve[{(1.0/alpha) D[T[x, y, z, t], t] ==
D[T[x, y, z, t], {x, 2}] + D[T[x, y, z, t], {y, 2}] +
D[T[x, y, z, t], {z, 2}],T[x, y, z, 0] == T0, -lambda D[T[x, y, 0, t], z] ==Piecewise[{{qmax Exp[-c ((x - xs0 - vx t)^2 + y^2)/r0^2], (x - xs0 - vx t)^2 + y^2 <=r0^2}, {-h (T[x, y, 0, t] - Tinf), (x - xs0 - vx t)^2 + y^2 >r0^2}}],lambda D[T[0, y, z, t], x] == h (T[0, y, z, t] - Tinf), -lambda D[T[L, y, z, t], x] == h (T[L, y, z, t] - Tinf), D[T[x, 0, z, t], y] == 0, -lambda D[T[x, W/2, z, t], y] == h (T[x, W/2, z, t] - Tinf), -lambda D[T[x, y, H, t], z] == h (T[x, y, H, t] - Tinf)}, T, {x, 0, L}, {y, 0, W}, {z, 0, H}, {t, 0, tf}]


Could someone help? Full code here: https://pastebin.com/U0qhNSeJ

• You need to define all parameters tf, alpha,lamda,... before NDSolve is applied. Oct 22 '19 at 14:20
• You might increase the chance to get a helpful answer if you provide complete executable code! Together with some information about the errors. Oct 22 '19 at 14:37
• Full code here: pastebin.com/U0qhNSeJ Oct 22 '19 at 14:48

We must change everything to meters. I updated the code and added a step function f[x], changed the boundary condition for y = 0 to Automatic (=NeumannValue[0, y==0]). Now all the pictures are in real time (sec) and in meters.

Needs["NDSolveFEM"]
L = 1000(*scale*);
(*Plate dimensions*)Ls = 250/L;
W = 200/L;
H = 30/L; reg = Cuboid[{0, 0, 0}, {Ls, W, H}]; mesh =
ToElementMesh[reg, MaxCellMeasure -> .0000005];
mesh["Wireframe"]
(*Material properties*)
rhocp = 4898556;
lambda = 36;
alpha = 1; ts = (lambda/rhocp)^(-1)(*time scale*);

T0 = 293(*Initial temperature*); Tinf = 293(*Ambient temperature*); Q \
= 4256(*Source power*); xs0 =
6/L(*Initial position of moving source*); xsf = (Ls -
xs0)(*final position of moving source*); r0 =
3/L(*radius of source*); c = 1(*source constant parameter*); vx =
ts 1.61/L(*moving source velocity x direction. Velocity *); qmax =
c Q/(Pi r0^2)/lambda(*source term*); h =
10/lambda(*heat convection coeffcient*); tf = (Ls - 2 xs0)/
vx(*Final time of calculation*);(*eq={(1.0/alpha) D[T[x,y,z,t],t]\
\[Equal]D[T[x,y,z,t],{x,2}]+D[T[x,y,z,t],{y,2}]+D[T[x,y,z,t],{z,2}],,-\
lambda D[T[x,y,0,t],z]\[Equal]Piecewise[{{qmax Exp[-c ((x-xs0-vx \
t)^2+y^2)/r0^2],(x-xs0-vx t)^2+y^2\[LessEqual]r0^2},{-h \
(T[x,y,0,t]-Tinf),(x-xs0-vx t)^2+y^2>r0^2}}],lambda D[T[0,y,z,t],x]\
\[Equal]h (T[0,y,z,t]-Tinf),-lambda D[T[L,y,z,t],x]\[Equal]h \
(T[L,y,z,t]-Tinf),D[T[x,0,z,t],y]\[Equal]0,-lambda D[T[x,W/2,z,t],y]\
\[Equal]h (T[x,W/2,z,t]-Tinf),-lambda D[T[x,y,H,t],z]\[Equal]h \
(T[x,y,H,t]-Tinf)}*)
f[x_] := (1 + Tanh[10000 x])/2
eq = D[T[x, y, z, t],
t] - (D[T[x, y, z, t], {x, 2}] + D[T[x, y, z, t], {y, 2}] +
D[T[x, y, z, t], {z, 2}]);

ic = T[x, y, z, 0] == T0; bc =
NeumannValue[
qmax Exp[-c ((x - xs0 - vx t)^2 + y^2)/r0^2] f[
r0^2 - (x - xs0 - vx t)^2 - y^2] -
h (T[x, y, z, t] - Tinf) f[((x - xs0 - vx t)^2 + y^2) - r0^2],
z == 0] +
NeumannValue[-h (T[x, y, z, t] - Tinf),
x == 0 || x == Ls || y == W || z == H];

sol = NDSolve[{eq == bc, ic},
T, {t, 0, tf}, {x, y, z} \[Element] mesh]
Table[DensityPlot[
Evaluate[T[x, y, 0, t] /. sol], {x, 0, Ls}, {y, 0, W},
PlotLegends -> Automatic, ColorFunction -> "Rainbow",
FrameLabel -> Automatic, PlotLabel -> Row[{"t =", t ts}],
PlotRange -> All], {t, .2 tf, tf, .2 tf}]

Table[DensityPlot[
Evaluate[T[x, 0, z, t] /. sol], {x, 0, Ls}, {z, 0, H},
PlotLegends -> Automatic, PlotRange -> All,
ColorFunction -> "Rainbow", FrameLabel -> Automatic], {t, .2 tf,
tf, .2 tf}]


We have nice pictures at z=H

Table[DensityPlot[
Evaluate[T[x, y, H, t] /. sol], {x, 0, Ls}, {y, 0, W},
PlotLegends -> Automatic, ColorFunction -> "Rainbow",
FrameLabel -> Automatic, PlotLabel -> Row[{"t =", t ts}],
PlotRange -> All], {t, .2 tf, tf, .2 tf}]


• Dear Alex, thank you very much for the solution! I'm new with Mathematica and I think I haven't understood all the steps you did. For example, on eq there is only one minus sign on the x derivative, should it also be negative for the y and z derivatives? I don't understand the boundary conditions... does one NeumannValue[h (T[x, y, z, t] - Tinf), True] accounts the convection on all surfaces? For the heat source term I couldn't follow how you did. Could you please help me? Thank you! Oct 23 '19 at 10:44
• 1. eq has a form D[T[x, y, z, t]-(…), and bracket () has the same meaning as in standard calculus. 2. I usually use True in such problems. But you can replace it with bc=NeumannValue[…,z==0]+NeumannValue[-h (T[x, y, z, t] - Tinf), x==0||y==0||x==Ls||y==W||z==H] and compare. 3. There is a thin source of radius r0 and a mesh with a mesh size of r0. Effectively, the source heats one cell. On the next cell we have Exp[-8]=0.000335463, there the cooling is turned on as (1-Exp[-8]). Oct 23 '19 at 11:48
• Is it possible to use the Piecewise[...] function to define the moving source as I did or NDSolve will raise error with it? Another problem is the symmetry bc at y=0, since the moving source is located at the center, i.e. at position W/2, for this case is it possible to also set .NeumannValue[0, y==0]? This is a Welding simulation and I'm also trying to solve it with finite difference method using ADI Douglas-Gunn method + TDMA. Oct 23 '19 at 14:49
• @JoãoVitor If this is welding, then explain what material it is and in what units the parameters are expressed? Oct 23 '19 at 16:23
• L, W, H, xs0, xsf and r0 in mm. tf in seconds. vx in mm/s. T0 and Tinf in °C. qmax in W. Q in J. c is dimensionless. lambda in W/m.K h in W/m².K rhocp in J/m³.K Oct 23 '19 at 17:16