Partial differential equation with mixed boundary conditions

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.

• You need to edit a bit. Do you mean $u_1$ and $u_2$ and are they constant? Switching to polar coordinates looks natural first, then possibly a Laplace transform in the $r$ variable.
– Paul
Apr 16 '21 at 9:43
• You are right, maybe it would be better stated as: $$\nabla^2u=k\text{}\frac{\partial u}{\partial x}$$ with $u(x,y,z)$ and boundary conditions: $$u(x,y,z)=u_1 \quad \text{ for } \quad x^2+y^2+z^2=r_1^2$$ $$u(x,y,z)=u_2 \quad \text{ for } \quad x^2+y^2+z^2\to \infty$$ with $u_1,u_2, \text{ and }r_1 \text{ constants}$
– umby
Apr 16 '21 at 10:21
• In spherical coordinate it is: $$k\left(\sin \theta \cos \phi \frac{\partial u}{\partial r}+\frac{\cos \theta \cos \phi }{r}\frac{\partial u}{\partial \theta }\right)=\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)$$ with boundary conditions: $$u=u_1\text{, }r=r_1$$ $$u=u_2\text{, }r\to \infty$$
– umby
Apr 16 '21 at 11:03

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$$).

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.

• Dear Juan Pablo Vesga, thank you very much for your answer. My problem is the mathematical formulation of a physical problem. In Cartesian coordinate, k is the constant x-component of a vector. Further, $u$ does not depend on $\phi$, so that the final equation should be the one I gave below.
– umby
Apr 16 '21 at 13:24
• Could you kindly guide me to what you think is the best approach? Thanks so much.
– umby
Apr 16 '21 at 14:34
• Hi! A couple of points: 1. Could you state the physical problem you are trying to solve? I feel it's something related with fluids, but it would be good to have some context 2. The fact that $u$ does not depend on $\phi$ is something you should have stated in your starting question. A full statement of the problem (where, for example, some kind of rotational symmetry becomes evident) would be very good. Let me know that and I'll try to give you some more hints Apr 16 '21 at 15:07
• Hi, the problem regards the conductive transfer of heat from an isotherm sphere in a moving un-viscid fluid. The flow field has just one constant component, along x for example; steady-state condition.
– umby
Apr 17 '21 at 9:29
• The problem in spherical coordinate is: $$\sin \theta \cos \phi \frac{\partial u}{\partial r}+\frac{\cos \theta \cos \phi }{r}\frac{\partial u}{\partial \theta }=\frac{\alpha}{v} \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 ]$$ with boundary conditions: $$u=u_1\text{, }r\leqslant r_1$$ $$u=u_2\text{, }r\to \infty$$ $\alpha$ and $v$ are constants. $u$ s expected to be independent of $\phi$
– umby
Apr 17 '21 at 13:44