Stuck on solving a non homogenous (linear) PDE

[fluidistic]

I've got a PDE that I derived from a physical problem, so I suppose it has a solution and that it is unique. I am solving for streamlines in the region having a quarter of an annulus shape, so $\theta$ ranges between $0$ and $\pi/2$ and $r$ ranges between $r_i$ and $r_o$. The boundary conditions are that the streamlines have to be constant over the edges of the region, and this constant is arbitrary so we can take it as $0$ to simplify things.

The PDE is
$\frac{\partial ^2 \psi}{\partial \theta^2}+r^2\frac{\partial^2\psi}{\partial r^2}+r\frac{\partial \psi}{r}=A\frac{\cos(2\theta)}{r^2}$. EDIT: I made a typo, there is no 1/r^2 factor. As is, I cannot apply separation of variables due to the inhomogeneous term (the RHS). I think that the general solution is the one of the corresponding homogeneous PDE plus the one of the particular PDE involving this homogeneous term. The problem is that, I believe, that the solution to the homogeneous PDE is... 0, i.e. a vanishing solution. If this is true, then this means that I have to solve the PDE in one shot, i.e. an Ansatz or something like that, unless I switch to another solving method.
For information, I do know the solution to this PDE except that the RHS is $\frac{A\sin(2\theta)}{r^2}$. EDIT: same typo, no such 1/r^2 term here either. In that particular case, $\psi(x,y)= \psi_0 \left[1-\left( \frac{1}{r_i^2+r_o^2} \right) \left( r^2+\frac{r_i^2+r_o^2}{r^2} \right) \right] \sin(2\theta)$.In

My attempt at solving the homogeneous PDE using separation of variables leads to separable solutions of the form $\psi = \Theta(\theta)R(r)$ where the equation for $\Theta(\theta)$ is $\frac{\Theta ''}{\Theta}=-n^2$ with solutions $\Theta_m(\theta) = a_m \sin(2m\theta)$.
This makes the equation for: $r^2R''+rR'-4m^2R=0$ with solutions $R_m(r)=c_{1,m}\cos[2m\ln(r)]+ c_{2,m}\sin[2m\ln(r)]$. We can already see this does not match the solution to the non homogeneous PDE with the sine term, so I suppose $c_{1,m}=c_{2,m}=0$.
I tried the obvious Ansatz $\psi(x,y)= \psi_0 \left[1-\left( \frac{1}{r_i^2+r_o^2} \right) \left( r^2+\frac{r_i^2+r_o^2}{r^2} \right) \right] \cos(2\theta)$ but unfortunately this does not vanish at $\theta=0$ and $\theta=\pi/2$) while it should. Hmm, therefore I am running out of ideas to solve the original PDE. Any idea is greatly welcome.

 

[pasmith]

