How Can Mathematica Be Used to Plot Transient Heating in a Slab?

[Dustinsfl]

You can find problems with downloadable notebooks now http://www.mathhelpboards.com/f49/engineering-analysis-notes-2882/.

If one boundary is insulated and the other is subjected to and held at a temperature of unity, we wish to determine the solution for the transient heating of the slab.
The governing equation is the usual 1-D heat equation and the boundary conditions (mixed) are given by
\begin{alignat*}{3}
T_x(0,t) & = & 0\\
T(L,t) & = & 1
\end{alignat*}
with initial conditions
$$
T(x,0) = 0.
$$
Obtain the solution to this problem.
For the special case $\alpha = 1$ and $L = \pi$ plot a sequence of the temperature profiles between the initial state and the steady-state, construct a contour plot, density plot, and 3D plot.

Once we solve the problem, we obtain the solution as
$$
T(x,t) = 1 + \frac{4}{\pi}\sum_{n = 1}^{\infty}\frac{(-1)^n}{2n - 1}\cos\left[\left(n - \frac{1}{2}\right)x\right]\exp\left[-\left(n - \frac{1}{2}\right)^2t\right]
$$
Next, we construct all the plots using Mathematica

Nmax = 40;
L = Pi;
\[Lambda] = Table[(n - 1/2)*Pi/L, {n, 1, Nmax}];
\[Alpha] = 1;
MyTime = Table[t, {t, 0.0001, 1, .05}];
f[x_] = -1;

A = Table[2/L*Integrate[f[x]*Cos[\[Lambda][[n]]*x], {x, 0, L}], {n, 1,Nmax}];

u[x_, t_] = 1+Sum[A[[n]]*Cos[\[Lambda][[n]]*x]*E^{-\[Alpha]*\[Lambda][[n]]^2*t}, {n, 1, Nmax}];

Plot[u[x, MyTime], {x, 0, L}, PlotStyle -> {Red}]

This will produce the graph
https://www.physicsforums.com/attachments/390._xfImport

Adding in this piece of code will produce an animation between the different time profiles.

Animate[Plot[u[x, t], {x, 0, L}, PlotRange -> {0, 1.1}, GridLines -> Automatic, Frame -> True, PlotStyle -> {Thick, Red}], {t,0, 1, 0.02}, 
 AnimationRunning -> False]

We can produce the contour plot by adding in

ContourPlot[u[x, y], {x, 0, L}, {y, 0, L}, PlotRange -> All, ColorFunction -> "Rainbow"]

The end result is

https://www.physicsforums.com/attachments/393._xfImport

The density plot

DensityPlot[u[x, y], {x, 0, L}, {y, 0, L}, PlotRange -> All, ColorFunction -> "Rainbow"]

The plot is
https://www.physicsforums.com/attachments/392._xfImport

And lastly the 3d plot

Plot3D[u[x, y], {x, 0, L}, {y, 0, L}, PlotRange -> All, Boxed -> False, ColorFunction -> "Rainbow"]

https://www.physicsforums.com/attachments/391._xfImport
Here we can see the Gibbs Phenomenon occurring.

 

[Dustinsfl]

Consider the transient diffusion of heat in a 1-D slab initially with a temperature $f(x)$ that is then subjected to symmetric convective cooling of equal strength at its two boundaries $x = -L$ and $x = L$.
Recognizing that the problem will have a symmetry plane about $x = 0$, the problem may be then treated more simply as one with an insulated condition at $x = 0$ and a convection condition at $x = L$.
With proper non-dimensionalization and assuming for convenience that the ratio of convection and conduction parameters $h/k = 1$, the cooling process is described by the following equations:
\begin{alignat*}{3}
T_t & = & T_{xx}\\
T_x(0,t) & = & 0\\
T_x(1,t) + T(1,t) & = & 0\\
T(x,0) & = & f(x)
\end{alignat*}
For the particular case of $f(x) = 1$, numerically determine the series coefficients $A_n$ and construct a series representation for $T(x,t)$.
Using this series approximation, plot the temperature profiles on the same set of axes for a range of time steps in the dimensionless time interval $t\in [0,4]$ in order to visualize the process.
Robin boundary conditions.
With this problem, let's identify the first 40 eigenvalues and construct the same plots.
In solving this problem, we have that the eigenvalues are
$$
\tan\lambda_n = \frac{1}{\lambda_n}.
$$

Nmax = 40;
MyTime = Table[t, {t, 0.0001, 4, .2}];
f[x_] = 1;
Plot[{Tan[x] , 1/x}, {x, 0, 5 Pi}]

Here we can look at the plot of the solution to $\tan\lambda_n = \frac{1}{\lambda_n}$.
View attachment 367
Now, let's find the first 40 eigenvalues numerically.

\[Lambda] = Flatten[Table[FindRoot[Tan[x] - 1/x == 0, {x, (n - 1/2)*\[Pi] - 0.0001}], {n, Nmax}]][[All, 2]]

{0.860334, 3.42562, 6.4373, 9.52933, 12.6453, 15.7713, 18.9024, 
22.0365, 25.1724, 28.3096, 31.4477, 34.5864, 37.7256, 40.8652, 
44.005, 47.1451, 50.2854, 53.4258, 56.5663, 59.707, 62.8478, 65.9886, 
69.1295, 72.2705, 75.4115, 78.5525, 81.6936, 84.8348, 87.976, 
91.1172, 94.2584, 97.3996, 100.541, 103.682, 106.824, 109.965, 
113.106, 116.248, 119.389, 122.53}

Let's setup the rest of our problem now so we can produces the plots and obtain the coefficients of the first 40 terms.

