Fdtxjr

The governing differential equation is

$$nabla^2 T=0 tag A$$

The boundary conditions for this problem are as foll0ws:

$$T(0,y,z)=T_hi tag 1A$$

$$T(L,y,z) = T(x,0,z) = T(x,l,z) = T(x,y,0)= T(x,y,mu) = 0 tag 1B$$

I need to solve for the distribution $T(x,y,z)$. I include my attempt here.

Solution attempt

Homogeneous Dirichlet type B.C. on $y$ and $z$ faces allow us to write the following form of preliminary temperature distribution:

$$T(x,y,z) = sum_n,m=1^inftysinbigg(fracnpi ylbigg)sinbigg(fracmpi zmubigg)T_nm(x) tag 1C$$

We substitute $(1mathrmC)$ in $mathrmA$ and apply the orthogonality properties of $sinbigg(fracnpi ylbigg)$ and $sinbigg(fracmpi zmubigg)$ to obtain the following:

$$fracmathrmd^2T_nm(x)mathrmdx^2 - gamma_nm^2 T_nm(x) = 0 tag 1D$$
where $gamma_nm = sqrt(fracnpil)^2+(fracmpimu)^2$

The general solution of $(1mathrmD)$ is of the form:

$$T_nm(x) = A_nme^gamma x + B_nme^-gamma x tag 1E$$

$A_nm,B_nm$ are the unknown Fourier coefficients that need to be determined. Applying $T(L,y,z) = 0$ on $(1mathrmD)$ and using $(1mathrmE)$, we arrive at:

$$A_nme^gamma L + B_nme^-gamma L = 0$$

$$bfB_nm = - A_nme^2gamma L tag 1F$$

$$T_nm(x) = A_nm(e^gamma x - e^2gamma L - gamma x) tag 1G$$

Using the B.C. $(1mathrmA)$:

$$T_hi = sum_n,m=1^inftysinbigg(fracnpi ylbigg)sinbigg(fracmpi zmubigg)A_nm(1 - e^2gamma L) tag 1H$$

Multiplying both sides of $(1mathrmH)$ by $int_0^lsinbigg(frackpi ylbigg)mathrmdy$ and $int_0^musinbigg(fracjpi zmubigg)mathrmdz$ and using the principle of orthogonality, we arrive at:

$$T_hiint_0^l int_0^mu sinbigg(frackpi ylbigg)sinbigg(fracjpi zmubigg) mathrmdymathrmdz = fraclmu4 A_kj (1 - e^2gamma L) tag 1I$$

Solving the definite integrals involved , we arrive at (for any arbitrary integer $n$ and $m$):

$$bfA_nm = frac4 T_hinmpi^2(1-e^2gamma L) (1-cos(npi))(1-cos(mpi)) tag 1J$$

$$T(x,y,z) = sum_n,m=1^infty(A_nme^gamma x + B_nme^-gamma x)sinbigg(fracnpi ylbigg)sinbigg(fracmpi zmubigg) tag 1K$$

On coding this solution in MATLAB. and substituting $x=0$, I find that the answer is not equal to $T_hi$.
For parameter values $L=0.5,l=0.5,mu=0.05,T_hi=50$. When i evaluate for $x=0,y=l/2,z=mu /2$, the answer should be $T=50$, but it evaluates to $96$ using the first four terms in the series. The series surely converges.

Is there something wrong, in the way, I am doing the problem ?Any help is greatly appreciated.

An observation

When I take $z=mu /2, y=l/2$ and $x=0$, and consider only odd values for $n$ and $m$, the solution can be written as:

$$Tbigg(0,fracl2,fracmu2bigg) = frac16 T_hipi^2underbracebigg[1 + frac19 + frac125 + frac149 + ........bigg]_pi^2 /8 = 2T_hi tag 1L$$

So, is there something wrong with my analytical solution ?

Are you sure you typed it in correctly? I checked your solution against mine and they're the same. I'll leave my answer below in case it helps anyways.

