0
$\begingroup$

Good afternoon!

I'm trying to solve the following heat equation:

enter image description here

with the following boundary conditions and initial value: enter image description here enter image description here enter image description here

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

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

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["NDSolve`FEM`"]
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}]

Figure 1

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}]

Figure 2

$\endgroup$
8
  • $\begingroup$ 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! $\endgroup$ Oct 23 '19 at 10:44
  • $\begingroup$ 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]). $\endgroup$ Oct 23 '19 at 11:48
  • $\begingroup$ 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. $\endgroup$ Oct 23 '19 at 14:49
  • $\begingroup$ @JoãoVitor If this is welding, then explain what material it is and in what units the parameters are expressed? $\endgroup$ Oct 23 '19 at 16:23
  • $\begingroup$ 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 $\endgroup$ Oct 23 '19 at 17:16

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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