A = Table[2*Integrate[f[x]*Cos[\[Lambda][[n]]*x], {x, 0, 1}], {n, 1, Nmax}]

{1.76225, -0.163604, 0.0476919, -0.0219042, 0.0124686, -0.00802462, 
0.0055897, -0.00411432, 0.00315382, -0.00249397, 0.00202131, 
-0.00167123, 0.00140477, -0.00119727, 0.00103256, -0.00089962, 
0.00079079, -0.000700571, 0.000624951, -0.000560943, 0.000506285, 
-0.000459243, 0.000418464, -0.000382884, 0.000351655, -0.000324096, 
0.000299655, -0.000277877, 0.000258389, -0.000240882, 0.000225095, 
-0.000210811, 0.000197844, -0.000186038, 0.000175258, -0.000165388, 
0.000156329, -0.000147995, 0.000140309, -0.000133208}

u[x_, t_] = Sum[A[[n]]*Cos[\[Lambda][[n]]*x]*E^{-\[Lambda][[n]]^2*t}, {n, 1, Nmax}];

Plot[u[x, MyTime], {x, 0, L}, PlotStyle -> {Red}]

View attachment 368
Let's produce our density plot

DensityPlot[u[x, y], {x, 0, 1}, {y, 0, 1}, PlotRange -> All, ColorFunction -> "Rainbow"]

View attachment 369
Next is our 3d plot

Plot3D[u[x, y], {x, 0, 1}, {y, 0, 1}, PlotRange -> All, Boxed -> False, ColorFunction -> "Rainbow"]

View attachment 370
Finally, we have the contour plot

ContourPlot[u[x, y], {x, 0, 1}, {y, 0, 1}, PlotRange -> All, ColorFunction -> "Rainbow"]

View attachment 371

 

[Dustinsfl]

Consider Laplace's equation $\nabla^2u = 0$ on the rectangle with the following boundary conditions:
$$
u_y(x,0) = f(x)\quad u(L,y) = 0\quad u(x,H) = 0\quad u(0,y) = g(y).
$$
Evaluate the series coefficients for the particular case of $f(x) = g(y) = 1$ and plot a contour map of $u(x,y)$ assuming $L = H = 1$.
So the solution is
$$
u(x,y) = \frac{4}{\pi}\sum_{n = 1}^{\infty}\left[\frac{\sin(2n - 1)\pi x\sinh\left[\pi(2n - 1) (1 - y)\right]}{\pi(2n - 1)^2\sinh\left[(2n - 1)\pi\right]} + \frac{\sin(2n - 1)\pi y\sinh\left[(2n - 1)\pi(1 - x)\right]}{(2n - 1)\sinh\left[(2n - 1)\pi\right]}\right].
$$
Let's use Mathematica to produce the desired plots.

Nmax = 40;
L = 1;
H = 1;
\[Lambda] = Table[(n*Pi)/L, {n, 1, Nmax}];

f[x_] = 1;
g[y_] = 1;A = Table[2/(L*\[Lambda][[n]]*Cosh[\[Lambda][[n]]*H])*Integrate[f[x]*Sin[\[Lambda][[n]]*x], {x, 0, L}], {n, 1, Nmax}];
B = Table[2/(L*Sinh[\[Lambda][[n]]*H])*Integrate[g[y]*Sin[\[Lambda][[n]]*y], {y, 0, H}], {n, 1, Nmax}];

u[x_, y_] = Sum[A[[n]]*Sin[\[Lambda][[n]]*x]*Sinh[\[Lambda][[n]]*(H - y)] + B[[n]]*Sin[\[Lambda][[n]]*y]*Sinh[\[Lambda][[n]]*(L - x)], {n, 1, Nmax}];

DensityPlot[u[x, y], {x, 0, L}, {y, 0, H}, PlotRange -> All, ColorFunction -> "Rainbow"]

View attachment 372

Plot3D[u[x, y], {x, 0, L}, {y, 0, H}, PlotRange -> All, Boxed -> False, ColorFunction -> "Rainbow"]

View attachment 373

ContourPlot[u[x, y], {x, 0, L}, {y, 0, H}, PlotRange -> All, ColorFunction -> "Rainbow"]

View attachment 374

 

[Dustinsfl]

Let's consider the problem we just did. What if $f(x)$ was defined

Nmax = 40;
L = 1;
H = 1;
\[Lambda] = Table[(n*Pi)/L, {n, 1, Nmax}];

f[x_] = Piecewise[{{1, x < L/2}, {0, x > L/2}}];
g[y_] = 1;A = Table[2/(L*\[Lambda][[n]]*Cosh[\[Lambda][[n]]*H])*Integrate[f[x]*Sin[\[Lambda][[n]]*x], {x, 0, L}], {n, 1, Nmax}];
B = Table[2/(L*Sinh[\[Lambda][[n]]*H])*Integrate[g[y]*Sin[\[Lambda][[n]]*y], {y, 0, H}], {n, 1, Nmax}];

u[x_, y_] = Sum[A[[n]]*Sin[\[Lambda][[n]]*x]*Sinh[\[Lambda][[n]]*(H - y)] + B[[n]]*Sin[\[Lambda][[n]]*y]*Sinh[\[Lambda][[n]]*(L - x)], {n, 1, Nmax}];

DensityPlot[u[x, y], {x, 0, L}, {y, 0, H}, PlotRange -> All, ColorFunction -> "Rainbow"]

View attachment 375

Plot3D[u[x, y], {x, 0, L}, {y, 0, H}, PlotRange -> All, Boxed -> False, ColorFunction -> "Rainbow"]

View attachment 376

ContourPlot[u[x, y], {x, 0, L}, {y, 0, H}, PlotRange -> All, ColorFunction -> "Rainbow"]

View attachment 377