Partial differential equation with mixed boundary conditions

Question

I am trying to find the steady-state temperature field in a semi-infinite solid on whose surface there is an isotherm spherical cap sunken by a length $p$ (the sphere from which the cap is derived has radius $R$). The temperature of the cap is $T_1$, while temperature far from the surface is $T_0$.

The solid is beneath the surface shown in the figure above. Except for the cap, the rest of the surface is adiabatic. Normalizing, length dimensions were re-scaled by means of $R$, this means that $p$, the depth of the cap below the surface, can vary from 0 to 1 (0 = no cap, 1 = hemisphere). The dimensionless temperature is: $\left (T-T_0 \right )/\left (T_1-T_0 \right )$, then $T_1$ in dimensionless form is $1$ and $T_0$ becomes $0$. In a spherical coordinate system (physical convention, polar axis perpendicular to the surface) centered at the "center" of the cap, the dimensionless equations should be: $$\frac{\partial }{\partial r}\left(r^2\frac{\partial T}{\partial r}\right)+\frac{1}{\sin \theta }\frac{\partial }{\partial \theta }\left(\sin \theta \frac{\partial T}{\partial \theta }\right)=0$$ B.C

$T=0$, $r\to \infty$ ($T$ equal to the initial one far from the cap),

$T=1$, $r=\frac{1}{2} \left(\sqrt{2} \sqrt{-p^2+(1-p)^2 \cos (2 \theta )+2 p+1}+2 (1-p) \cos (\theta )\right)$ ($T$ on the cap, this is the equation of the cap in the chosen reference system),

$\frac{\partial T}{\partial \theta}\bigg| _{\theta=\pi/2}=0$ (adiabatic condition),

$\frac{\partial T}{\partial \theta}\bigg| _{\theta=\pi}=0$ (symmetry conditions). The problem is quite important in spot welding applications. Following the advice of @Juan Pablo Vesga I tried to separate the equation imposing $T=f(r)g(\theta)$. This yields: $$\frac{1}{f}\frac{\partial }{\partial r}\left ( r^{2} \frac{\partial f}{\partial r}\right )=-k$$ and $$\frac{1}{sin\theta}\frac{1}{g}\frac{\partial }{\partial \theta}\left ( sin\theta \frac{\partial g}{\partial \theta}\right )=k$$ Applying the BC is quite complicate at this point, however. I am quite stack and any help will be appreciated.

Answer 1

Your problem, as other people commenting have noted, seems very much suited for using spherical coordinates. Defining these as (notice that I use the "physicist convention" for the polar angle)

$$ x = r\cos \phi \sin \theta, \\ y = r\sin \phi \sin \theta,\\ z = r\cos \theta $$

You can write the Laplacian (you can look into the derivation here conversion of laplacian from cartesian to spherical coordinates) as

$${\nabla}^2u = \frac{1}{r^2} \frac{\partial}{\partial r}\left(r^2\frac{\partial u}{\partial r} \right) + \frac{1}{r^2 \sin\theta}\frac{\partial}{\partial \theta}\left( \sin \theta \frac{\partial u}{\partial \theta}\right) + \frac{1}{r^2\sin^2 \theta} \frac{\partial^2 u}{\partial \phi^2}$$

But, in this case, you also need to rewrite the partial derivative with respect to $x$ in terms of the partial derivatives of $u$ with respect to $r,\theta,\phi$. Using the chain rule,

$$\frac{\partial u}{\partial x} = \frac{\partial r}{\partial x}\frac{\partial u}{\partial r} + \frac{\partial \theta}{\partial x}\frac{\partial u}{\partial \theta} + \frac{\partial \phi}{\partial x}\frac{\partial u}{\partial \phi}$$

You can invert the relations defining the spherical coordinates to obtain

$$ r = \sqrt{x^2 + y^2 + z^2}, \\ \theta = \arccos \left(\frac{z}{\sqrt{x^2+y^2+z^2}}\right), \\ \phi = \arctan \left(\frac{y}{x}\right) $$

to compute the partial derivatives (which you should check for yourself)

$$ \frac{\partial r}{\partial x} = \sin \theta \cos \phi, \\ \frac{\partial \theta}{\partial x} = \frac{\cot \theta}{r}, \\ \frac{\partial \phi}{\partial x} = -\frac{\sin \theta}{r\sin \theta} $$

and then, you obtain the PDE in spherical coordinates:

