I've been trying to solve this problem with separation of variables where T is a function of r and z
$$\frac{1}{r}\frac{\partial }{\partial r}(r\frac{\partial T}{\partial r})+\frac{\partial^2T}{\partial z^2}=0$$
$$\left.-k\frac{\partial T}{\partial r}\right\rvert_{r=R}=h[T(R,z)-T_{\infty}]$$
$$\left.-k\frac{\partial T}{\partial z}\right\rvert_{z=H}+h[T(r,H)-T_{\infty}]=q_s$$
$$T(r,0)=T_0$$
$$\left.-k\frac{\partial T}{\partial r}\right\rvert_{x=0}=0$$
$where \ T_0 \ ,T_{\infty} \ and \ q_s \ are \ constants$
I did a change of variable $\theta(r,z) = T(r,z)-T_{\infty}$ which resulted in
$$\frac{1}{r}\frac{\partial }{\partial r}(r\frac{\partial \theta}{\partial r})+\frac{\partial^2 \theta}{\partial z^2}=0$$
$$\left.-k\frac{\partial \theta}{\partial r}\right\rvert_{r=R}=h\theta$$
$$\left.-k\frac{\partial \theta}{\partial z}\right\rvert_{z=H}+h\theta=q_s$$
$$\theta(r,0)=\theta_0$$
$$\left.-k\frac{\partial \theta}{\partial r}\right\rvert_{x=0}=0$$
I've been trying to solve it but the non-homogeneous BC, the one that equals to $q_s$ has been giving me problems. If I use $\theta (r,z)=G(r)M(z)$ I get
$$G(r)=C_1J_0(\lambda r)$$ and $$M(z)=C_3\cosh(\lambda z)+C_4\sinh(\lambda z)$$
To find $C_3$ and $C_4$ becomes a problem, if $q_s$ wasn't there I could express $C_3$ in terms of $\lambda , C_4$ and other functions. I've solved problems where there is a similar condition as the one with $q_s$ but T is a function of x and time $T(x,t)$ and I can apply something called "transient temperature" which is more or less a change of variable that looks like $T(x,t)=u(x,t)+v(x)$. Here I'm just stuck.
In case it is important, the equation and boundary conditions come from a two dimensional heat flow problem, where a cylinder with a radius R and a height H gets applied a uniform heat flux $q_s$ on the top. At the same time the cylinder undergoes convection and one of its faces is in contact with a surface at $T_0$, i.e. the cylinder has one of its "circular faces" resting against a surface while all the uncovered faces experience convection, and the other uncovered "circular face" has a uniform heat flux applied to it.
I have to find the temperature distribution i.e. $\theta (r,z)$
Hope I'm making myself clear, I'm just trying to get help.
Best Answer
Your general idea is good. Let's start with a statement of the problem.
Solve
$$\begin{cases} \frac{1}{r} \partial_r (r \partial_r T) + \partial_{z z} T = 0 & \text{subject to} \\ T (r, z = 0) = T_0 \\ -k \frac{\partial T}{\partial z} \big|_{z = H} = q_S'' + h (T (r, H) - T_\infty) \\ \frac{\partial T}{\partial r} \big|_{r = 0} = 0 \\ -k \frac{\partial T}{\partial r} \big|_{r = R} = h (T (R, z) - T_\infty) \end{cases}$$
Like you've done, the first step is to define a nondimensional temperature difference $\Theta = \frac{T - T_\infty}{T_0 - T_\infty}$ which helps to make the convection condition homogeneous. While we're at it, let's define non-dimensional variables $\hat{r} = \frac{r}{R} \in [0, 1]$, $\hat{z} = \frac{z}{H} \in [0, 1]$, $\hat{q}_S'' = \frac{q_S'' H}{k (T_0 - T_\infty)}$, ${\rm Bi} = \frac{h R}{k}$, and $\Pi = \frac{H}{R}$. Our system becomes
$$\begin{cases} \frac{1}{\hat{r}} \partial_{\hat{r}} (\hat{r} \partial_{\hat{r}} \Theta) + \frac{1}{\Pi^2} \partial_{\hat{z} \hat{z}} \Theta = 0 & \text{subject to} \\ \Theta (\hat{r}, \hat{z} = 0) = 1 \\ \frac{\partial \Theta}{\partial \hat{z}} \big|_{\hat{z} = 1} = -(\hat{q}_S'' + \Pi {\rm Bi} \Theta (\hat{r}, \hat{z} = 1)) \\ \frac{\partial \Theta}{\partial \hat{r}} \big|_{\hat{r} = 0} = 0 \\ \frac{\partial \Theta}{\partial \hat{r}} \big|_{\hat{r} = 1} = -{\rm Bi} \Theta (\hat{r} = 1, \hat{z}) \end{cases}$$
The boundary conditions in $\hat{r}$ here are homogeneous. Therefore, we may consider a solution by separation of variables:
$$\Theta = \sum_n g_n (\hat{z}) J_0 (\lambda_n \hat{r})$$
where $J_0$ are Bessel functions of the first kind. Because we know that $\frac{d }{d x} J_0 (x) = -J_1 (x)$, the eigenvalues satisfy
$$\lambda_n \frac{J_1 (\lambda_n)}{J_0 (\lambda_n)} = {\rm Bi}$$
Substituting this solution into the homogeneous governing equation, we find that for each mode $n$,
$$\left( -\lambda_n^2 g_n + \frac{1}{\Pi^2} g_n'' \right) J_0 (\lambda_n \hat{r}) = 0$$
which means that
$$g_n'' - (\lambda_n \Pi)^2 g_n = 0 \implies g_n = C_{1, n} \cosh (\lambda_n \Pi \hat{z}) + C_{2, n} \sinh (\lambda_n \Pi \hat{z})$$
where $C_{1, n}, C_{2, n}$ are unknown constants determined by the remaining boundary conditions in $\hat{z}$. Decomposing into the Bessel basis, we have that constant functions have coefficients
$$\langle J_0 (\lambda_n \hat{r}) | 1 \rangle = \frac{2}{\lambda_n} \frac{J_1 (\lambda_n)}{J_0^2 (\lambda_n) + J_1^2 (\lambda_n)}$$
Therefore, the boundary condition at $\hat{z} = 0$ gives
$$C_{1, n} = 2 \left( \frac{1}{\lambda_n} \right) \left( \frac{J_1 (\lambda_n)}{J_0^2 (\lambda_n) + J_1^2 (\lambda_n)} \right)$$
The boundary condition at $\hat{z} = 1$ can be rewritten as $\frac{\partial \Theta}{\partial \hat{z}} \big|_{\hat{z} = 1} + \Pi {\rm Bi} \Theta (\hat{r}, \hat{z} = 1) = -\hat{q}_S''$.
$$C_{1, n} (\lambda_n \Pi \sinh (\lambda_n \Pi) + \Pi {\rm Bi} \cosh (\lambda_n \Pi)) + C_{2, n} (\lambda_n \Pi \cosh (\lambda_n \Pi) + \Pi {\rm Bi} \sinh (\lambda_n \Pi)) = -\frac{2 \hat{q}_S''}{\lambda_n} \frac{J_1 (\lambda_n)}{J_0^2 (\lambda_n) + J_1^2 (\lambda_n)}$$
$$\implies C_{2, n} = -2 \left( \frac{1}{\lambda_n} \right) \left( \frac{\hat{q}_S'' + \lambda_n \Pi \sinh (\lambda_n \Pi) + \Pi {\rm Bi} \cosh (\lambda_n \Pi)}{\lambda_n \Pi \cosh (\lambda_n \Pi) + \Pi {\rm Bi} \sinh (\lambda_n \Pi)} \right) \left( \frac{J_1 (\lambda_n)}{J_0^2 (\lambda_n) + J_1^2 (\lambda_n)} \right)$$
The solutions above unfortunately aren't very numerically stable, but they do work while keeping a moderate number of modes. I plot below the solution for ${\rm Bi} = 1.0$, $\Pi = 1.0$, and $\hat{q}_S'' = -3.0$ (prescribed heat flux heats the system):
and for ${\rm Bi} = 1.0$, $\Pi = 1.0$, and $\hat{q}_S'' = 3.0$ (prescribed heat flux cools the system):