Often times I find it better to start from the beginning when I'm stuck to remove doubt of errors made in assumptions.
begincases
beginalign
Delta u(x,y,z) &= 0 \
u(0,y,z) &= u_0 \
uvert_partialOmegabackslashx=0 &= 0.
endalign
endcases
You know about separation of variables already, so we may take the eigenfunctions in the $y$ and $z$ directions as $sin fracnpi yell$ and $sin fracmpi zmu$. The equation in $x$ is known as $a_mne^gamma_mnx + b_mne^-gamma_mnx.$ Thus $u$ is the resulting infinite linear combination
$$
u(x,y,z) = sum_m,n geq 1left(a_mne^gamma_mnx + b_mne^-gamma_mnxright)sinfracnpi yellsinfracmpi zmu.
$$
The boundary conditions on $y$ and $z$ are already satisfied so we look towards those for $x$.
begincases
beginalign
u_0 &= sum_m,n geq 1(a_mn + b_mn)sinfracnpi yellsinfracmpi zmu \
0 &= sum_m,n geq 1 left(a_mne^gamma_mnL + b_mne^-gamma_mnLright)sinfracnpi yellsinfracmpi zmu
endalign
endcases
Note that using orthogonality is not multiplying each side by the respective integrals; it is multiplying each side by the eigenfunctions with dummy indices and then integrating. Continuing with this in mind,
begincases
beginalign
int_0^ell int_0^mu u_0 sinfracnpi yell sinfracmpi zmu , mathrmdz ,mathrmdy &= fracmuell4(a_mn + b_mn) \
0 &= a_mne^gamma_mnL + b_mne^-gamma_mnL
endalign
endcases
The first integral is easily found and we have the system of equations
$$
beginpmatrix
1 & 1 \
e^gamma_mnL & e^-gamma_mnL
endpmatrix
beginpmatrix
a_mn \
b_mn
endpmatrix
=
beginpmatrix
frac4big((-1)^m - 1big) big((-1)^n - 1big)pi^2mnu_0 \
0
endpmatrix.
$$
Inverting to solve for $a_mn$ and $b_mn$ then,
beginalign
beginpmatrix
a_mn \
b_mn
endpmatrix
&=
-fracoperatornamecsch(gamma_mnL)2
beginpmatrix
e^-gamma_mnL & -1 \
-e^gamma_mnL & 1
endpmatrix
beginpmatrix
frac4big((-1)^m - 1big) big((-1)^n - 1big)pi^2mn \
0
endpmatrix
\
&= -frac2u_0pi^2fracbig((-1)^m - 1big) big((-1)^n - 1big)sinh(gamma_mnL) mn
beginpmatrix
e^-gamma_mnL \
-e^gamma_mnL
endpmatrix.
endalign
Thus $u$ is written
beginequation
u(x,y,z) = -frac4u_0pi^2 sum_m,n geq 1 Gamma_mn fracsinhbig(gamma_mn(x-L)big)sinh(gamma_mnL) sinfracn pi yell sin fracmpi zmu
endequation
with
$$
Gamma_mn = fracbig((-1)^m - 1big) big((-1)^n - 1big)mn.
$$

$hskip 1 in$

Mathematica code - it should be copy-pastable - enjoy.

ClearAll["Global`*"]
[Gamma]mn = Sqrt[((n [Pi])/[ScriptL])^2 + ((m [Pi])/[Mu])^2];
(* u0 not specified yet *)
L = 1;
[ScriptL] = 1;
[Mu] = 1;
(* Keep q low (~10-15 or so) as the number of terms grows as q^2 *)
(* i.e. higher q gives more accuracy but takes longer to run in effect *)
q = 5;

(*The coefficients amn and bmn amn,bmn=a,b/.Solve[a+b[Equal]4/
[Pi]^2(((-1)^m-1)((-1)^n-1))/(m n)u0,a Exp[[Gamma]mn L]+b Exp[-
[Gamma]mn L][Equal]0,a,b][[1]];*)

u[x_, y_, 
 z_] = (-4 u0)/[Pi]^2 Sum[
 Sum[(((-1)^m - 1) ((-1)^n - 1))/(m n) Sinh[[Gamma]mn (x - L)]/
 Sinh[[Gamma]mn L] Sin[n [Pi] y/[ScriptL]] Sin[
 m [Pi] z/[Mu]], n, 1, q], m, 1, q];

Plot3D[u[0, y, z]/u0, y, 0, [ScriptL], z, 0, [Mu], 
 PlotLabel -> 
 "Laplace's Equation Solutionnat x=0 withn[ScriptL]=1, [Mu]=1, 
25 terms", 
 AxesLabel -> "y", "z", "!(*FractionBox[(u), (u0)])", 
 Boxed -> False]

(* Something extra - it's computationally expensive (3d), at least for my 
computer, so let it run for a bit *)
Manipulate[
 ContourPlot[u[x, y, z]/u0, x, 0, L, y, 0, [ScriptL], 
 ColorFunction -> "DarkRainbow"], z, 0, [Mu]]

(* alt+.to quit if it takes too long & reduce q *)
ContourPlot3D[
 u[x, y, z]/u0, x, 0, L, y, 0, [ScriptL], z, 0, [Mu], 
 Mesh -> False, ColorFunction -> "Rainbow", 
 PlotLabel -> "Contours of Temperature", AxesLabel -> "x", "y", "z"]