Set $\psi(r, \theta) = f(r) \cos 2\theta$. Then $$\frac{\partial^2 \psi}{\partial \theta^2} + r \frac{\partial \psi}{\partial r} + r^2 \frac{\partial^2 \psi}{\partial r^2} = (r^2 f'' + rf' - 4f)\cos 2\theta = Ar^{-2} \cos 2\theta$$ and hence $$r^2 f'' + rf' - 4f = Ar^{-2}.$$ To find the complementary functions, set $f(r) = r^\alpha$, giving $$(\alpha^2 - 4)r^\alpha = 0.$$ Thus $\alpha = \pm 2$. Since the right hand side of the non-homogenous problem is one of these complementary functions, we must modify it to find a particular solution. In this case we modify it by multiplying by $\ln r$, so that $$f(r) = Dr^2 + Er^{-2} + Cr^{-2} \ln r$$ where $C$ is chosen to satisfy the nonhomogenous equation and $D$ and $E$ are chosen to satisfy the boundary conditions.

I suspect, however, that you are trying to solve $$\nabla^2 \psi = \frac{\partial^2 \psi}{\partial r^2} + \frac 1r \frac{\partial \psi}{\partial r} + \frac{1}{r^2} \frac{\partial^2 \psi}{\partial \theta^2} = \frac{A}{r^2} \cos 2\theta$$ which after multipliying by $r^2$ becomes $$r^2 \frac{\partial^2 \psi}{\partial r^2} + r \frac{\partial \psi}{\partial r} + \frac{\partial^2 \psi}{\partial \theta^2} = A\cos 2\theta$$ leading to $$r^2 f'' + rf' - 4f = A.$$

 

[fluidistic]

Thanks a lot pasmith. I am verry sorry, I have made a typo mistake in the RHS, there should not be the $1/r^2$ factor. I am following your method but I realize I already tried it and got stuck.

Here it goes: assuming $\psi=f(r)\cos(2\theta)$, one reaches $f(r)=\frac{c_1}{r^2}+c_2r^2-\frac{A}{4}$. This is terrible because this implies what I wrote in the OP, namely that $\psi(r, \theta)=\left( \frac{c_1}{r^2}+c_2r^2-\frac{A}{4}\right)\cos(2\theta)$. Impossible to satisfy the boundary conditions. For example at $\theta=0$, psi must vanish for all r, meaning psi does not depend on r. And for r equals to rinner and router, psi must vanish for all theta, meaning psi doesn't depend on theta and that it's trivially 0 everywhere.

Do I miss something?

 

[renormalize]

Do I miss something?

Are you certain that you want $\psi=\text{const.}$ at the walls? Since the fluid velocity is $\vec{v}=\nabla\psi$ for potential flow, the normal (pun intended) boundary condition is that the normal-component of $\vec{v}$ vanishes at the walls. In your case that gives$$\left.\frac{\partial\psi}{\partial\theta}\right|_{\theta=0}=\left.\frac{\partial\psi}{\partial\theta}\right|_{\theta=\frac{\pi}{2}}=0$$and your solution$$\psi(r, \theta)=\left( \frac{c_1}{r^2}+c_2r^2-\frac{A}{4}\right)\cos(2\theta)$$satisfies this condition.

 

[fluidistic]

Are you certain that you want $\psi=\text{const.}$ at the walls? Since the fluid velocity is $\vec{v}=\nabla\psi$ for potential flow, the normal (pun intended) boundary condition is that the normal-component of $\vec{v}$ vanishes at the walls. In your case that gives$$\left.\frac{\partial\psi}{\partial\theta}\right|_{\theta=0}=\left.\frac{\partial\psi}{\partial\theta}\right|_{\theta=\frac{\pi}{2}}=0$$and your solution$$\psi(r, \theta)=\left( \frac{c_1}{r^2}+c_2r^2-\frac{A}{4}\right)\cos(2\theta)$$satisfies this condition.

Yeah I am sure.
Using your nomenclature, $v_r =\frac{1}{r}\frac{\partial \psi}{\partial \theta}$ and $v_\theta = - \frac{\partial \psi}{\partial r}$. When $\theta = 0$, $v_\theta = 0$, therefore it is $\frac{\partial \psi}{\partial r}$ that must vanish, not $\frac{\partial \psi}{\partial \theta}$. My solution does not satisfy this criterion.

 

[fluidistic]

...

Now I am wondering something. The curl of v changes sign at theta equals pi/4, so v has 2 swirls in opposite directions, along the sides of the quarter of the annulus. Now the streamlines ought to be constant over all sides, but I wonder if there should be a jump of this constant at theta equals pi/4. I think not, but maybe that's the only way to find a solution?
Any other idea?

 

[renormalize]

Using your nomenclature, $v_r =\frac{1}{r}\frac{\partial \psi}{\partial \theta}$ and $v_\theta = - \frac{\partial \psi}{\partial r}$. When $\theta = 0$, $v_\theta = 0$, therefore it is $\frac{\partial \psi}{\partial r}$ that must vanish, not $\frac{\partial \psi}{\partial \theta}$. My solution does not satisfy this criterion.

OK, my apologies, I confused your use of the fluid stream function $\psi$ with the potential function $\phi$.
So as I now understand it, you're interested in the 2D vorticity vector $\vec{\omega}\equiv\omega\left(r,\theta,A\right)\hat{e}_{z}\equiv\nabla\times\vec{v}=-\nabla^{2}\psi\,\hat{e}_{z}\,$, where $\omega\left(r,\theta,A\right)\equiv-\frac{A\,\cos2\theta}{r^{2}}$. Since it's easy to verify that $\nabla^{2}\omega\left(r,\theta,A\right)=0$, your stream function $\psi(r, \theta)=\left( \frac{c_1}{r^2}+c_2r^2-\frac{A}{4}\right)\cos2\theta$ actually satisfies the homogeneous 4th-order biharmonic equation $\nabla^{4}\psi=0$. This puts you squarely into the realm of Stokes flow, which requires two boundary conditions on each surface (http://brennen.caltech.edu/FLUIDBOOK/basicfluiddynamics/Stokesflow/equationsofstokesflow.pdf):

Because of this, satisfying your boundary conditions is possible, but it unfortunately requires a substantially more complicated stream function than the one shown above (http://brennen.caltech.edu/fluidbook/basicfluiddynamics/Stokesflow/planarflow.pdf):

Another useful reference to these equations is https://personalpages.manchester.ac...il/Lectures/Fluids/Material_2016/Chapter8.pdf:

Hope you find this helpful.

 

[fluidistic]

OK, my apologies, I confused your use of the fluid stream function $\psi$ with the potential function $\phi$.
So as I now understand it, you're interested in the 2D vorticity vector $\vec{\omega}\equiv\omega\left(r,\theta,A\right)\hat{e}_{z}\equiv\nabla\times\vec{v}=-\nabla^{2}\psi\,\hat{e}_{z}\,$, where $\omega\left(r,\theta,A\right)\equiv-\frac{A\,\cos2\theta}{r^{2}}$. Since it's easy to verify that $\nabla^{2}\omega\left(r,\theta,A\right)=0$, your stream function $\psi(r, \theta)=\left( \frac{c_1}{r^2}+c_2r^2-\frac{A}{4}\right)\cos2\theta$ actually satisfies the homogeneous 4th-order biharmonic equation $\nabla^{4}\psi=0$. This puts you squarely into the realm of Stokes flow, which requires two boundary conditions on each surface (http://brennen.caltech.edu/FLUIDBOOK/basicfluiddynamics/Stokesflow/equationsofstokesflow.pdf):
View attachment 350977
Because of this, satisfying your boundary conditions is possible, but it unfortunately requires a substantially more complicated stream function than the one shown above (http://brennen.caltech.edu/fluidbook/basicfluiddynamics/Stokesflow/planarflow.pdf):
View attachment 350978
Another useful reference to these equations is https://personalpages.manchester.ac...il/Lectures/Fluids/Material_2016/Chapter8.pdf:
View attachment 350979
Hope you find this helpful.

I see... I will study your reply. Numerically though, I certainly do not get a vanishing $\vec v$ on the boundaries. I am skeptical this problem falls under this category.

It's also a bit crazy to see the difference in complications if the $\cos (2\theta)$ is replaced by $\sin(2\theta)$ which also makes sense physically but corresponds to applying a different boundary conditions to a linked PDE. There, $\vec v$ is also non vanishing on the boundaries.

 

[fluidistic]

@renormalize OK, that was very informative! Yes, "my" $\psi$ apparently does satisfy the biharmonic equation. That's really nice, to learn that the current density (my $\vec v$) behaves like a fluid under Stokes conditions.

Ok, and thanks a lot for the general form of $\psi(r, \theta)$. I couldn't have asked for any better.
I do, however, face big concerns. When I apply the $\psi(r, \theta=0)$ boundary condition, I get that a bunch of constants have to vanish, and similarly for the condition at $\theta = \pi/2$. This leaves me with $\psi(r, \theta)=B_0r^2[1-\cos(2\theta)-r^2\sin(2\theta)]+Cr^3[\cos(\theta)-\cos(3\theta)]$. I verified that it indeed satisfies the $\theta$ boundary conditions.

Now the real killer, the boundary conditions at $r=r_i$ and $r=r_o$. When $r=r_i$, nothing must depend on theta as psi must vanish. Boom. The constant $B_0$ and $C$ have to vanish, making psi vanishing itself everywhere, as usual, i.e. same as I got ever since I opened this thread.

Do I miss something?

 

[renormalize]

Do I miss something?

Can you please clarify for me: do you seek solutions of the biharmonic eq. $\nabla^{4}\psi=0$ that satisfy $\psi=0$ on the boundaries (Dirichlet BC), or $\nabla_{\Vert}\psi=\nabla_{\bot}\psi=0$ on the boundaries (Stokes flow BCs), or possibly some other boundary condition?

 

[fluidistic]

Can you please clarify for me: do you seek solutions of the biharmonic eq. $\nabla^{4}\psi=0$ that satisfy $\psi=0$ on the boundaries (Dirichlet BC), or $\nabla_{\Vert}\psi=\nabla_{\bot}\psi=0$ on the boundaries (Stokes flow BCs), or possibly some other boundary condition?

I haven't thought thoroughly about all the b.c. I need to apply. I started by applying the one I know to hold, namely $\psi=0$ on the boundaries. This already killed the whole function $\psi$.

 

[renormalize]

I started by applying the one I know to hold, namely $\psi=0$ on the boundaries.

How do you know that BC holds? What do you actually require of the components of $\vec{v}$ on the boundary? Or, indeed, of the vorticity $\omega\equiv -\nabla^{2}\psi$ on the boundary? Maybe you should solve $0=\nabla^{4}\psi=-\nabla^{2}\omega$ with the condition that $\omega=0$ on the boundary. Then, just as with the electrostatic potential inside a conducting shell, you get the solution $\omega=-\nabla^{2}\psi=\text{const.}$ everywhere inside the boundary.

 

[fluidistic]

How do you know that BC holds? What do you actually require of the components of $\vec{v}$ on the boundary? Or, indeed, of the vorticity $\omega\equiv -\nabla^{2}\psi$ on the boundary? Maybe you should solve $0=\nabla^{4}\psi=-\nabla^{2}\omega$ with the condition that $\omega=0$ on the boundary. Then, just as with the electrostatic potential inside a conducting shell, you get the solution $\omega=-\nabla^{2}\psi=\text{const.}$ everywhere inside the boundary.

As far as I understand, $\psi=0$ on the boundaries is the only way to ensure the normal component of $\vec v$ vanish there, which is my physical requirement. For example, if I take the original expression I had above, the one involving a $\cos (2\theta)$, $\psi$ does not satisfy this boundary condition. I retrieved the corresponding $\vec v$ to make sure it is incorrect, and it is indeed incorrect in that $\vec v$ enters/leave the material. In my case, $\vec v$ is actually the current density, and one of the conditions it has to satisfy is that no current enters/leave the material, therefore $\vec v$ has to be parallel to the sides/boundaries, there can be no perpendicular component.

I haven't been able to find a single expression for $\psi$ that satisfy this condition as well as the PDE in the OP, so far at least.

 

[fluidistic]

Now that I think about it, having a curl of v proportional to cos(2 theta) is terrible because it means that this curl does not vanish at theta equals 0 and pi/2, meaning v cannot be parallel to those sides, meaning v has to have a non zero normal component if v is non zero there.

This explains my troubles, I guess. What do you think?
Then I will need to explain how I got a cos(2 theta) term in the first place...


Edit: nevermind this comment. It's perfectly possible to have a non vanishing curl along a side.and have the field to be parallel to that side. So, ok, I still have hopes to find a solution that satisfies the vanishing b.c.s. for psi.

 

[renormalize]

Edit: nevermind this comment. It's perfectly possible to have a non vanishing curl along a side.and have the field to be parallel to that side. So, ok, I still have hopes to find a solution that satisfies the vanishing b.c.s. for psi.

I think you need to precisely pin down your boundary conditions based on the physics underlying your problem. Perhaps I can help you do that by enumerating some of the possibilities.
Start from the basic definitions of the 2D velocity-field (or "current-field" in your parlance) $\vec{v}$ and the vorticity-field $\omega$ in terms of the stream function $\psi\,$:$$
\vec{v}\equiv\nabla\times\psi\,\hat{e}_{z}\,,\;\vec{\omega}\equiv\nabla\times\vec{v}\equiv\omega\,\hat{e}_{z}=-\nabla^{2}\psi\,\hat{e}_{z}\tag{1a,1b}$$Then for incompressible, non-viscous flow, the 2D Navier-Stokes equation reduces to a single 4th-order differential equation for the stream function:$$\nabla^{4}\psi=-\nabla^{2}\omega=0\tag{2}$$(For the argument I will give below, it's important that eq.(2) should hold everywhere in the domain $A$ of your problem; i.e., there can be no point sources or sinks of vorticity anywhere within $A$ or on its boundary $\partial A$. Please let me know if this assumption is incorrect!) Now into the identity:$$\left(\nabla^{2}\psi\right)^{2}\equiv\nabla\cdot\left[\nabla\psi\left(\nabla^{2}\psi\right)-\psi\nabla\left(\nabla^{2}\psi\right)\right]+\psi\nabla^{4}\psi\tag{3}$$I substitute eqs.(1,2) to find:$$\omega^{2}=\nabla\cdot\left(\psi\nabla\omega-\omega\nabla\psi\right)\tag{4}$$and then integrate over $A$ and apply the divergence theorem to ultimately obtain:$$\intop_{A}\omega^{2}\,d^{2}x=\oint_{\partial A}\left(\psi\nabla\omega-\omega\nabla\psi\right)\cdot d\vec{s}=\oint_{\partial A}\left(\psi\frac{\partial\omega}{\partial n}-\omega\frac{\partial\psi}{\partial n}\right)ds\tag{5}$$This is a useful result because it relates the square-of-vorticity averaged over the domain to the fields and their normal-derivatives evaluated along the boundary of the domain.
On this boundary, denoting the velocity components parallel and perpendicular to the boundary wall by $v_{\Vert}\,,v_{\bot}\,$, standard fluid dynamics demands that $v_{\bot}=0$ (no fluid penetrates into the walls) and $v_{\Vert}=0$ relative to the wall (no fluid slips at the walls). In terms of the stream function, these translate to $\psi=\text{const.}$ and $\partial\psi /\partial n=\text{wall tangential-velocity}$ on the boundary. So in particular, no-penetration and stationary walls compel setting $\psi=\partial\psi /\partial n=0$ in eq.(5), yielding $\intop_{A}\omega^{2}d^{2}x=0$. Because the integrand $\omega^2$ is positive definite, this forces $\omega=0$ everywhere; i.e., the flow within stationary walls is irrotational and thus $\vec{v}=\nabla\phi$. Indeed, due to fluid incompressibility ($\nabla\cdot\vec{v}=0$), $\vec{v}\equiv 0$ everywhere unless balanced sources and sinks of fluid are present to allow potential flow. (As a consequence, simple CFD test problems intended to demonstrate vorticity employ moving surfaces, like the top of a box sliding horizontally, or the rotation of the inner cylinder of an annular tank, to induce vorticity within the contained fluid.)
Now regarding your specific problem: the most basic condition on your fluid is that it must be contained: $v_{\bot}=0$, so on the walls you set $\psi=\text{const.}=0$ without loss of generality. It then follows by eq.(5) that non-vanishing vorticity $\omega$ is possible only if the product $\omega\,\partial\psi /\partial n=\left(\hat{e}_z\cdot\nabla\times\vec{v}\right)\,v_{\Vert}$ is nonzero over some finite length along the boundary. In other words, to obtain a nontrivial solution, either your problem includes a moving wall or you must drop the "no slip" condition (i.e., $v_{\Vert}\neq 0$) somewhere along your stationary boundary and the vorticity $\omega$ must also be nonzero over at least a finite part of that same boundary region.
Are these various boundary restrictions compatible with the physics of your problem?

 

[fluidistic]

I think you need to precisely pin down your boundary conditions based on the physics underlying your problem. Perhaps I can help you do that by enumerating some of the possibilities.
Start from the basic definitions of the 2D velocity-field (or "current-field" in your parlance) $\vec{v}$ and the vorticity-field $\omega$ in terms of the stream function $\psi\,$:$$
\vec{v}\equiv\nabla\times\psi\,\hat{e}_{z}\,,\;\vec{\omega}\equiv\nabla\times\vec{v}\equiv\omega\,\hat{e}_{z}=-\nabla^{2}\psi\,\hat{e}_{z}\tag{1a,1b}$$Then for incompressible, non-viscous flow, the 2D Navier-Stokes equation reduces to a single 4th-order differential equation for the stream function:$$\nabla^{4}\psi=-\nabla^{2}\omega=0\tag{2}$$(For the argument I will give below, it's important that eq.(2) should hold everywhere in the domain $A$ of your problem; i.e., there can be no point sources or sinks of vorticity anywhere within $A$ or on its boundary $\partial A$. Please let me know if this assumption is incorrect!) Now into the identity:$$\left(\nabla^{2}\psi\right)^{2}\equiv\nabla\cdot\left[\nabla\psi\left(\nabla^{2}\psi\right)-\psi\nabla\left(\nabla^{2}\psi\right)\right]+\psi\nabla^{4}\psi\tag{3}$$I substitute eqs.(1,2) to find:$$\omega^{2}=\nabla\cdot\left(\psi\nabla\omega-\omega\nabla\psi\right)\tag{4}$$and then integrate over $A$ and apply the divergence theorem to ultimately obtain:$$\intop_{A}\omega^{2}\,d^{2}x=\oint_{\partial A}\left(\psi\nabla\omega-\omega\nabla\psi\right)\cdot d\vec{s}=\oint_{\partial A}\left(\psi\frac{\partial\omega}{\partial n}-\omega\frac{\partial\psi}{\partial n}\right)ds\tag{5}$$This is a useful result because it relates the square-of-vorticity averaged over the domain to the fields and their normal-derivatives evaluated along the boundary of the domain.
On this boundary, denoting the velocity components parallel and perpendicular to the boundary wall by $v_{\Vert}\,,v_{\bot}\,$, standard fluid dynamics demands that $v_{\bot}=0$ (no fluid penetrates into the walls) and $v_{\Vert}=0$ relative to the wall (no fluid slips at the walls). In terms of the stream function, this translates to $\psi=\text{const.}$ and $\partial\psi /\partial n=\text{wall-tangential-velocity}$ on the boundary. So in particular, no-penetration and stationary walls compel setting $\psi=\partial\psi /\partial n=0$ in eq.(5), yielding $\intop_{A}\omega^{2}d^{2}x=0$. Because the integrand $\omega^2$ is positive definite, this forces $\omega=0$ everywhere; i.e., the flow within stationary walls is irrotational and thus $\vec{v}=\nabla\phi$. Indeed, due to fluid incompressibility ($\nabla\cdot\vec{v}=0$), $\vec{v}= 0$ everywhere if there are no balanced sources and sinks of fluid to allow potential flow. (As a consequence, CFD test problems intended to demonstrate vorticity employ moving surfaces, like the top of a box sliding horizontally, or the rotation of the inner cylinder of an annular tank, to induce vorticity within the contained fluid.)
Now regarding your specific problem: the most basic condition on your fluid is that it must be contained ($v_{\bot}=0$), so on the walls you set $\psi=\text{const.}=0$ without loss of generality. To obtain non-vanishing vorticity $\omega$, per the line-integral on the right of eq.(5), this in turn requires that the product $\omega\,\partial\psi /\partial n=\left|\nabla\times\vec{v}\right|\,v_{\Vert}$ must be nonzero over some finite length along the boundary. So to obtain a nontrivial solution, either your problem includes a moving wall or you must drop the "no slip" condition somewhere along your stationary boundary. And overlapping that same boundary region, the vorticity must also be nonzero for a finite length. Are these various boundary restrictions compatible with the physics of your problem?

First of all, thanks again for your detailed answer, I will need some time to follow you,entirely. But I,can already respond to some of your inqueries. Yes, in my,case.there are vortices, point vortices inside the domain, so I do not know how this changes things regarding your conclusions. There is 1 vortex in each half of.the quarter of annulus. And since.the curl of J (v using your notation), they,have opposite vorticity. So I guess that psi is also constant and vanishing at.theta equals pi/4, i.e.at the middle of the quarter of an annulus but this shouldn't be a boundary,condition to enforce, rather something that one should see if one has the solution.


Unless I can effectively split.the problem in 2 and impose.psi equal to 0 in the eighth of an annulus. I already tried but I wasn't able.to get a solution either. I am sure the solution exists though, it's just.that I,am not able.to get it.

Physically v must.be tangential to all boundaries, and.the curl has to satisfy what I wrote, those are.the only requirement for v.

 

[renormalize]

There is 1 vortex in each half of.the quarter of annulus. And since.the curl of J (v using your notation), they,have opposite vorticity.

This I don't understand. In your OP you wrote the equation $\nabla^{2}\psi=\frac{A\,\cos2\theta}{r^{2}}$. If your problem has point vortices, shouldn't the right side contain some delta functions at specific locations? Or do you mean that you're actually trying to solve something like $\nabla^{2}\omega=-\nabla^{4}\psi=q\left(\delta^{3}\left(x_{1}\right)-\delta\left(x_{2}\right)\right)$ in the quarter-anulus?

 

[pasmith]

@renormalize OK, that was very informative! Yes, "my" $\psi$ apparently does satisfy the biharmonic equation. That's really nice, to learn that the current density (my $\vec v$) behaves like a fluid under Stokes conditions.

You say that $v$ (in your notation) is a current density. @renormalize is trying to treat it as a fluid velocity, which may not be appropriate. Please check the physics.

 

[fluidistic]

This I don't understand. In your OP you wrote the equation $\nabla^{2}\psi=\frac{A\,\cos2\theta}{r^{2}}$. If your problem has point vortices, shouldn't the right side contain some delta functions at specific locations? Or do you mean that you're actually trying to solve something like $\nabla^{2}\omega=-\nabla^{4}\psi=q\left(\delta^{3}\left(x_{1}\right)-\delta\left(x_{2}\right)\right)$ in the quarter-anulus?

No, no Dirac deltas. If I have the time tonight I might show you the numerical solution. Or take the analytical solution I posted for.the case involving the sine case. In that case there is only one vortex in the middle of the quarter of the annulus. And there is no delta function in the curl of v (or J to be more specific). In that case the streamlines are concentric around the central vortex. This can be unserstood because.sine of 2 theta does not change sign in the material, and J or v is contained, which means the flux can rotate around a single point in the xy plane.


In my,case.with cos 2 theta, this expression changes sign at.the theta equals pi/4, meaning the spinning reverse when you vary theta and cross.that line. Numerically I see a vortex in each one of those halves, so each half is a bit similar to the case of sin 2 theta, but possibly a bit different, I,am not sure yet.

 

[fluidistic]

You say that $v$ (in your notation) is a current density. @renormalize is trying to treat it as a fluid velocity, which may not be appropriate. Please check the physics.

Indeed. I,already specified all the requirements my ''J'' or v has to satisfy, as well as the boundary conditions. I also verified that psi, its streamlines function, satisfies the biharmonic equation, although I am unable to get a psi that satisfies vanishing b.c.s on all,the boundaries, but this condition is of,utmost importance.

 

[renormalize]

No, no Dirac deltas. If I have the time tonight I might show you the numerical solution. Or take the analytical solution I posted for.the case involving the sine case. In that case there is only one vortex in the middle of the quarter of the annulus. And there is no delta function in the curl of v (or J to be more specific). In that case the streamlines are concentric around the central vortex. This can be unserstood because.sine of 2 theta does not change sign in the material, and J or v is contained, which means the flux can rotate around a single point in the xy plane.


In my,case.with cos 2 theta, this expression changes sign at.the theta equals pi/4, meaning the spinning reverse when you vary theta and cross.that line. Numerically I see a vortex in each one of those halves, so each half is a bit similar to the case of sin 2 theta, but possibly a bit different, I,am not sure yet.

OK, I'm even more confused. So can you please:

• post a simple diagram of your annulus that shows the boundaries and the location(s) of the point vortex/vortices.
• clearly write down the differential equation that you are trying to solve analytically, including the manifest form of the point sources (which you say are not Dirac deltas).

Thanks!

 

[fluidistic]

OK, I'm even more confused. So can you please:

• post a simple diagram of your annulus that shows the boundaries and the location(s) of the point vortex/vortices.
• clearly write down the differential equation that you are trying to solve analytically, including the manifest form of the point sources (which you say are not Dirac deltas).

Thanks!

Sure, I can possibly do that within 7 hours from now. I just tested and I am able to display the streamlines from the analytical expression for the problem involving the sine term rather than the cosine one. They match my numerical results at first glance.

I can show you my numerical results for.the case involving the cosine term, this should give us an indication of how the analytical solution should look like.

But the math is already in this thread. The analytical psi I posted above for the sine problem satisfies its Poisson equation and boundary conditions, and it does contain a vortex. If you plot its constant streamlines, they are circling around a single point and covering the whole quarter of an annulus for the ones of lower intensity.

 

[fluidistic]

Alrighty, here I go (on a computer so I shouldn't be as sloppy as in some of my previous replies).
I will not go into details in how I derived the PDEs, but I have to solve for $\vec J$ in several coupled PDEs. There is, for instance, $\nabla \cdot \vec J =0$ as well as $\nabla \times \vec J \propto \cos(2\theta)/r^2$. The region in which I must solve for $\vec J$ is a quarter of an annulus, $\theta$ ranging from $0$ to $\pi/2$ and $r$ ranging from $r_i$ to $r_o$. The material being isolated, cannot have any entering nor leaving current density at any point. This means that no normal component of $\vec J$ can exist.

If $\vec J(r, \theta)=J_r(r, \theta) \hat r + J_\theta(r, \theta) \hat \theta$, then it means that $J_r(r=r_i)=J_r(r=r_o)=0$ and $J_\theta(\theta=0)=J_\theta(\theta=\pi/2)=0$. Those are the boundary conditions. I wrote down the divergence and the curl in polar coordinates, I did 3 pages of algebra to manipulate the coupled PDE's to finally have a separable form for $J_\theta$ and an uglier form for $J_r$. Unfortunately I was not able to find any solution that would make the boundary conditions to hold. I gave up and tried the streamlines approach. Constant streamlines share the same direction than the current density.

Using this approach, the vector problem becomes a "simple" scalar PDE, where the condition of the vanishing divergence of $\vec J$ is automatically enforced. The problem of the curl is simplified to a Poisson problem, just like in electrostatics. The problem becomes $\nabla ^2 \psi(r, \theta)\propto \cos(2\theta)/r^2$. The boundary conditions on $\psi$ are such that it must vanish on all boundaries (in fact it must be constant, but this constant is arbitrary and taken as 0 here for simplicity). Therefore I should be able to use any method taught to solve this PDE, namely with a Green function (the problem reduces, I suppose, then, to solve a really ugly integral), or separation of variables solving first the homogeneous problem and then the particular one. But, as evidenced above, I could not, for the life of me, find a coherent answer. I suspect I made a mistake and not all the terms cancel out to 0 in the PDF screenshots posted above, but I haven't found out where my mistake is.
Ok, and we have $J_r =-\frac{1}{r}\frac{\partial \psi}{\partial \theta}$ and $J_\theta = \frac{\partial \psi}{\partial r}$, meaning that solving for $\psi(r, \theta)$ yields $\vec J$ which is what we're after.

Before I display the numerical results, I wanted to mention that a similar problem is much simpler to solve. If you replace the $\cos(2\theta)$ by $\sin(2\theta)$, then $\psi(r, \theta)=\psi_0 \left[ 1-\left(\frac{1}{r_i^2+r_o^2} \right) \left( r^2+\frac{r_i^2r_o^2}{r^2} \right) \right]\sin(2\theta)$, it satisfies the same boundary condition as in my problem. In fact, it turns out that this change of problem still corresponds to another physical problem, so it's not just a mathematical curiosity. I can then compare my numerical result to this exact solution. There is a slight difference in that my numerical problem treats a slightly more complicated problem and the analytical solution is a solution to a slightly simplified problem.

Here's a picture of $\vec J$ for this similar problem. Coloring is based on the magnitude of $\vec J$, the arrows only indicate the direction of $\vec J$. Not sure if you can see it, but $\vec J$ has a single vortex, the current density swirls around it. It's not clear to me whether $\vec J$ vanishes there, it seems to occur right at the center to me, so we would have to retrieve $\vec J$ from the analytical formula and check if $\theta =\pi/4$, $r=(r_o+r_i)/2$ makes $\vec J=\vec 0$.

And here's a zoom on the center, to see the vortex more closely.

And here's a plot of several constant streamlines, taken from the analytical solution.

To me, visually at least, they seem to be in good agreement with the current density computed numerically.

Finally, here is my numerical result of the original problem, the one I am still looking for its analytical solution:

I am not sure if you can see it, but I see at least 2 vortices. At the bottom, a vortex with a clockwise swirl. At the top, counter clockwise swirl, which is inline with the curl changing sign about $\pi/4$, where it vanishes (and $\vec J$ seems to be following an inward straight line there, which makes sense). To me, it looks like the swirls do not occur exactly at $\theta = \pi/8$, they seem to be closer to the bottom and top edges than to the center. The magnitude of $\vec J$ seems to decrease as $r$ increases. Not sure if we can get much more from this.

 

[renormalize]

And here's a plot of several constant streamlines, taken from the analytical solution.

Thanks for posting all the detailed info. After seeing the streamlines from your analytical solution, a light bulb lit up! Of course I haven't seen your derivations, but I simply cannot reconcile your desired equation $\nabla ^2 \psi(r, \theta)\propto \cos(2\theta)/r^2$ with your required boundary conditions. So let me offer a different analysis of your problem.
Consider the general series-solution of $\nabla^{4}\psi=0$ exhibited in the references shown previously in my post #7. In terms of the polar angle $\theta$, that solution is a Fourier series involving $\text{cos}(m\theta),\text{sin}(m\theta)$, where $m=0,1,2,\ldots\,$. But as you've already realized, it's not possible to match your BCs using cosines, so I zero-out all the terms with $\text{cos}(m\theta)$, and drop as well those with $\text{sin}(m\theta)$ whenever $m$ is an odd-integer. That leaves me with the specific series solution:$$\psi\left(r,\theta\right)=\sum_{n=1}^{\infty}\left(a_{n}r^{2n}+b_{n}r^{-2n}+c_{n}r^{2n+2}+d_{n}r^{2-2n}\right)\sin\left(2n\theta\right)\equiv\sum_{n=1}^{\infty}\psi_{n}\left(r,\theta\right)\tag{1}$$By construction, each term $\psi_{n}$ in this series is a solution of $\nabla^{4}\psi_n=0$ and also manifestly satisfies your condition $\psi_n\left(r,0\right)=\psi_n\left(r,\frac{\pi}{2}\right)=0$. Moreover, 2 of the 4 arbitrary constants in each term can be chosen so that $\psi_n\left(r_{i},\theta\right)=\psi_n\left(r_{o},\theta\right)=0$, thereby resulting in:
\begin{align}
\psi_{n}\left(r,\theta\right) & =\left\{ \left[\frac{\left(c_{n}r_{i}^{4n}r_{o}^{4n}\left(r_{i}^{2}-r_{o}^{2}\right)+d_{n}\left(r_{o}^{4n}r_{i}^{2}-r_{i}^{4n}r_{o}^{2}\right)\right)}{r_{i}^{4n}-r_{o}^{4n}}\right]r^{-2n}\right.\nonumber \\
&\quad \left.-\left[\frac{\left(c_{n}\left(r_{i}^{4n+2}-r_{o}^{4n+2}\right)+d_n\left(r_{i}^{2}-r_{o}^{2}\right)\right)}{r_{i}^{4n}-r_{o}^{4n}}\right]r^{2n}+c_{n}r^{2n+2}+d_{n}r^{2-2n}\right\} \sin\left(2n\theta\right) \nonumber\tag{2}
\end{align}
The set of all these $\psi_n$ are the "normal modes" for your problem because they can be superposed with the weightings $c_n,d_n$ to form the most general stream function $\psi$ that vanishes everywhere along the boundary of your annulus.
I turn next to the vorticity $\omega\left(r,\theta\right)\equiv\sum_{n=1}^{\infty}\omega_{n}\left(r,\theta\right)$, where:$$\omega_{n}\left(r,\theta\right)\equiv-\nabla^{2}\psi_{n}\left(r,\theta\right)=-4\left[\left(2n+1\right)c_{n}r^{2n}-\left(2n-1\right)d_{n}r^{-2n}\right]\sin\left(2n\theta\right)\tag{3}$$From this it's clear that the constants $c_n,d_n$ can be adjusted to set specified boundary values for $\omega_n$ on the inner and outer cylinder walls. Alternatively, these constants can be tied to the tangential velocities $v_{\Vert}=\frac{\partial \psi}{\partial r}$ of the fluid on those same walls. Picking these values requires a knowledge of the physics you're trying to model (which I lack), so I leave selecting $c_n,d_n$ to you.
That said, it's easy to use Mathematica to evaluate and plot the stream-function normal-modes $\psi_n$ from eq.(2) if I choose arbitrary values for $c_n,d_n$. For example, the values $r_{i}=0.3,r_{o}=0.4,c_{n}=1,d_{n}=1,n=1,2,3,4$ produce:

Clearly $n$ determines the number of distinct vortexes and, as you observed, each vortex is adjacent to another that rotates in the opposite direction. The similarity of your graphs to my $n=1,2$ plots above suggests to me that your CFD code is computing something akin to $\nabla ^2 \psi_n(r, \theta)\propto \sin(2n\theta)$ rather than anything involving a cosine function.
And for completeness, I also show the vorticity normal-modes $\omega_n$ of eq.(3) for the same conditions:

Please let me know your thoughts.

 

[fluidistic]

Thanks for posting all the detailed info. After seeing the streamlines from your analytical solution, a light bulb lit up! Of course I haven't seen your derivations, but I simply cannot reconcile your desired equation $\nabla ^2 \psi(r, \theta)\propto \cos(2\theta)/r^2$ with your required boundary conditions. So let me offer a different analysis of your problem.
Consider the general series-solution of $\nabla^{4}\psi=0$ exhibited in the references shown previously in my post #7. In terms of the polar angle $\theta$, that solution is a Fourier series involving $\text{cos}(m\theta),\text{sin}(m\theta)$, where $m=0,1,2,\ldots\,$. But as you've already realized, it's not possible to match your BCs using cosines, so I zero-out all the terms with $\text{cos}(m\theta)$, and drop as well those with $\text{sin}(m\theta)$ whenever $m$ is an odd-integer. That leaves me with the specific series solution:$$\psi\left(r,\theta\right)=\sum_{n=1}^{\infty}\left(a_{n}r^{2n}+b_{n}r^{-2n}+c_{n}r^{2n+2}+d_{n}r^{2-2n}\right)\sin\left(2n\theta\right)\equiv\sum_{n=1}^{\infty}\psi_{n}\left(r,\theta\right)\tag{1}$$By construction, each term $\psi_{n}$ in this series is a solution of $\nabla^{4}\psi_n=0$ and also manifestly satisfies your condition $\psi_n\left(r,0\right)=\psi_n\left(r,\frac{\pi}{2}\right)=0$. Moreover, 2 of the 4 arbitrary constants in each term can be chosen so that $\psi_n\left(r_{i},\theta\right)=\psi_n\left(r_{o},\theta\right)=0$, thereby resulting in:
\begin{align}
\psi_{n}\left(r,\theta\right) & =\left\{ \left[\frac{\left(c_{n}r_{i}^{4n}r_{o}^{4n}\left(r_{i}^{2}-r_{o}^{2}\right)+d_{n}\left(r_{o}^{4n}r_{i}^{2}-r_{i}^{4n}r_{o}^{2}\right)\right)}{r_{i}^{4n}-r_{o}^{4n}}\right]r^{-2n}\right.\nonumber \\
&\quad \left.-\left[\frac{\left(c_{n}\left(r_{i}^{4n+2}-r_{o}^{4n+2}\right)+d_n\left(r_{i}^{2}-r_{o}^{2}\right)\right)}{r_{i}^{4n}-r_{o}^{4n}}\right]r^{2n}+c_{n}r^{2n+2}+d_{n}r^{2-2n}\right\} \sin\left(2n\theta\right) \nonumber\tag{2}
\end{align}
The set of all these $\psi_n$ are the "normal modes" for your problem because they can be superposed with the weightings $c_n,d_n$ to form the most general stream function $\psi$ that vanishes everywhere along the boundary of your annulus.
I turn next to the vorticity $\omega\left(r,\theta\right)\equiv\sum_{n=1}^{\infty}\omega_{n}\left(r,\theta\right)$, where:$$\omega_{n}\left(r,\theta\right)\equiv-\nabla^{2}\psi_{n}\left(r,\theta\right)=-4\left[\left(2n+1\right)c_{n}r^{2n}-\left(2n-1\right)d_{n}r^{-2n}\right]\sin\left(2n\theta\right)\tag{3}$$From this it's clear that the constants $c_n,d_n$ can be adjusted to set specified boundary values for $\omega_n$ on the inner and outer cylinder walls. Alternatively, these constants can be tied to the tangential velocities $v_{\Vert}=\frac{\partial \psi}{\partial r}$ of the fluid on those same walls. Picking these values requires a knowledge of the physics you're trying to model (which I lack), so I leave selecting $c_n,d_n$ to you.
That said, it's easy to use Mathematica to evaluate and plot the stream-function normal-modes $\psi_n$ from eq.(2) if I choose arbitrary values for $c_n,d_n$. For example, the values $r_{i}=0.3,r_{o}=0.4,c_{n}=1,d_{n}=1,n=1,2,3,4$ produce:
View attachment 351135
Clearly $n$ determines the number of distinct vortexes and, as you observed, each vortex is adjacent to another that rotates in the opposite direction. The similarity of your graphs to my $n=1,2$ plots above suggests to me that your CFD code is computing something akin to $\nabla ^2 \psi_n(r, \theta)\propto \sin(2n\theta)$ rather than anything involving a cosine function.
And for completeness, I also show the vorticity normal-modes $\omega_n$ of eq.(3) for the same conditions:
View attachment 351136
Please let me know your thoughts.

Thanks a million, once again for your interest, time and insights. I don't have the time right now to give a detailed reply nor to go through math, but I,believe there is a solution we miss. I am thinking maybe the angular part is cos(2 theta) plus a quadratic function of theta, which would cancel out the cos term at theta equal to 0 and pi/2 as well as satisfying the Poisson equation. The quqdratic part arises when solving the angular part when n equals 0. What do you think?

I have good reasons to believe a curl of J should be proportional to cos 2 theta. This holds if the temperature profile in the annulus is linear in theta, which corresponds to applying a Delta T at theta equals 0 and pi/2.
In the case of.the curl proportional to sin(2 theta), this arises when the Delta T is applied between ri and ro, i.e. T is purely radial, equal to a constant plus ln(r).

 

[fluidistic]

I have a few more minutes. I got the $\cos(2\theta)/r^2$ curl expression from the fact that $\nabla \times \vec J \propto \frac{\partial ^2T}{\partial y \partial x}$ where $T(x,y)$ is the temperature distribution in the material. If the sides $\theta = 0$ and $\theta = \pi/2$ are kept at fixed temperature differing by $\Delta T$, the temperature distribution is linear in $\theta$: $T(r, \theta)=T_0 + \frac{2}{\pi}\Delta T \theta$.
In order to compute $\frac{\partial ^2T}{\partial y \partial x}$, I replaced $\theta$ by $\arctan (y/x)$, performed the derivatives and converted back to polar, yielding the $\cos(2\theta)/r^2$ term.
I redid those computations about 10 times, I haven't found any mistake. But maybe I have done a mistake... I wish I did. I wish I would find something like $\sin(4\theta)/r^2$. Perhaps I goof when I deal with the $\arctan$ function? But I don't think so...

 

[renormalize]

$\nabla \times \vec J \propto \frac{\partial ^2T}{\partial y \partial x}$ where $T(x,y)$ is the temperature distribution in the material.

This expression declares that the vector $\nabla \times \vec J$ is proportional to an off-diagonal component $\frac{\partial ^2T}{\partial y \partial x}$ of the symmetric rank-2 tensor $\nabla(\nabla T)$. How do I make sense of this as a conventional tensor relationship defined so as to be covariant under rotations? In other words, can you rewrite it in standard tensor-component form, something like $\epsilon_{ijk}J_{j,i}\propto\ldots\,$?

 

[fluidistic]

This expression declares that the vector $\nabla \times \vec J$ is proportional to an off-diagonal component $\frac{\partial ^2T}{\partial y \partial x}$ of the symmetric rank-2 tensor $\nabla(\nabla T)$. How do I make sense of this as a conventional tensor relationship defined so as to be covariant under rotations? In other words, can you rewrite it in standard tensor-component form, something like $\epsilon_{ijk}J_{j,i}\propto\ldots\,$?

By adding a $\hat k$ that I forgot: $\vec J \propto \frac{\partial ^2T}{\partial y \partial x} \hat k$.

 

[renormalize]

By adding a $\hat k$ that I forgot: $\vec J \propto \frac{\partial ^2T}{\partial y \partial x} \hat k$.

Did you mean to write $\nabla\times\vec J \propto \frac{\partial ^2T}{\partial y \partial x} \hat k\,$?
If so, from the basic equations we've been using in this thread, your proportionality states:$$\vec{J}\equiv\nabla\times\psi\,\vec{k}\Rightarrow\nabla\times\vec{J}=\nabla\times\nabla\times\psi\,\vec{k}=-\nabla^{2}\psi\,\vec{k}\propto\frac{\partial^{2}T}{\partial x\partial y}\,\vec{k}$$or, dropping the unit vectors, $-\nabla^{2}\psi\propto\frac{\partial^{2}T}{\partial x\partial y}\,$. My problem with this relation is that the quantity on the left is a manifest scalar under 2D rotations, whereas the right-side is the off-diagonal element of the 2x2 tensor $\nabla(\nabla T)$ and therefore is not a scalar under 2D rotations. So I'm asking you: how can the scalar field $\psi$ be sourced by the non-scalar quantity $\frac{\partial^{2}T}{\partial x\partial y}$? (I'm asking this because I suspect it's the source of your original difficulty.)

 

[fluidistic]

Did you mean to write $\nabla\times\vec J \propto \frac{\partial ^2T}{\partial y \partial x} \hat k\,$?
If so, from the basic equations we've been using in this thread, your proportionality states:$$\vec{J}\equiv\nabla\times\psi\,\vec{k}\Rightarrow\nabla\times\vec{J}=\nabla\times\nabla\times\psi\,\vec{k}=-\nabla^{2}\psi\,\vec{k}\propto\frac{\partial^{2}T}{\partial x\partial y}\,\vec{k}$$or, dropping the unit vectors, $-\nabla^{2}\psi\propto\frac{\partial^{2}T}{\partial x\partial y}\,$. My problem with this relation is that the quantity on the left is a manifest scalar under 2D rotations, whereas the right-side is the off-diagonal element of the 2x2 tensor $\nabla(\nabla T)$ and therefore is not a scalar under 2D rotations. So I'm asking you: how can the scalar field $\psi$ be sourced by the non-scalar quantity $\frac{\partial^{2}T}{\partial x\partial y}$? (I'm asking this because I suspect it's the source of your original difficulty.)

Thanks again, this is a bit over my head. I hope the following gives you an answer.
We start with $\vec J = - \sigma \nabla V-\sigma S \nabla T$, where $V=V(x,y)$ and $T=T(x,y)$ are scalar fields, $\sigma$ is a 2x2 matrix (or tensor of order 1 if I am not mistaken), but it has the same components on its diagonal and the off-diagonal elements are zero, so I reduce it to a scalar. $S$ is like sigma, but it has different diagonal components, making the physical property "S" anisotropic in the x-y plane. Its off-diagonal components are zeroes.
When I compute the curl of $\vec J$, I get:
$\nabla \times \vec J = -\sigma \left(S_{xx}\frac{\partial ^2 T}{\partial x \partial y} -S_{yy}\frac{\partial ^2 T}{\partial y \partial x} \right) \hat k$. Since the crossed derivatives are equal to one another, I can factor them out. I am left with the curl of J being proportional to those crossed derivatives.
Is there anything wrong with this?

 

[renormalize]

Is there anything wrong with this?

Thanks for this clarification! I think I can now resolve your problem. As I learned from a modicum of googling, your current:$$\vec{J}\equiv\nabla\times\psi\,\hat{k}=\sigma\left(-\nabla V\right)+\sigma\overleftrightarrow{S}\left(-\nabla T\right)\tag{1}$$is the sum of the usual ohmic electric current, induced by a voltage gradient in a medium characterized by a scalar conductivity $\sigma$, with a thermoelectric current induced by a thermal gradient, as characterized by a symmetric tensor $\overleftrightarrow{S}$ of Seebeck coefficients. With (1) now properly expressed as a 3-vector equation, I can use standard coordinate-transformation rules to convert it to polar coordinates.
Since your problem is independent of the $z$-coordinate, I begin by writing the general 3x3 Seebeck tensor in matrix form as $\overleftrightarrow{S}=\left(\begin{array}{cc} S & 0\\ 0 & 0 \end{array}\right)$ where $S$ is a 2D symmetric, constant matrix expressible in cartesian form as:$$
S\equiv\left(\begin{array}{cc}
S_{xx} & S_{xy}\\
S_{xy} & S_{yy}
\end{array}\right)=\left(\begin{array}{cc}
\frac{1}{2}\left(t+\left(t^{2}-4d\right)^{\frac{1}{2}}\cos2\rho\right) & -\frac{1}{2}\left(t^{2}-4d\right)^{\frac{1}{2}}\sin2\rho\\
-\frac{1}{2}\left(t^{2}-4d\right)^{\frac{1}{2}}\sin2\rho & \frac{1}{2}\left(t-\left(t^{2}-4d\right)^{\frac{1}{2}}\cos2\rho\right)
\end{array}\right)
\tag{2}$$The right-most expression in eq.(2) displays the matrix in terms of its two rotational-invariants $t\equiv\text{Tr}S,d\equiv\text{Det}S$ and a rotation-angle $\rho$ that parametrizes the orientation of the cartesian axes $x,y$ to which $S$ is referenced. Of particular interest are the two specific values of the rotation angle $\rho=0$ (the "Diagonal" case) and $\rho=-\frac{\pi}{4}$ (the "Alternative" case):$$
\left.S\right|_{\rho=0}\equiv S_{D}=\left(\begin{array}{cc}
\frac{1}{2}t+\frac{1}{2}\left(t^{2}-4d\right)^{\frac{1}{2}} & 0\\
0 & \frac{1}{2}t-\frac{1}{2}\left(t^{2}-4d\right)^{\frac{1}{2}}
\end{array}\right),\;\left.S\right|_{\rho=-\frac{\pi}{4}}\equiv S_{A}=\left(\begin{array}{cc}
\frac{1}{2}t & \frac{1}{2}\left(t^{2}-4d\right)^{\frac{1}{2}}\\
\frac{1}{2}\left(t^{2}-4d\right)^{\frac{1}{2}} & \frac{1}{2}t
\end{array}\right)
\tag{3a,b}$$$S_D$ is the most common representation of $S$, in which the orientation of the $x,y$ axes has been chosen to zero-out the off-diagonal term (this is the form you used and it's the source of your problem). $S_A$ is a less common (alternative) form of $S$ but it contains exactly the same information as $S_D$ and using it resolves your issue, as I now demonstrate.
Taking the curl of eq.(1) and projecting its $z$-component gives:$$\hat{k}\cdot\nabla\times\vec{J}\equiv\hat{k}\cdot\left(\nabla\times\nabla\times\psi\,\hat{k}\right)=-\nabla^{2}\psi=\hat{k}\cdot\left(\sigma\nabla\times\overleftrightarrow{S}\left(-\nabla T\right)\right)\equiv\omega\left(T\right)\tag{4}$$I refrain here from displaying the general form of the source-term $\omega(T)$, but I should at least mention that in polar coordinates it becomes a linear combination of all 5 first and second derivatives of $T$ w.r.t. $r,\theta$. Instead, I specialize to your simple temperature profile $T\left(r,\theta\right)=T_{0}+\left(\frac{2\Delta T}{\pi}\right)\theta$ (which obviously satisfies both the heat equation $\nabla^2T=0$ and your boundary conditions). The resulting equation for $\psi$ is:$$\nabla^{2}\psi=\left(\frac{2\sigma\Delta T}{\pi r^{2}}\right)\left[\left(s_{xx}-s_{yy}\right)\cos\left(2\theta\right)+2s_{xy}\sin\left(2\theta\right)\right]=\left(\frac{2\sigma\Delta T}{\pi r^{2}}\right)\left(t^{2}-4d\right)^{\frac{1}{2}}\cos\left(2\left(\theta+\rho\right)\right)\tag{5}$$This equation is now manifestly invariant under 2D rotations: the scalar $\psi$ is sourced by the scalar term $\left(t^{2}-4d\right)/r^2$ built from the Seeback tensor. And eq.(5) verifies your work: if, as you assumed, $s_{xy}$ vanishes (i.e., $\rho=0$), it agrees exactly with the result of your derivation, but of course it's also inconsistent with your boundary conditions for $\psi$, which leads to your impasse. So instead, I encourage you to "go alternative" and adopt the form $S_A$ (i.e., $\rho=-\frac{\pi}{4}$) for the Seebeck matrix so you may solve the consistent equation:$$\nabla^{2}\psi=\left(\frac{2\sigma\Delta T}{\pi r^{2}}\right)\left(t^{2}-4d\right)^{\frac{1}{2}}\sin\left(2\theta\right)\tag{6}$$

 

[fluidistic]

Thanks for this clarification! I think I can now resolve your problem. As I learned from a modicum of googling, your current:$$\vec{J}\equiv\nabla\times\psi\,\hat{k}=\sigma\left(-\nabla V\right)+\sigma\overleftrightarrow{S}\left(-\nabla T\right)\tag{1}$$is the sum of the usual ohmic electric current, induced by a voltage gradient in a medium characterized by a scalar conductivity $\sigma$, with a thermoelectric current induced by a thermal gradient, as characterized by a symmetric tensor $\overleftrightarrow{S}$ of Seebeck coefficients. With (1) now properly expressed as a 3-vector equation, I can use standard coordinate-transformation rules to convert it to polar coordinates.
Since your problem is independent of the $z$-coordinate, I begin by writing the general 3x3 Seebeck tensor in matrix form as $\overleftrightarrow{S}=\left(\begin{array}{cc} S & 0\\ 0 & 0 \end{array}\right)$ where $S$ is a 2D symmetric, constant matrix expressible in cartesian form as:$$
S\equiv\left(\begin{array}{cc}
S_{xx} & S_{xy}\\
S_{xy} & S_{yy}
\end{array}\right)=\left(\begin{array}{cc}
\frac{1}{2}\left(t+\left(t^{2}-4d\right)^{\frac{1}{2}}\cos2\rho\right) & -\frac{1}{2}\left(t^{2}-4d\right)^{\frac{1}{2}}\sin2\rho\\
-\frac{1}{2}\left(t^{2}-4d\right)^{\frac{1}{2}}\sin2\rho & \frac{1}{2}\left(t-\left(t^{2}-4d\right)^{\frac{1}{2}}\cos2\rho\right)
\end{array}\right)
\tag{2}$$The right-most expression in eq.(2) displays the matrix in terms of its two rotational-invariants $t\equiv\text{Tr}S,d\equiv\text{Det}S$ and a rotation-angle $\rho$ that parametrizes the orientation of the cartesian axes $x,y$ to which $S$ is referenced. Of particular interest are the two specific values of the rotation angle $\rho=0$ (the "Diagonal" case) and $\rho=-\frac{\pi}{4}$ (the "Alternative" case):$$
\left.S\right|_{\rho=0}\equiv S_{D}=\left(\begin{array}{cc}
\frac{1}{2}t+\frac{1}{2}\left(t^{2}-4d\right)^{\frac{1}{2}} & 0\\
0 & \frac{1}{2}t-\frac{1}{2}\left(t^{2}-4d\right)^{\frac{1}{2}}
\end{array}\right),\;\left.S\right|_{\rho=-\frac{\pi}{4}}\equiv S_{A}=\left(\begin{array}{cc}
\frac{1}{2}t & \frac{1}{2}\left(t^{2}-4d\right)^{\frac{1}{2}}\\
\frac{1}{2}\left(t^{2}-4d\right)^{\frac{1}{2}} & \frac{1}{2}t
\end{array}\right)
\tag{3a,b}$$$S_D$ is the most common representation of $S$, in which the orientation of the $x,y$ axes has been chosen to zero-out the off-diagonal term (this is the form you used and it's the source of your problem). $S_A$ is a less common (alternative) form of $S$ but it contains exactly the same information as $S_D$ and using it resolves your issue, as I now demonstrate.
Taking the curl of eq.(1) and projecting its $z$-component gives:$$\hat{k}\cdot\nabla\times\vec{J}\equiv\hat{k}\cdot\left(\nabla\times\nabla\times\psi\,\hat{k}\right)=-\nabla^{2}\psi=\hat{k}\cdot\left(\sigma\nabla\times\overleftrightarrow{S}\left(-\nabla T\right)\right)\equiv-\lambda\left(T\right)\tag{4}$$I refrain here from displaying the general form of the source-term $\lambda(T)$, but I should at least mention that in polar coordinates it becomes a linear combination of all 5 first and second derivatives of $T$ w.r.t. $r,\theta$. Instead, I specialize to your simple temperature profile $T\left(r,\theta\right)=T_{0}+\left(\frac{2\Delta T}{\pi}\right)\theta$ (which obviously satisfies both the heat equation $\nabla^2T=0$ and your boundary conditions). The resulting equation for $\psi$ is:$$\nabla^{2}\psi=\left(\frac{2\sigma\Delta T}{\pi r^{2}}\right)\left[\left(s_{xx}-s_{yy}\right)\cos\left(2\theta\right)+2s_{xy}\sin\left(2\theta\right)\right]=\left(\frac{2\sigma\Delta T}{\pi r^{2}}\right)\left(t^{2}-4d\right)^{\frac{1}{2}}\cos\left(2\left(\theta+\rho\right)\right)\tag{5}$$This equation is now manifestly invariant under 2D rotations: the scalar $\psi$ is sourced by the scalar term $\left(t^{2}-4d\right)/r^2$ built from the Seeback tensor. And eq.(5) verifies your work: if, as you assumed, $s_{xy}$ vanishes (i.e., $\rho=0$), it agrees exactly with the result of your derivation, but of course it's also inconsistent with your boundary conditions for $\psi$, which leads to your impasse. So instead, I encourage you to "go alternative" and adopt the form $S_A$ (i.e., $\rho=-\frac{\pi}{4}$) for the Seebeck matrix so you may solve the consistent equation:$$\nabla^{2}\psi=\left(\frac{2\sigma\Delta T}{\pi r^{2}}\right)\left(t^{2}-4d\right)^{\frac{1}{2}}\sin\left(2\theta\right)\tag{6}$$

I thank you in the most honest way I can. I will try to go through your derivations, including the ones you didn't write down explicitely here.

I have 2 questions, please tell me what you think.
1: Does choosing rho equals pi/4 corresponds to having the crystallographic axis rotated by pi/4 rad compared to my rho equals 0 case? I am 99.99 percent sure it does, I,would like your confirmation. This indeed changes the problem quite a bit.
2: Does this imply that the assumption of a linear (in theta) temperature profile is not consistent with having a J that is contained in the quarter of an annulus and which satisfies the cos(2 theta)/r^2 expression? I believe so. What I didn't tell you is that the real heat equation is more complex than the one assumed here. It contains 2 other terms involving J and its spatial derivative. However, there is a parameter involving S_xx -S_yy, which I believe can be considered a very small one so that the heat eq PDE could be tackled by perturbation method. It turns out that this doesn't seem to work here, I am very surprised. Physically it would mean that applying the Delta T on the 2 sides generates a contained current, but it is able to be like so solely because T departs from a linear T(theta). The departure from linearity should ''fix'' the mathematical problem and this should result in a T(r, theta) that creates a consistent J that satisfies the curl relationship as well as the boundary conditions. I am still skeptical about this conclusion. I will write down the full heat equation. If you could tell me your thoughts, I would be enlightenend in my understanding of this problem.

$\nabla \cdot (\kappa \nabla T) +\frac{|\vec J|^2}{\sigma}-T\left ( S_{xx}\frac{\partial J_y}{\partial x}+S_{yy}\frac{\partial J_y}{\partial y}\right ) =0$.

I wanted to obtain an analytical expression for J (an approximation was the goal because the full problem is too complex) to get some insights and check whether some relations hold true. I can already test if they hold true in the case where the inner and outer radii are kept at T fixed.

This would go against a conclusion in what's written in a paper but I should check the details and assumptions. I suspect that the second term of.the heat equation should be equal in magnitude but opposite in sign to the 3rd term, when integrated over the whole material (therefore surface). To me, this is the only way for.the 1st law of thermodynamics to hold, but a paper by Lukosz (1964) claims that the 3rd term is dominant.


Edit: I should also mention that visually, my numerical solution yields a T distribution that is ''very close'' to a linear one. This is why I am quite surprised that the departure from linearity is what impacts J in such a way as to make it consistent with the curl expression as well as the b.c.s.

 

[renormalize]

1: Does choosing rho equals pi/4 corresponds to having the crystallographic axis rotated by pi/4 rad compared to my rho equals 0 case? I am 99.99 percent sure it does, I,would like your confirmation. This indeed changes the problem quite a bit.
2: Does this imply that the assumption of a linear (in theta) temperature profile is not consistent with having a J that is contained in the quarter of an annulus and which satisfies the cos(2 theta)/r^2 expression? I believe so.

Yes, choosing $\rho=-\pi/4$ means rotating the crystallographic axis by $45^{\circ}$, which may not make physical sense in your problem. So how about rotating your annular region instead? In other words, rotate your quarter-annulus by $45^{\circ}$ and impose your boundary conditions at $\theta=\pi/4$ and $\theta=3\pi/4$.

 

[fluidistic]

Yes, choosing $\rho=-\pi/4$ means rotating the crystallographic axis by $45^{\circ}$, which may not make physical sense in your problem. So how about rotating your annular region instead? In other words, rotate your quarter-annulus by $45^{\circ}$ and impose your boundary conditions at $\theta=\pi/4$ and $\theta=3\pi/4$.

I mean, it does make physical sense, but it's a different problem. Physically, it means preparing a sample cut differently with regards to its crystallographic axis, and placing it in the same position and shape as in my original problem. This certainly impacts the current density created in it, but it makes physical sense.

About a pure rotation of my original problem, I am missing to see the point. Isn't the problem totally equivalent to the original one? I do not see how the math "gets fixed" at all.

I am thinking to take your $\nabla^{2}\psi=\left(\frac{2\sigma\Delta T}{\pi r^{2}}\right)\left(t^{2}-4d\right)^{\frac{1}{2}}\cos\left(2\left(\theta+\rho\right)\right)\tag{5}$ expression, and solve it for the general case with $\rho >0$, I am not sure if it a feasible task.
For instance, I do not get a nice simplification when $\theta \neq \pi/4$.

I wasn't able yet to reach your expressions, and I am not satisfied with the $\lambda(T)$ I get. I get those 5 terms, but the problem is that my expression vanishes when $T$ is linear in $\theta$, so I probably goofed somewhere.

 

[renormalize]

I wasn't able yet to reach your expressions, and I am not satisfied with the $\lambda(T)$ I get. I get those 5 terms, but the problem is that my expression vanishes when $T$ is linear in $\theta$, so I probably goofed somewhere.

I fear this discussion has become rather narrowly focused and probably lacks broad interest. I suggest you ask a moderator to close this thread and then we can continue our discussion via PF messaging.