Wrong solution from multi materials FEM NDSolve
With a small modification of code we have
Needs["NDSolve`FEM`"]
g = {0.250, 0.114, 0.040};(*thickness*)gw = Total[g];
λ = {8, 1.8, 44};
ρ = {3100, 2100, 7800};
cp = {1050, 1100, 540};
dc = Table[λ[[i]]/(ρ[[i]]*cp[[i]])/10^-5, {i,
Length[cp]}];
a[x_] := Piecewise[{{dc[[1]], 0 <= x < g[[1]]}, {dc[[2]],
g[[1]] <= x < g[[2]] + g[[1]]}, {dc[[3]], True}}]
σ =
QuantityMagnitude[
UnitConvert[Quantity["StefanBoltzmannConstant"]]] // N;
Trob = 1700.;
Tamb = 297;
h = 10;
ε = 0.85;
bc1 = DirichletCondition[
T[t, x] == Exp[-1000 t] + Trob/Tamb (1 - Exp[-1000 t]), x == 0.];
bc2 = 10^5/(ρ[[3]] cp[[3]]) NeumannValue[
h*(1 - T[t, x]) + ε*σ*Tamb^3 (1 - T[t, x]^4),
x == gw];
bc2rad = NeumannValue[ε*σ*Tamb^3 (1 - T[t, x]^4),
x == gw];
ic1 = T[0, x] == 1;
pde = D[T[t, x], t] - a[x]*D[T[t, x], x, x];
mesh = ToElementMesh[Line[{{0.}, {gw}}], MaxCellMeasure -> gw/404,
PrecisionGoal -> 5, AccuracyGoal -> 5]
sol = NDSolveValue[{pde == bc2, bc1, ic1}, T, {t, 0, .36},
x ∈ mesh, Method -> {"FiniteElement"}]
(*Visualization *)
{Plot[a[x]/10^5, {x, 0, gw}, PlotRange -> All, Frame -> True,
AxesOrigin -> {0, 0}, Filling -> Axis],
Plot3D[Tamb sol[10^-5 t, x], {t, 0, 36000}, {x, 0., gw},
AxesLabel -> Automatic, ColorFunction -> "Rainbow", Mesh -> None],
Plot[Table[Tamb sol[10^-5 t, x], {t, 2000, 36000, 2000}], {x, 0.,
gw}, ColorFunction -> "Rainbow"]}
I have not checked @Alex Trounev's answer, but this answer shows that there is good agreement between Mathematica and COMSOL Multiphysics.
Since you have a variety of thicknesses, I create a little routine so that I could mesh each region with the same number of elements (100 each).
Needs["NDSolve`FEM`"]
(* User Supplied Parameters *)
g = {0.25, 0.114, 0.04};(*thickness*)
gw = {0}~Join~Accumulate[g];
λ = {8, 1.8, 44};
ρ = {3100, 2100, 7800};
cp = {1050, 1100, 540};
(* Create a Multiregion Mesh *)
ClearAll[seg, appendCrdRight]
seg[thick_, nelm_, marker_] := Module[{crd, inc, marks},
crd = Subdivide[0, thick, nelm];
inc = Partition[Range[crd // Length], 2, 1];
marks = ConstantArray[marker, inc // Length];
<|"c" -> crd, "i" -> inc, "m" -> marks|>
]
appendCrdRight[a1_, a2_] := Module[{crd, inc, marks, len, lcrd},
len = a1["c"] // Length;
lcrd = a1["c"] // Last;
inc = Join[a1["i"], a2["i"] + len - 1];
crd = Join[a1["c"], Rest[a2["c"] + lcrd]];
marks = Join[a1["m"], a2["m"]];
<|"c" -> crd, "i" -> inc, "m" -> marks|>]
a = Fold[appendCrdRight, MapIndexed[seg[#1, 100, First[#2]] &, g]];
mesh = ToElementMesh["Coordinates" -> Partition[a["c"], 1],
"MeshElements" -> {LineElement[a["i"], a["m"]]},
"BoundaryElements" -> {PointElement[{{1}, {a["c"] // Length}}, {1,
2}]}];
Show[mesh["Wireframe"["MeshElementStyle" -> {Red, Green, Blue}]],
PlotRange -> {-0.01, 0.01}]
Now, we can set up the PDE system and solve it on our newly created mesh.
σ = First[UnitConvert[Quantity["StefanBoltzmannConstant"]]];
Trob = 1700;
Tamb = 297;
h = 10;
ε = 0.85;
rhocp = Evaluate[
Piecewise[{{ρ[[1]] cp[[1]], gw[[1]] <= x <= gw[[2]]},
{ρ[[2]] cp[[2]], gw[[2]] <= x <= gw[[3]]},
{ρ[[3]] cp[[3]], gw[[3]] <= x <= gw[[4]]}}]];
k = Evaluate[Piecewise[{{λ[[1]], gw[[1]] <= x <= gw[[2]]},
{λ[[2]], gw[[2]] <= x <= gw[[3]]},
{λ[[3]], gw[[3]] <= x <= gw[[4]]}}]];
bc1 = DirichletCondition[T[t, x] == Trob, x == 0];
bc2conv = NeumannValue[h*(Tamb - T[t, x]), x == Last@gw];
bc2rad = NeumannValue[ε*σ*(Tamb^4 - T[t, x]^4),
x == Last@gw];
ic1 = T[0, x] == Tamb;
op = Inactive[Div][{{-k}}.Inactive[Grad][T[t, x], {x}], {x}] +
rhocp*Derivative[1, 0][T][t, x];
pde = op == bc2conv + bc2rad;
sol = NDSolveValue[{pde, bc1, ic1},
T, {t, 0, 36000}, {x} ∈ mesh, StartingStepSize -> 0.01];
The model I set up in COMSOL Multiphysics (v 5.5) shows similar results to those shown in the OP.
For comparison purposes, I extracted the temperature data at each phase boundary point in COMSOL.
I exported these data to compare versus the Mathematica solution.
data = {{0, 1000, 2000, 3000, 4000, 5000, 6000, 7000, 8000, 9000,
10000, 11000, 12000, 13000, 14000, 15000, 16000, 17000, 18000,
19000, 20000, 21000, 22000, 23000, 24000, 25000, 26000, 27000,
28000, 29000, 30000, 31000, 32000, 33000, 34000, 35000,
36000}, {1700, 1700, 1700, 1700, 1700, 1700, 1700, 1700, 1700,
1700, 1700, 1700, 1700, 1700, 1700, 1700, 1700, 1700, 1700, 1700,
1700, 1700, 1700, 1700, 1700, 1700, 1700, 1700, 1700, 1700, 1700,
1700, 1700, 1700, 1700, 1700, 1700}, {297, 297.9169787`,
320.0562147`, 374.4552427`, 444.9013611`, 517.6131837`,
587.4876631`, 652.6604327`, 712.3644688`, 766.9603206`,
816.5391802`, 861.866491`, 902.8730203`, 940.4564489`,
974.5556695`, 1005.867455`, 1034.417079`, 1060.665637`,
1084.866141`, 1107.411419`, 1128.099762`, 1146.931167`,
1164.637928`, 1180.832645`, 1195.499525`, 1208.917884`,
1221.536363`, 1233.003818`, 1243.320249`, 1252.972747`,
1261.872597`, 1269.909554`, 1277.155111`, 1284.007597`,
1290.216067`, 1295.780522`, 1300.901468`}, {297, 297.0000101`,
297.0108185`, 297.2403045`, 298.3422144`, 301.3296677`,
306.8304462`, 315.0786727`, 326.0187665`, 339.0198185`,
353.9950315`, 370.1369655`, 387.5159699`, 405.1722292`,
423.1836315`, 440.8382141`, 458.14222`, 474.6735528`,
490.3439464`, 504.9171794`, 518.5145476`, 531.1360512`,
542.7808248`, 553.4493263`, 563.1415743`, 571.9455027`,
580.0023514`, 587.2015743`, 593.5431713`, 599.3724133`,
604.6264161`, 609.2270331`, 613.2390417`, 617.0233547`,
620.3526001`, 623.2267777`, 625.8287217`}, {297, 297.0000065`,
297.0084849`, 297.2058139`, 298.1991325`, 300.9831864`,
306.2034638`, 314.1201414`, 324.7019404`, 337.3400768`,
351.9481631`, 367.722907`, 384.7337123`, 402.0228897`,
419.6676093`, 436.9560503`, 453.8952359`, 470.0643493`,
485.3780489`, 499.6031165`, 512.8593059`, 525.1466173`,
536.4765686`, 546.8430665`, 556.2458626`, 564.7760878`,
572.5801167`, 579.5433842`, 585.6658902`, 591.2927421`,
596.3610853`, 600.7928104`, 604.6517643`, 608.293677`,
611.4944415`, 614.2540579`, 616.7511966`}};
Show[Plot[Evaluate[sol[t, #] & /@ gw], {t, 0, 36000}],
ListPlot[data[[2 ;; -1]], DataRange -> {0, 36000}]]
As you can see, there is very little difference between COMSOL (dots) and Mathematica (solid lines).
Update to Include the Basic Form
@AlexTrounev requested a comparison of the basic form to COMSOL as defined by:
$$\rho {{\hat C}_p}\frac{{\partial T}}{{\partial t}} - \lambda \frac{{{\partial ^2}T}}{{\partial {x^2}}} = 0$$
To use the FEM method, I recommend to cast your equations into coefficient form as shown FEM Tutorial.
$$\frac{{{\partial ^2}}}{{\partial {t^2}}}u + d\frac{\partial }{{\partial t}}u + \nabla \cdot\left( { - c\nabla u - \alpha u + \gamma } \right) + \beta \cdot\nabla u + au - f = 0$$
I find it easier to make comparisons of commercial solver (such as COMSOL) results to Mathematica results.
As shown with the following workflow, the Alex's basic form also matches COMSOL quite closely. I also included a case where I tried to thermal diffusivity in coefficient form and it fails to match COMSOL. Finally, it may be interesting to note that COMSOL's Laplace Equation Interface does not contain a Laplacian, rather:
$$\nabla \cdot \left( { - \nabla u} \right) = 0$$
(* User Supplied Parameters *)
g = {0.25, 0.114, 0.04};(*thickness*)
gw = {0}~Join~Accumulate[g];
λ = {8, 1.8, 44};
ρ = {3100, 2100, 7800};
cp = {1050, 1100, 540};
σ = First[UnitConvert[Quantity["StefanBoltzmannConstant"]]];
Trob = 1700;
Tamb = 297;
h = 10;
ε = 0.85;
bmesh = ToBoundaryMesh["Coordinates" -> Partition[gw, 1],
"BoundaryElements" -> {PointElement[{{1}, {2}, {3}, {4}}]}]; nrEle \
= 100; pt = Partition[gw, 2, 1]; mesh =
ToElementMesh[bmesh,
"RegionMarker" ->
Transpose[{Partition[(Mean /@ pt), 1], {1, 2, 3},
Abs[Subtract @@@ pt]/nrEle}]]
Show[mesh["Wireframe"["MeshElementStyle" -> {Red, Green, Blue}]],
PlotRange -> {-0.01, 0.01}]
rhocp = Evaluate[
Piecewise[{{ρ[[1]] cp[[1]], gw[[1]] <= x <= gw[[2]]},
{ρ[[2]] cp[[2]], gw[[2]] <= x <= gw[[3]]},
{ρ[[3]] cp[[3]], gw[[3]] <= x <= gw[[4]]}}]];
k = Evaluate[Piecewise[{{λ[[1]], gw[[1]] <= x <= gw[[2]]},
{λ[[2]], gw[[2]] <= x <= gw[[3]]},
{λ[[3]], gw[[3]] <= x <= gw[[4]]}}]];
bc1 = DirichletCondition[T[t, x] == Trob, x == 0];
bc2conv = NeumannValue[h*(Tamb - T[t, x]), x == Last@gw];
bc2rad = NeumannValue[ε*σ*(Tamb^4 - T[t, x]^4),
x == Last@gw];
ic1 = T[0, x] == Tamb;
(* Coefficient Form *)
op = Inactive[Div][{{-k}}.Inactive[Grad][T[t, x], {x}], {x}] +
rhocp*Derivative[1, 0][T][t, x];
pde = op == bc2conv + bc2rad;
Tcoef = NDSolveValue[{pde, bc1, ic1},
T, {t, 0, 36000}, {x} ∈ mesh, StartingStepSize -> 0.01];
(* Alex's "Basic Form" *)
op = rhocp*D[T[t, x], t] - k D[T[t, x], x, x];
pde = op == bc2conv + bc2rad;
Tbasic = NDSolveValue[{pde, bc1, ic1},
T, {t, 0, 36000}, {x} ∈ mesh, StartingStepSize -> 0.01];
(* Coefficient form with thermal diffusivity *)
bc2conv = NeumannValue[h*(Tamb - T[t, x])/rhocp, x == Last@gw];
bc2rad = NeumannValue[ε*σ*(Tamb^4 - T[t, x]^4)/
rhocp, x == Last@gw];
op = Inactive[Div][{{-k/rhocp}}.Inactive[Grad][T[t, x], {x}], {x}] +
Derivative[1, 0][T][t, x];
pde = op == bc2conv + bc2rad;
Talphainside =
NDSolveValue[{pde, bc1, ic1}, T, {t, 0, 36000}, {x} ∈ mesh,
StartingStepSize -> 0.01];
(* Plot Alex's "Basic Form" *)
Show[Plot[Evaluate[Tbasic[t, #] & /@ gw], {t, 0, 36000}],
ListPlot[data[[2 ;; -1]], DataRange -> {0, 36000}]]
(* Comparison of Methods *)
Show[Plot[Evaluate[Tcoef[t, #] & /@ gw], {t, 0, 36000},
PlotStyle -> ConstantArray[{Opacity[0.2], Thickness[0.015]}, 4]],
Plot[Evaluate[Talphainside[t, #] & /@ gw], {t, 0, 36000},
PlotStyle -> Dashed],
Plot[Evaluate[Tbasic[t, #] & /@ gw], {t, 0, 36000},
PlotStyle -> DotDashed]]
Once again, thank you to everyone who decided to help me in this calculation. As I wrote I have Mathematica since February 2020. I'm learning, but sometimes it's better to ask professionals.
Below is a solution that is based on MMA tutorials. Especially:
https://reference.wolfram.com/language/PDEModels/tutorial/HeatTransfer/HeatTransfer.html https://reference.wolfram.com/language/PDEModels/tutorial/HeatTransfer/ModelCollection/ShrinkFitting.html
I also used the elegant way of creating a 1D mesh given by @user21.
It should work.
Clear["Global`*"]
Needs["NDSolve`FEM`"]
HeatTransferModel[T_, X_List, k_, ρ_, Cp_, Velocity_, Source_] :=
Module[{V, Q, a = k},
V = If[Velocity === "NoFlow",
0, ρ*Cp*Velocity.Inactive[Grad][T, X]];
Q = If[Source === "NoSource", 0, Source];
If[FreeQ[a, _?VectorQ], a = a*IdentityMatrix[Length[X]]];
If[VectorQ[a], a = DiagonalMatrix[a]];
(*Note the-sign in the operator*)
a = PiecewiseExpand[Piecewise[{{-a, True}}]];
Inactive[Div][a.Inactive[Grad][T, X], X] + V - Q]
TimeHeatTransferModel[T_, TimeVar_, X_List, k_, ρ_, Cp_,
Velocity_, Source_] := ρ*Cp*D[T, {TimeVar, 1}] +
HeatTransferModel[T, X, k, ρ, Cp, Velocity, Source]
g = {0.25, 0.114, 0.04};
gw = {0}~Join~Accumulate[g];
bmesh = ToBoundaryMesh["Coordinates" -> Partition[gw, 1],
"BoundaryElements" -> {PointElement[{{1}, {2}, {3}, {4}}]}]; nrEle \
= 10; pt = Partition[gw, 2, 1]; mesh =
ToElementMesh[bmesh,
"RegionMarker" ->
Transpose[{Partition[(Mean /@ pt), 1], {1, 2, 3},
Abs[Subtract @@@ pt]/nrEle}]];
ρ1 = 3100;
Cp1 = 1050;
k1 = 8;
ρ2 = 2100;
Cp2 = 1100;
k2 = 1.8;
ρ3 = 7800;
Cp3 = 540;
k3 = 44;
parameters = {ρ ->
Piecewise[{{ρ1, ElementMarker == 1}, {ρ2,
ElementMarker == 2}, {ρ3, ElementMarker == 3}}],
Cp -> Piecewise[{{Cp1, ElementMarker == 1}, {Cp2,
ElementMarker == 2}, {Cp3, ElementMarker == 3}}],
k -> Piecewise[{{k1, ElementMarker == 1}, {k2,
ElementMarker == 2}, {k3, ElementMarker == 3}}]};
σ = First[UnitConvert[Quantity["StefanBoltzmannConstant"]]];
Tamb = 297;
h = 10;
Trob = 1700;
bc2conv = NeumannValue[h*(Tamb - T[t, x]), x == 0.404];
bc2rad = NeumannValue[0.85*σ*(297^4 - T[t, x]^4), x == 0.404];
ic1 = {T[0, x] == Tamb};
bc1 = DirichletCondition[T[t, x] == Trob, x == 0];
pde = {TimeHeatTransferModel[T[t, x], t, {x}, k, ρ, Cp, "NoFlow",
"NoSource"] == bc2conv + bc2rad, bc1, ic1} /. parameters;
sol = NDSolveValue[pde, T, {t, 0, 36000}, x ∈ mesh]
sol[36000, 0.404]
Plot[Table[sol[t, x], {t, 3600, 36000, 1800}], {x, 0, 0.404},
PlotRange -> {{0, 0.404}, {290, 1700}}, PlotTheme -> "Scientific",
ColorFunction -> "Rainbow"]
Likzew