# A really difficult heat-diffusion problem

A very long beam with a quadratic cross-section $$0\le x\le L$$, $$0\le y\le L$$ lies on a plane with a temperature of 0 degrees. The other surfaces of the beam are in contact with air which holds 10 degrees temperature. Determine the stationary temperature distribution in the beam, assuming it is of infinite length $$-\infty\le z\le \infty$$.

The solution:

I am sure that there are two ways to solve this, one using the heat equation, and one using the Laplace equation.

The heat equation requires a time-dimension, which is not given, however they term "stationary" says it all. $$\frac{\partial}{\partial t}=0$$. So the heat equation for this should be:

$$\nabla^2 u=0$$

which in fact becomes a Laplace equation, where $$u=u(x,y,z)$$.

There is a problem with ansatz, it is not easy to make an ansatz here, but I shall try. Consider that the x and y dimensions at the origin are both in contact with a zero-temperature surface and a 10-degree temperature air, then we could consider both to be non-zero at the origin, if we use Kelvin degrees. So we assume both are cosinoid, and therefore have the ansatz:

$$u(x,y,z)=\cos\frac{n\pi}{L}x\cos\frac{m\pi}{L}y\cdot u(z)$$

The operator is then:

$$\nabla^2\bigg(\cos\frac{n\pi}{L}x\cos\frac{m\pi}{L}y\cdot u(z)\bigg):=-\frac{\cos\frac{m\pi x}{L}\cos\frac{n\pi x}{L}\big(\pi^2(m^2+n^2)\big)}{L^2}u(z) -u_{zz}$$ This gives

Insert in the PDE:

$$u_{zz} +\frac{\cos\frac{m\pi x}{L}\cos\frac{n\pi x}{L}(\pi^2(m^2+n^2)\big)}{L^2}u(z) =0$$

Let $$f(x,y)=\frac{\cos\frac{m\pi x}{L}\cos\frac{n\pi x}{L}\big(\pi^2(m^2+n^2)\big)}{L^2}$$, then we have the "solution":

$$u(z)=C_1\cos (f(x,y)\cdot z)+C_2\sin(f(x,y)\cdot z)$$

Insert in the original ansatz and expand:

