Finish and plot this separation of variables problem

[Daniel Sellers]

TL;DR Summary

I have found general solutions to a PDE which I need help applying Dirichlet conditions to so I can plot a final solution.

I have a PDE which I have solved numerically using a guass-seidel method, but I want to compare it to the analytical solution. I have used separation of variables to get the general solution, but I need help applying it.

The PDE is

(1/s)⋅(∂/∂s)[(s/ρ)(∂ψ/∂s)] + (1/s2)⋅(∂/∂Φ)[(1/ρ)(∂ψ/∂Φ)] - 2Ω + ρ(c0 + c1ψ) = 0

and applying separation of variables gives (I think):

ψ = (1/ρc₁)⋅[AJ_k(ρs√c₁) + BY_k(ρs√c₁)]⋅[Ccos(kΦ) + Dsin(kΦ)] + (2Ω/ρc₁) + (c₀/c₁)

Where ρ, Ω, c₀, c₁ are known constants.

I am not sure what to do from here. I want to apply Dirichlet boundary conditions on an inner boundary s₀ and outer boundary s₁. Should these be series solutions? Like summations over k =1,2 ... and if so then how do I handle the coefficients (A,B,C,D)?

The domain is a polar annulus, so k must be an integer, but that's all I can readily figure out.

 

[archaic]

I suggest you write the math in LaTeX for legibility purposes.
https://www.physicsforums.com/help/latexhelp/

 

[Daniel Sellers]

Fair enough!:
The PDE:

$$ \frac{1}{s}⋅\frac{∂}{∂s}(\frac{s}{ρ}\frac{∂ψ}{∂s}) + \frac{1}{s^2}⋅\frac{∂}{∂Φ}(\frac{1}{ρ}\frac{∂ψ}{∂Φ}) - 2Ω + ρc_0 + ρc_1ψ = 0$$

What I think the general solution is:

$$ ψ = \frac{1}{ρc_1}[AJ_k(ρs\sqrt{c_1}) + BY_k(ρs\sqrt{c_1})]⋅[Ccos(kΦ) + Dsin(kΦ)] + \frac{2Ω}{ρc_1} + \frac{c_0}{c_1}$$

I am inexperienced with using bessel functions, and I admit I have forgotten a lot about solving PDEs, so I could really use some help/advice here.

 

[Daniel Sellers]

I will also elaborate on my recent attempts to solve this:

If I want the solution to be axisymmetric (which physically, I think I do), then it seems I need to require k = 0, then I can solve for A and B using inner and outer boundary conditions $(ψ(s_0) = ψ_{inner} , ψ(s_1) = ψ_{outer})$ with order zero for the bessel functions. This gives a solution, but I'm not confident that it is correct, and it requires me to choose A arbitrarily.

Or do I need to solve the Φ dependent part as Fourier series in order to determine A?

Still very confused, any help is appreciated

 

[archaic]

I will also elaborate on my recent attempts to solve this:

If I want the solution to be axisymmetric (which physically, I think I do), then it seems I need to require k = 0, then I can solve for A and B using inner and outer boundary conditions $(ψ(s_0) = ψ_{inner} , ψ(s_1) = ψ_{outer})$ with order zero for the bessel functions. This gives a solution, but I'm not confident that it is correct, and it requires me to choose A arbitrarily.

Or do I need to solve the Φ dependent part as Fourier series in order to determine A?

Still very confused, any help is appreciated

I will also advise you to ask this on math.stackexchange, people answer relatively quickly there.

 

[Fred Wright]

I find your problem interesting because I see a sort of analog with flux pinning in type II superconductors. For the constant $\rho$ problem I suggest you treat the problem analogous to the wave equation of a circular drum head, for example here: http://ramanujan.math.trinity.edu/rdaileda/teach/s14/m3357/lectures/lecture_4_3_slides.pdf. I treat the problem as being axially symmetric with a constant radial distribution fuction to evaluate the constants in the PDE. I get,$$ \psi(r)= 2 \rho_0 \sum_{n=0}^\infty \frac {J_0 (\lambda_{0m} r)} {\lambda_{0m} J_1 (\lambda_{0m})} +\rho_ 0 (\rho_0 c_0 - 2\Omega)$$ where $$\lambda_{mn}=\frac {\alpha_{mn} \rho_0^2 c_1} {R_{static}}$$ and $\alpha_{mn}$ is the $n^{th}$ positive zero of $J_m$, $R_{static}$ is the radius where the stream field goes to zero (approximately) and $\rho_0$ is the static constant mass density.

I have some questions about the physical interpretation of the PDE. I rewrite your PDE as $$\frac {1} {\rho} \nabla^2 \psi =2\Omega - \rho c_0 + \vec \nabla \frac {1} {\rho} \cdot \vec \nabla \psi - c_1 \psi =? \omega$$ where $\omega$ is (possibly?) the vorticity pointing along the z axis. Is $\Omega$ an externally applied angular velocity (like from a centrifuge motor or a constant magnetic field penetrating a type II superconductor)? Is this a vorticity wave equation for density coupled into the stream field? What functional form of $\rho(r,\theta)$ are you using? How was this equation derived? I didn't take fluid dynamics in college and I am confused by the physics behind the PDE:-)

 

