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$
  • $\begingroup$ You need to define all parameters tf, alpha,lamda,... before NDSolve is applied. $\endgroup$ – Ulrich Neumann 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$ – Ulrich Neumann Oct 22 '19 at 14:37
  • $\begingroup$ Full code here: pastebin.com/U0qhNSeJ $\endgroup$ – João Vitor 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

|improve this answer|||||
$\endgroup$
  • $\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$ – João Vitor 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$ – Alex Trounev 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$ – João Vitor 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$ – Alex Trounev 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$ – João Vitor 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.