Solving a PDE (2D Laplacian) coupled with an ODE

Question

I have asked this question before, but this is my new attempt and so instead of cluttering the previous one, I am making a new post. I am trying to analytically solve a PDE ($\nabla^2 T(x,y)=0$) coupled with an ODE. The PDE is subjected to the following boundary conditions:

$$\frac{\partial T(0,y)}{\partial x}=\frac{\partial T(L,y)}{\partial x}=0 \tag 1$$

$$\frac{\partial T(x,0)}{\partial y}=\gamma \tag 2$$

$$\frac{\partial T(x,l)}{\partial y}=\beta (T(x,l)-t) \tag 3$$

where $t$ is governed by the ODE:

$$\frac{\partial t}{\partial x}+\alpha(t-T(x,l))=0 \tag 4$$

subjected to $t(x=0)=0$. I am trying separation of variables. I manipulated $(4)$ to express $t$ as $t=\alpha e^{-\alpha x}\Bigg(\int_0^x e^{\alpha s }T(s,l)\mathrm{d}s\Bigg)$ and substituted in $(3)$ while applying the 3rd b.c.

My attempt is (I must acknowledge Bill Watts here as I have used methods I learned from his answer on MMA SE):

pde = D[T[x, y], x, x] + D[T[x, y], y, y] == 0

(*product form*)
T[x_, y_] = X[x] Y[y]
pde/T[x, y] // Expand
xeq = X''[x]/X[x] == -a^2
DSolve[xeq, X[x], x] // Flatten
X[x_] = X[x] /. % /. {C[1] -> c1, C[2] -> c2}
yeq = Y''[y]/Y[y] == a^2
DSolve[yeq, Y[y], y] // Flatten
Y[y_] = (Y[y] /. % /. {C[1] -> c3, C[2] -> c4})

(*addition form*)
T[x_, y_] = Xp[x] + Yp[y]
xpeq = Xp''[x] == b
DSolve[xpeq, Xp[x], x] // Flatten
Xp[x_] = Xp[x] /. % /. {C[1] -> c5, C[2] -> c6}
ypeq = Yp''[y] + b == 0
DSolve[ypeq, Yp[y], y] // Flatten
Yp[y_] = Yp[y] /. % /. {C[1] -> 0, C[2] -> c7}

T[x_, y_] = X[x] Y[y] + Xp[x] + Yp[y]
pde // FullSimplify

(*Applying the first and second b.c.*)
(D[T[x, y], x] /. x -> 0) == 0
c6 = 0
c2 = 0
c1 = 1

(D[T[x, y], x] /. x -> L) == 0
b = 0
a = (n π)/L
$Assumptions = n ∈ Integers

(*Applying the third b.c.*)
(D[T[x, y], y] /. y -> 0) == γ
c4 = c4 /. Solve[Coefficient[%[[1]], Cos[(π n x)/L]] == 0, c4][[1]]
c7 = c7 /. Solve[c7 == γ, c7][[1]]
T[x, y] // Collect[#, c3] &

(*now splitting T[x,y] into two parts*)
T[x, y] /. n -> 0
T0[x_, y_] = 2 c3 + c5 + y γ /. c5 -> 0
Tn[x_, y_] = T[x, y] - T0[x, y] // Simplify

(*applying the fourth b.c. to each part individually and using orthogonality*)

bcfn0 = (D[T0[x, y], y] /. y -> l) == β (T0[x, l] - α E^(-α x) Integrate[E^(α s) T0[s, l], {s, 0, x}])
Integrate[bcfn0[[1]], {x, 0, L}] == Integrate[bcfn0[[2]], {x, 0, L}]
Solve[%, c3]
c3 = c3 /. %[[1]]

bcfn = (D[Tn[x, y], y] /. y -> l) == β (Tn[x, l] - α E^(-α x) Integrate[E^(α s) Tn[s, l], {s, 0, x}])
Solve[Integrate[bcfn[[1]]*Cos[(n*Pi*x)/L], {x, 0, L}] == Integrate[bcfn[[2]]*Cos[(n*Pi*x)/L], {x, 0, L}], c5];
c5 = c5 /. %[[1]];//FullSimplify
T0[x_, y_] = T0[x, y] // Simplify
Tn[x_, y_] = Tn[x, y] // Simplify

Now we declare some constants and compile the functions

α = 62.9/2;
β = 1807/390;
γ = 3091.67/390;
L = 0.060;
l = 0.003;

T[x_, y_, mm_] := T0[x, y] + Sum[Tn[x, y], {n, 1, mm}]

Plot[{Evaluate[T[x, 0, 10]], Evaluate[T[x, l/2, 10]], Evaluate[T[x, l, 10]]}, {x, 0, L}]

The plot results are extremely ambiguous. The solution is not even converging (as I increase the number of terms , the T value keeps increasing). I cannot figure out what I have done wrong. Since the $T$ results are completely out, I have not calculated $t$. I cannot figure out what I have done wrong.

Answer 1

I can fix the problem of your solution increasing with increasing n, but that will not give you a solution. Rather than copy your entire solution I will start where I think the problem starts.

You have

T0[x_, y_] = 2 c3 + c5 + y γ /. c5 -> 0

Change that to

T0[x_, y_] = 2 c3 + c5 + y γ /. c3 -> 0
(*c5 + γ y*)

Then

Tn[x_, y_] = T[x, y] - T0[x, y] // FullSimplify
(*2 c3 Cos[(π n x)/L] Cosh[(π n y)/L]*)

In your case you were carrying an extra constant term c5 with Tn which was being added for every term in your sum which is why your solution increased with each term. In my case I carry c5 as the constant term, but only with T0. The changes below will require changing solving for c5 with bcf0 and solving for c3 with bcfn.

This next problem I fear is insurmountable with the calculation of bcfn0.

bcfn0 = (D[T0[x, y], y] /. y -> l) == β (T0[x, l] - α E^(-α x) Integrate[E^(α s) T0[s, l], {s, 0, x}]) // FullSimplify
(*γ E^(α x) == β (c5 + γ l)*)

Examining this result, it is obvious that there is no constant value c5 can take to satisfy this equation.

Furthermore, with the new Tn the orthogonality equation will result in c3 = 0. This means that T will have no x dependence, which when you think about it makes sense, if T is to satisfy Laplace eq and have x derivatives equal to zero at both ends in the x direction.

If T has no x dependence, then its derivatives can also have no x dependence, but with the y derivative of T depending on t which has x dependence, we have a problem.