$$u(x,y,z)=\cos\frac{n\pi}{L}x\cos\frac{m\pi}{L}y\cdot \big(C_1\cos \bigg(\frac{\cos\frac{m\pi x}{L}\cos\frac{n\pi x}{L}\big(\pi^2(m^2+n^2)\big)}{L^2}\cdot z\bigg)+C_2\sin\bigg(\frac{\cos\frac{m\pi x}{L}\cos\frac{n\pi x}{L}\big(\pi^2(m^2+n^2)\big)}{L^2}\cdot z\bigg)$$

Looking only at a level curve of $$u(x,y,z)$$, with $$z=0$$, it looks like this

But the function $$u(x,y,z)$$ is weird, because each time a level curve is made with x or y dimension equal to zero, then the z-term becomes 1.

In fact, considering a similar ODE problem:

$$y''+\cos xy'=0$$ gives a rather exotic solution, so solving the given PDE must be impossible, analytically.

Does anyone have an idea on how to solve this problem?

Thanks

• I think you should solve the equation including time with boundary condition at time $t = 0$ of 0 degree on the surface and 10 degrees in the air and in the final solution let $t \rightarrow \infty$. You are right in the sense $\partial/\partial t = 0$ but the boundary conditions given are at time $t=0$, you need to incorporate that ? So i think you need to solve the full heat equation ? or somehow deduce/find boundary conditions at $t = \infty$ and solve the laplace equation. Dec 5, 2022 at 10:22
• Further the problem statement seems symmetric w.r.t $z$, so shouldn't $u(z) = constant$ ? Dec 5, 2022 at 10:28
• The plane on which the beam lies is the $x$-$z$ plane or the $y$-$z$ plane. In either case, the solution has no $z$ dependence in the steady state, so you are solving $\nabla^2 u(x,y)=0$ in the square region $0<x<L$, $0<y<L$. The boundary conditions are $u=u_0$ along one side of the square and $u=u_1$ on the other three sides. Let $u(x,y)=v(x,y)+u_1$ so that $v$ satisfies Laplace's equation with BCs $v=u_0-u_1$ on one side and $v=0$ on the three others. Make the separation ansatz $u(x,y)=X(x)Y(y)$, then the rest should be 'standard'
– Sal
Dec 5, 2022 at 16:24
• I suppose so- the transient behavior could be $z$ dependent if the initial condition is but since you're asking about the steady state that's not necessary at all
– Sal
Dec 5, 2022 at 19:25
• @Sal see my suggestion Dec 6, 2022 at 8:56

The BVP is well posed and the solution may be represented as an infinite sum. The mistake is in the ansatz $$u(x,y)=u(y)\sin(n\pi x/L)$$, because the BCs at $$x=0$$ and $$x=L$$ are non-homogeneous.

Let $$u(x,y)=v(x,y)+C$$, with $$C=10$$. Now $$v$$ also satisfies Laplace's equation in the square but with the simpler BCs

$$\tag{1} v(x,0)=-C\qquad,\qquad v(x,L)=v(0,y)=v(L,y)=0$$

Solving (1) by separation of variables $$v(x,y)=X(x)Y(y)$$ is straightforward see eg here . The result is

$$\tag{2} v(x,y)=\sum\limits_{n=1}^\infty c_n\sin(n\pi x/L)\sinh(n\pi-n\pi y/L)\\ c_n=\frac{2}{L\sinh(n\pi)}\int\limits_0^L dx \ \sin(n\pi x/L)v(x,0)$$

With $$v(x,0)=-C$$ we find

$$\tag{3} c_n=-\frac{2C}{n\pi \sinh(n\pi)}\left(1-(-1)^n\right)$$

Inserting (3) into (2) we have

$$\tag{4} v(x,y)=-4C\sum\limits_{n=1,3,5\cdots}\frac{\sinh(n\pi-n\pi y/L)}{n\pi\sinh(n\pi)}\sin(n\pi x/L)$$

Before plotting, it's useful to switch to dimensionless variables. Define $$\chi=x/L$$, $$\eta=y/L$$, $$\phi=u/C$$. The BVP is now posed on $$0<\chi<1$$, $$0<\eta<1$$, with BCs $$\phi=1$$ and $$\phi=0$$. The solution is

$$\tag{5} \phi(\chi,\eta)=1-\sum\limits_{n=1,3,5,\cdots} \frac{4 \sinh(n\pi-n\pi \eta)}{n\pi \sinh(n\pi)}\sin(n\pi \chi)$$

Here is a plot of the $$n=1$$ term (left) and the partial sum to $$11$$ terms (right)

• reg 2, should it not b $\sinh((n\pi-n\pi) y/L)\\$? tnks Sal Dec 7, 2022 at 19:05
• PS the plot with N=1 and L=1 is exactly the same as mine , I add link a the end of your answer Dec 7, 2022 at 19:10
• @Luthier415Hz I think the argument of $\sinh$ should be dimensionless, so I don't think there is a typo there. Following the link it would be something like $\sinh\left(n\pi (H-y)/L\right)$ which reduces to (2) when $L=H$
– Sal
Dec 8, 2022 at 13:05
• Hi Sal thanks. Did you try to plot out your result with n=L=1? It is the same as my plot! Dec 8, 2022 at 13:06

By Sals suggestion, I add the proposed solution.

Let the boundaries of the heat-diffusion problem be defined as

We ignore the z-dimension.

The PDE is given therefore by:

$$\nabla^2 u=0 \ \ \ \ \ \ 0\le x\le L, \ \ \ \ 0\le y\le L\\ u(0,y)=10,\ \ \ \ u(L,y)=10 \ \ \ \ \ (1)\\ u(x,0)=0, \ \ \ \ u(x,L)=10 \ \ \ \ \ \ (2)$$

The IC in (1) are Dirichlet homogeneous on x, thus a suitable ansatz is $$u(x,y)=\sin\frac{n\pi x}{L}u(y)$$

We insert in the PDE and get:

$$-\bigg(\frac{n\pi}{L}\bigg)^2\sin\frac{n\pi x}{L}u(y)+\sin\frac{n\pi x}{L}u_{yy}=0$$

which gives the ODE:

$$u(y)_{yy}-\bigg(\frac{n\pi}{L}\bigg)^2u(y)=0$$

This gives the solution $$u(y)=C_1\cosh\frac{n\pi}{L}y+C_2\sinh\frac{n\pi}{L}y$$

Use I.C. from (2) and obtain:

$$u(y)=\frac{10}{\sinh n\pi}\sinh\frac{n\pi}{L}y$$

Then, the final solution is:

$$u(x,y)=\frac{10}{\sinh n\pi}\sinh\frac{n\pi}{L}y\sin\frac{n\pi}{L}x$$

Which has the plots

Any comment or eventual improvements are welcome!

• It is strange that it is ill-posed, because this is an exercise from a mathematical physics exam. How could the Prof miss so big? Dec 6, 2022 at 11:14
• (detypoed) Your model is ill-posed: What if $(x,y)$ is $(0,0)$ or $(L,0)$? The problem can be remodelled in a well-posed form: Cut off the problematic corners, for example by radiusing them with $\epsilon$-radius quarter-circles, along which the temperature rises from $0^\circ$ to $10^\circ$ (with $O(1/\epsilon)$ steepness!) according to some assumed profile. The temperature variation will get $(1/\epsilon)$-wild near those corners, and be highly dependent on the assumed temperature profile, but fairly stable elsewhere. I doubt that this has a tractable analytical solution. Dec 6, 2022 at 14:27
• Please try to provide a solution, I will certainly delete mine once yours is posted. Dec 6, 2022 at 15:45
• Sorry, but this is a question more for engineers than mathematicians. First you need to put in an intuitive profile for the rounded corners. Then run a computer software package on it. (I'm sure there is such a package for numerically solving heat-flow problems.) From a strictly mathematical perspective, the question is invalid. Actually, a physically realistic solution would be even harder, because passive contact with air at $10^\circ$ would not lift the surface of the beam to that temperature if the beam was cooled by continuous solid contact across one face. Dec 6, 2022 at 19:31
• Note that your solution violates equation 1: namely, the value of the function along the faces $x=0$ and $x=L$ is $0$, not $10$. Also, its value along $y=L$ violates equation 2b: it varies between $0$ and $10$ rather than holding the constant value $10$. It's always recommended, and usually easy, to check whether a proposed solution for a PDE satisfies the given boundary conditions. Dec 6, 2022 at 20:04