[Daniel Sellers]

I find your problem interesting because I see a sort of analog with flux pinning in type II superconductors. For the constant $\rho$ problem I suggest you treat the problem analogous to the wave equation of a circular drum head, for example here: http://ramanujan.math.trinity.edu/rdaileda/teach/s14/m3357/lectures/lecture_4_3_slides.pdf. I treat the problem as being axially symmetric with a constant radial distribution fuction to evaluate the constants in the PDE. I get,$$ \psi(r)= 2 \rho_0 \sum_{n=0}^\infty \frac {J_0 (\lambda_{0m} r)} {\lambda_{0m} J_1 (\lambda_{0m})} +\rho_ 0 (\rho_0 c_0 - 2\Omega)$$ where $$\lambda_{mn}=\frac {\alpha_{mn} \rho_0^2 c_1} {R_{static}}$$ and $\alpha_{mn}$ is the $n^{th}$ positive zero of $J_m$, $R_{static}$ is the radius where the stream field goes to zero (approximately) and $\rho_0$ is the static constant mass density.

I have some questions about the physical interpretation of the PDE. I rewrite your PDE as $$\frac {1} {\rho} \nabla^2 \psi =2\Omega - \rho c_0 + \vec \nabla \frac {1} {\rho} \cdot \vec \nabla \psi - c_1 \psi =? \omega$$ where $\omega$ is (possibly?) the vorticity pointing along the z axis. Is $\Omega$ an externally applied angular velocity (like from a centrifuge motor or a constant magnetic field penetrating a type II superconductor)? Is this a vorticity wave equation for density coupled into the stream field? What functional form of $\rho(r,\theta)$ are you using? How was this equation derived? I didn't take fluid dynamics in college and I am confused by the physics behind the PDE:-)

There is a lot going into this derivation but it is intended to give the stream function and therefore velocity at equilibrium of a rotating system. Ω is simply frame rotation (taken to be known for the purposes of this question). Indeed the more compact form of this PDE is:

$$ξ_z = ρf(ψ) - 2Ω$$

where $ξ_z$ is the system vorticity. We don't have a function for ρ (except that in this case I am calling it constant in space) but in general will be found numerically from a different part of the algorithm.

So... you're saying there are 2 indices to take into account representing double sums over $m,n$ except that we set n = 0 to make it axisymmetric?

Can you explain more about where my boundary conditions come into play? For instance what if I want to apply ψ = 0 at the inner boundary of the annulus and ψ = 1 at the outer boundary ($R_{static}$?)? Is this possible or does the solution you've suggested only work with certain boundary conditions? If that's the case, it may still be useful as a test.

I will study the link you posted and see if I can answer this myself as well

Thanks for all your help so far!

 

[Daniel Sellers]

As I said I am not practiced with using bessel functions in any capacity (besides a few homework problems over a year ago) but I am getting the impression that they are difficult to work with unless you can state a boundary R(a) = 0. This had something to do with the fact that the bessel function zeroes have to be determined numerically?

Am I misunderstanding this? Is there a method for applying arbitrary Dirichlet BCs?

 

[Daniel Sellers]

Ok after studying that link I might be a little closer to making sense of this.

My domain does not include the origin. Should I discard the $Y_k$ anyway because the function needs to remain physical anyway?

In any case I should treat the bessel functions the same way I would treat sine or cosine functions in a Fourier series, and apply fouriers trick (or something analogous for bessel functions, since they have the same orthogonality as trig functions) at the boundaries.

Still not sure specifically how to do what I just said, but does this seems like I'm on the right track?

 

[Fred Wright]

My domain does not include the origin. Should I discard the YkY_k anyway because the function needs to remain physical anyway?

Yes, I think you are on the right track. I had a thought about modeling $\rho(r',\theta)$, $r'=\frac {r}{R_{static}}$. The centrifugal potential is$$\phi_{cent}=-\frac {1} {2} \Omega^2 r'^2 sin(\theta)^2$$ and the potential for a "rigid" 2d rotating fluid is$$\phi_{rf}=\frac {1} {2} \rho \Omega^2 r'^2$$I propose (naively) an effective potential,$$\phi_{eff}= -\frac {1} {2}\Omega^2 r'^2 sin(\theta)^2 + \Omega^3 r'^2 r'^{\Omega}$$This potential would tend to increase the density towards $R_{static}$ as $\Omega$ increases, and the density function becomes $$\rho(r,\theta)= \rho_0 -sin(\theta)^2 + \Omega r'^{\Omega} $$Just a thought and probably wrong but I don't know how you can model $\rho(r,\theta)$ without some physics to go by.

 

[Daniel Sellers]

As I understand it the density is going to be calculated based on the enthalpy once psi is known for each iteration of a larger algorithm.
So setting m = 0 and B = 0, and choosing ψ (s₁) = ψ_outer = 0 as a boundary condition then I can state the R(r) solution as a sum, with λ_0,n defined in terms of the zeros of each bessel function J_k...
This approach doesn't take the inner boundary condition into account, and I'm not sure I can carry it out properly, but that's a start!
This is a lot of work just to check one very specific case of this equation!