$$ \frac{1}{r^2} \frac{\partial}{\partial r}\left(r^2\frac{\partial u}{\partial r} \right) + \frac{1}{r^2 \sin\theta}\frac{\partial}{\partial \theta}\left( \sin \theta \frac{\partial u}{\partial \theta}\right) + \frac{1}{r^2\sin^2 \theta} \frac{\partial^2 u}{\partial \phi^2} - k\left(\sin \theta \cos \phi \frac{\partial u}{\partial r} + \frac{\cot \theta}{r} \frac{\partial u}{\partial \theta} - \frac{\sin \theta}{r\sin \theta} \frac{\partial u}{\partial \phi}\right) = 0 $$

and the boundary conditions become $u(r = r_1) = u_1$ and $u(r \to \infty) = u_2 $. How would you go about finding the general solution to the differential equation? Is there any approach you would initially take?

EDIT:

Let's assume that the PDE you need to solve is indeed

$$ k\left[ \frac{1}{r^2} \frac{\partial}{\partial r}\left(r^2\frac{\partial u}{\partial r} \right) + \frac{1}{r^2 \sin\theta}\frac{\partial}{\partial \theta}\left( \sin \theta \frac{\partial u}{\partial \theta}\right) \right] = \cos \theta \frac{\partial u}{\partial r} - \frac{\sin \theta}{r} \frac{\partial u}{\partial \theta} $$

The next step you'd like to take is to see if the equation is separable, which I will illustrate with an example. Consider the equation

$$ \frac{\partial u}{\partial x} - \frac{\partial^2 y}{\partial y^2} = 0 $$

with boundary condition $u(x,0) = u(x,L) = 0$ and let's try to write the solution in the form $u(x,y) = f(x)g(y)$. The equation would then become

$$ g(y)\frac{df}{dx} - f(x)\frac{d^2g}{dy^2} = 0 $$

(notice that the partial derivatives become ordinary derivatives, since each function depends on a single variable). Now, divide through by $f(x)g(y)$ and put one term on the RHS of the equation; that is,

$$ \frac{1}{f(x)}\frac{df}{dx} = \frac{1}{g(y)}\frac{d^2g}{dy^2} $$

Now, you have an equation between two functions of different arguments, which must hold for any value of $x,y$. That means that if you, for example, held $x$ constant and started varying $y$, equality must still hold. However, if the expression $\frac{1}{g(y)}\frac{d^2g}{dy^2}$ varied as $y$ changes while we hold the LHS to a single, fixed value, then the equality would be lost. This means that the only way to satisfy this equation is if each side of the equation is equal to a constant! Let me pick it in the following way:

$$ \frac{1}{g(y)}\frac{d^2g}{dy^2} = -k^2 $$

which implies that also

$$ \frac{1}{f(x)}\frac{df}{dx} = -k^2 $$

(You can pick the constant to be anything you want since its actual value would be fixed by the boundary conditions. In this case, I picked it as $-k^2$ because I know what kind of solutions I expect). Now, These equations are easy to solve: The first one is given in terms of sines and cosines

$$ g(y) = A\sin{ky} + B\cos{ky} $$

while the second one is an exponential

$$ f(x) = Ce^{-k^2 x} $$

with $A, B, C$ constants of integration. The boundary condition $u(x,0) = 0$ implies $g(0) = 0$, which means that $B = 0$. The second boundary condition implies $g(L) = 0$, which gives the relation

$$ A\sin{kL} = 0 $$

Now, this either means that $A = 0$ or $\sin{kL} = 0$. The first solution is not very interesting, since it would simply mean that $u(x,y) = 0$. The second one, though, is more interesting: since the sine vanishes at any value of the form $n\pi$ where $n$ is any whole number, what we get is an infinite (but countable) set of values for k:

$$ k_n = \frac{n\pi}{L} $$

where $n= 1,2,3,...$ up to infinity. Thus, we have a set of "fundamental" solutions

$$ u_n(x,y) = A_n e^{-\frac{\pi^2 n^2}{L^2}x}\sin{\left(\frac{n\pi}{L}y\right)} $$

and, since our equation is linear, any linear combination of these solutions will itself be a solution. This means that

$$ u(x,y) = \sum_{n=1}^{\infty}A_n e^{-\frac{\pi^2 n^2}{L^2}x}\sin{\left(\frac{n\pi}{L}y\right)} $$

Know, the constants $A_n$ you fix using additional boundary conditions. Let's say that, for example, you had the condition $u(0,y) = h(y)$ for some arbitrary function of $y$. Then, you would have the relation

