I Stuck on solving a non homogenous (linear) PDE

  • Thread starter Thread starter fluidistic
  • Start date Start date
  • Tags Tags
    Partial Pde
  • #31
fluidistic said:
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}$$
 
Last edited:
Physics news on Phys.org
  • #32
renormalize said:
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.
 
Last edited:
  • #33
fluidistic said:
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##.
 
  • #34
renormalize said:
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.
 
  • #35
fluidistic said:
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.
 
  • #36
renormalize said:
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.
In this case can't we just continue in private without closing the thread, just in case someone else pops up with insights we could have missed? It will be closed after 2 years anyway.
 
  • #37
I've resolved the question posed by the OP so I'm posting my solution to benefit future readers of this thread.
@fluidistic seeks an analytic solution to the 2D Poisson equation:$$-\nabla^{2}\psi\equiv-\frac{\partial^{2}\psi}{\partial r^{2}}-\frac{1}{r}\frac{\partial\psi}{\partial r}-\frac{1}{r^{2}}\frac{\partial^{2}\psi}{\partial\theta^{2}}=\frac{A\cos2\theta}{r^{2}}\tag{1}$$where ##A## is a constant and the stream-function ##\psi\left(r,\theta\right)## must vanish on the boundary of the quarter-annulus; i.e.:$$\psi\left(r_{i},\theta\right)=\psi\left(r_{o},\theta\right)=\psi\left(r,0\right)=\psi\left(r,\frac{\pi}{2}\right)=0\tag{2a,b,c,d,}$$Evidently the particular solution ##\psi=\psi_{p}\equiv\frac{1}{4}A\cos 2\theta## satisfies (1) but, alas, it fails all of the above boundary conditions (2). But since we're free to superpose ##\psi_p## with any homogeneous solution ##H##, where ##\nabla^{2}H=0##, we choose to add ##H=h\left(\theta\right)\equiv\frac{1}{4}A\left(\frac{4\theta}{\pi}-1\right)## to define a new particular solution:$$\psi_{p^{\prime}}\equiv\psi_{p}+h\left(\theta\right)=\frac{1}{4}A\left(\cos2\theta+\frac{4\theta}{\pi}-1\right)\tag{3}$$that does satisfy the two angular conditions (2c,2d). Next, we Fourier expand (3) into a sine series to get:$$\psi_{p^{\prime}}=\frac{A}{2\pi}\sum_{n=1}^{\infty}\frac{\sin4n\theta}{n\left(4n^{2}-1\right)}\tag{4}$$Finally, noting that ##H=h_{4n}\left(r,\theta\right)\equiv\frac{\sin4n\theta\left(r^{8n}+r_{i}^{4n}r_{o}^{4n}\right)}{r^{4n}\left(r_{i}^{4n}+r_{o}^{4n}\right)}## is a homogeneous solution for any integer ##n##, we make the substitution:$$\sin4n\theta\rightarrow\sin4n\theta-h_{4n}\left(r,\theta\right)=\sin4n\theta-\frac{\sin4n\theta\left(r^{8n}+r_{i}^{4n}r_{o}^{4n}\right)}{r^{4n}\left(r_{i}^{4n}+r_{o}^{4n}\right)}=-\frac{\sin4n\theta\left(r^{4n}-r_{i}^{4n}\right)\left(r^{4n}-r_{o}^{4n}\right)}{r^{4n}\left(r_{i}^{4n}+r_{o}^{4n}\right)}\tag{5}$$into (4) to arrive at the full analytic solution for the stream-function:$$\psi\left(r,\theta\right)=-\frac{A}{2\pi}\sum_{n=1}^{\infty}\frac{\sin4n\theta\left(r^{4n}-r_{i}^{4n}\right)\left(r^{4n}-r_{o}^{4n}\right)}{n\left(4n^{2}-1\right)r^{4n}\left(r_{i}^{4n}+r_{o}^{4n}\right)}\tag{6}$$This manifestly satisfies all four radial and angular boundary conditions (2).
Because performing practical computations necessitates keeping only a finite number of terms, we denote by ##\psi_N## the truncation of the infinite series in (6) to order ##N## and examine the form of the analytic stream-function for modest values of ##N##. For example, setting ##N=4## yields the approximate analytic solution:
\begin{align}
\psi\left(r,\theta\right)\approx\psi_{4}\left(r,\theta\right) & = -\frac{A}{2\pi}\left[\frac{\sin4\theta\left(r^{4}-r_{i}^{4}\right)\left(r^{4}-r_{o}^{4}\right)}{3r^{4}\left(r_{i}^{4}+r_{o}^{4}\right)}+\frac{\sin8\theta\left(r^{8}-r_{i}^{8}\right)\left(r^{8}-r_{o}^{8}\right)}{30r^{8}\left(r_{i}^{8}+r_{o}^{8}\right)}\right. \nonumber \\
& \left.+\frac{\sin12\theta\left(r^{12}-r_{i}^{12}\right)\left(r^{12}-r_{o}^{12}\right)}{105r^{12}\left(r_{i}^{12}+r_{o}^{12}\right)}+\frac{\sin16\theta\left(r^{16}-r_{i}^{16}\right)\left(r^{16}-r_{o}^{16}\right)}{252r^{16}\left(r_{i}^{16}+r_{o}^{16}\right)}\right] \nonumber\tag{7}
\end{align}A comparison of the finite-element numerical solution ##\psi_{FE}## of eqs.(1,2) to the analytic-approximation ##\psi_4## looks like this:
1731711348095.png

Thus, ##\psi_4## offers a reasonably good analytic approximation to the full numerical stream-function ##\psi_{FE}##.
 
Last edited:
  • Like
  • Love
Likes WWGD, Spinnor and fluidistic

Similar threads

  • · Replies 3 ·
Replies
3
Views
4K
  • · Replies 4 ·
Replies
4
Views
4K
  • · Replies 4 ·
Replies
4
Views
2K
  • · Replies 1 ·
Replies
1
Views
1K
  • · Replies 3 ·
Replies
3
Views
3K
  • · Replies 5 ·
Replies
5
Views
2K
  • · Replies 5 ·
Replies
5
Views
3K
  • · Replies 13 ·
Replies
13
Views
3K
  • · Replies 2 ·
Replies
2
Views
3K
  • · Replies 3 ·
Replies
3
Views
816