Weird behavior in pdepe at origin

3 次查看(过去 30 天)
I am have fluid dynamics model in spherical coordinates that I am simulating pdepe. It works well, except that I get weird bevahior at the origin, where it seems like pdepe is not maintaining zero-flux boundary conditions.
This figure shows the mass flux in the fluid near the origin at various times in the simulation:
You can see it is not zero at the origin. And if you look closely at the purple and green lines, you can see that if the mass flux was zero at the origin, the slope between r=0 and r=1 would more closely match the slope between r=1 and r=2.
The figure shows the mass density in the fluid near the origin:
These curves should have a zero derivative at the origin, but instead suddenly dip down near r=0. And, because of this, the mass density becomes negative at the origin at 40 ns (the purple curve), which is unphysical.
My questions are:
1) Is this expected/known behavior for pdepe?
2) Does the non-zero flux and/or negative density at the origin affect the solution? I know pdepe is attempting to set zero-flux boundary conditions at the origin, so does the value at the origin actually affect the rest of the solution?
3) Is there any way to fix this?
Thanks for your help!
  3 个评论
Brian Cluggish
Brian Cluggish 2025-4-17
Fair enough. The equations and code are pretty complicated. I'll see if I can reproduce the effect on a simpler model.
William Rose
William Rose 2025-5-16
I agree that something is wrong when you have a negative concentration, and you can;t ignore it.
Maybe you can write the explicit forward-in-time difference equation, with centered spatial diferences on the right hand side. (See here.) You may have to use a quite small value for delta t to achieve stability. But it is probably the easiest difference equation to write and implement. You could at least test things out, even though it may be too ineficient for regular use.
When simulating diffusion-advection PDEs in cylindrical geometry, where v(r,t) the velocity of flow at radius r and time t, I use boundary conditions v(r,t)=0 when r=R (non-slip at the wall), and dv/dr=0 at r=0. And the radial flux=0 at r=0 and at r=R (at the center and at the wall).

请先登录,再进行评论。

回答(0 个)

类别

Help CenterFile Exchange 中查找有关 1-D Partial Differential Equations 的更多信息

产品


版本

R2024b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by