$$ h(y) = \sum_{n=1}^{\infty}A_n \sin{\left(\frac{n\pi}{L}y\right)} $$

which looks difficult, but it is enough to notice that this equation is giving you the Fourier sine expansion of $h$, and then finding $A_n$ would amount to computing the coefficients of the Fourier sine expansion of $h(y)$ (If you've never done this or are a bit rusty, you can look it up in this link https://tutorial.math.lamar.edu/classes/de/FourierSineSeries.aspx).

In the case of your equation, try to write $u(r,\theta) = f(r)g(\theta)$ and see if you can achieve a form where each side is dependent only on $r$ or $\theta$. The solution to each of the resulting equations will be more involved than the simple case I just wrote down, most likely involving Legendre polynomials in $\theta$, but the general idea is the same. If you get stuck, reply again and I'll try and help you a bit more.

Answer 2

Partial solution

We will find a fundamental solution of the equation (here $r=|\vec r|=\sqrt{x^2+y^2+z^2}$, $Z$ axix is chosen as polar, and $z=r\cos\theta$) $$\nabla^2G(\vec r)-p\text{}\frac{\partial G}{\partial z}(\vec r)=\delta(\vec r)$$ Using FT we can write: $$G(\vec r)=\frac{1}{(2\pi)^3}\int_{R_k^3} d^3\vec k \,\hat G(\vec k)\,e^{-i\vec k\vec r}$$ The equation for Fourier transform of the function is $$(-k_x^2-k_y^2-k_z^2+ipk_z)\,\hat G(\vec k)=1$$

$$G(\vec r)=-\frac{1}{(2\pi)^3}\int d^3\vec k\frac{e^{-i\vec k\vec r}}{k_x^2+k_y^2+k_z^2-ipk_z}=-\frac{1}{(2\pi)^3}\int_0^{\infty}ds\int d^3\vec ke^{-i\vec k\vec r}e^{-s(k_x^2+k_y^2+k_z^2-ipk_z)}$$ $$=-\frac{1}{(2\pi)^3}\int_0^{\infty}ds\int dk_ze^{-s(k_z^2-ik_z(sp+z)\frac{1}{s}-(sp+z)^2\frac{1}{4s^2}+(sp+z)^2\frac{1}{4s^2})}\int dk_xe^{-s(k_x^2-ik_x\frac{1}{s}-x^2\frac{1}{4s^2}+x^2\frac{1}{4s^2})}\int dk_y(..)$$ $$=-\frac{\pi^{\frac{3}{2}}}{(2\pi)^3}\int_0^{\infty}e^{-\frac{x^2}{4s}-\frac{y^2}{4s}-\frac{(sp+z)^2}{4s}}s^{-\frac{3}{2}}ds=-\frac{2\pi^{\frac{3}{2}}e^{-\frac{pz}{2}}}{(2\pi)^3}\int_0^{\infty}e^{-(\frac{x^2}{4}+\frac{y^2}{4}+\frac{z^2}{4})t^2}e^{-\frac{p^2}{4t^2}}dt$$ $$=-\frac{e^{-\frac{pz}{2}}}{4\pi^{\frac{3}{2}}}\int_0^{\infty}e^{-\frac{r^2}{4}t^2-\frac{p^2}{4t^2}}dt=-\frac{e^{-\frac{pz}{2}}}{2\pi^{\frac{3}{2}}r}\int_0^{\infty}e^{-t^2-\frac{c^2}{t^2}}dt\,,\text{where } c^2=\frac{p^2r^2}{16}$$

$$I=\int_0^{\infty}e^{-t^2-\frac{c^2}{t^2}}dt=\frac{e^{-2|c|}}{2}\sqrt\pi$$ $$G(\vec r)=-\frac{1}{4\pi}\frac{\exp\bigl(-\frac{|p|\,r}{2}-\frac{p\,r\cos\theta}{2}\bigr)}{r}$$ At $p=0$ we get the fundamental solution for Poisson's equation $G(\vec r)=-\frac{1}{4\pi\,r}$

Fundamental solution allows to find the solution of the equation with non-zero RHS:

if $$\nabla^2u-p\text{}\frac{\partial u}{\partial z}=f(\vec r)$$ $$u(\vec r)=\int_{R^3}G(\vec r-\vec r')f(\vec r')d^3\vec r'$$

It is not clear, though, how we can find a solution in the case if the the function on the boundary (on the surface of spheres) is constant (i.e $u(a,\theta)=u_1$ and $u(b,\theta)=u_2$).