convection diffusion with time-space dependent coefficients

[aao24]

Hello! I would like to solve a convection diffusion equation with a point source and time-and-space-dependent coefficients $D=D(x,t)$, u=(x,t), $x\in R^d$, $d=1,2,3.$

Is it possible to use ode45 or similar solver instead of doing loop in time if have time-dependent coefficients? Or would you recommend to use time-loop?

I tried to modify example diffusionODEtutorial but unfortunately it does not work. The example convectionODEexample does not work either...

Best regards, Anna

[simulkade]

Hi Anna,

Thanks for reporting this. The examples (and the code) are fixed now. ode45 uses adaptive time stepping which is superior than the Euler method that I use for time stepping, using transientTerm function. But implementing it in FVTool can be a bit dirty sometimes. Try to modify the examples and let me know if I can help.

[aao24]

Thanks a lot!

[aao24]

RefelctingBoundaries.pdf Ali, I am now solving convection (diffusion turned out to be neglectable) equation in 2D on a rectangular grid with the time-space dependent velocity field. I use TVD for convection term. It seem to work very good until I reach the boundary. Is there some simple way to address the reflecting boundaries in your toolbox? Basically I just want the tracer to freely flow away from the simulated region...

Thanks, Anna

[simulkade]

Hi Anna, I don't know anything about the reflecting boundaries. Do you have any simple formulation for it, using ghost cells?

[aao24]

Hi Ali, By "reflecting boundaries" I meant that the boundary start to behave like a "source" as soon as we reach them. This is opposite of what I want- I want the contaminant to flow freely out of the boundaries. I was trying to model this with various boundary conditions, but not successfully. Since I am quite new to this business, I wonder if you have any advice for me...

Best, Anna

[simulkade]

If you are solving an advection equation, the default Neumann BC should do the trick. I do not know what goes wrong here unless I see the code.

On Mon, Nov 6, 2017 at 2:34 PM, aao24 notifications@github.com wrote:

Hi Ali, By "reflecting boundaries" I meant that the boundary start to behave like a "source" as soon as we reach them. This is opposite of what I want- I want the contaminant to flow freely out of the boundaries. I was trying to model this with various boundary conditions, but not successfully. Since I am quite new to this business, I wonder if you have any advice for me...

Best, Anna

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/simulkade/FVTool/issues/15#issuecomment-342149731, or mute the thread https://github.com/notifications/unsubscribe-auth/ABy1xhMOewFCC4rGyTPMBEeL5eZi7nxUks5szwragaJpZM4P_CbN .

[dankirshner]

If I understand it correctly, the Neumann BC code sets a boundary to a specific concentration -- which I believe requires that a flux is allowed through that boundary.

To me, a "reflective boundary" is a zero-flux boundary -- which means one does NOT specify the concentration there.

I do want to specify some zero-flux boundaries. I've been looking through the examples, and haven't figured it out. "fluxLimiter" sounds like the right thing. Am I on the right track?

Thanks very much!

--Dan

[simulkade]

Hi Dan,

You can specify the flux with a Neumann bc. The concentration is set with Dirichlet. In fvtool, the default is Neumann with zero flux. You can also set manually. I'll add a sample code when from my computer.

On Nov 8, 2017 21:01, "dankirshner" notifications@github.com wrote:

If I understand it correctly, the Neumann BC code sets a boundary to a specific concentration -- which I believe requires that a flux is allowed through that boundary.

To me, a "reflective boundary" is a zero-flux boundary -- which means one does NOT specify the concentration there.

I do want to specify some zero-flux boundaries. I've been looking through the examples, and haven't figured it out. "fluxLimiter" sounds like the right thing. Am I on the right track?

Thanks very much!

--Dan

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/simulkade/FVTool/issues/15#issuecomment-342941465, or mute the thread https://github.com/notifications/unsubscribe-auth/ABy1xkkn7Ft7xKZ7HIsmz4e7QaZ5zKO8ks5s0giKgaJpZM4P_CbN .

[simulkade]

The Neumann bc (fixed flux) is defined by

bc = createBC(m)
bc.left.a(:) = 1.0 % transfer coefficient
bc.left.b(:) = 0.0 % must be zero for Neumann
bc.left.c(:) = 0.0 % flux value

[dankirshner]

Thanks very much indeed! I've discovered that the comments in FVTdemo.m are very helpful -- it would be nice to have a prominent pointer directing people to that file. It would also be really nice to have some collected documentation -- I have hopes that I could help with that. This is looking like a great tool, so I'd like to help!

In the meantime, I've been trying to get a 1D diffusion-advection case working with a fixed concentration on the right and a zero-flux condition on the left. It seems to work as expected for some (high) values for the (uniform) convection velocity, but not for other (lower) values.

The code is here: dk1D_conv.m.txt

I run a time-course of 1D concentration profiles, where the initial concentration is zero except that the right boundary condition fixed the concentration = 1.0. (The length is 1.0e-9 m - I'm doing diffusion of a small molecule into a protein pocket!)

Here, with velocity (u) set to 2.0, the concentration profile converges to a steady state (the red line at the bottom is the initial time step, subsequent time steps are above), where, with the convection "pulling" things to the right, I expect an exponential shape (per the constant velocity).

But, when I reduce the velocity to 1, I still expect to see an exponentially-shaped steady state. But instead the concentration profile converges to a constant 1.0. (This is with dt = 500e-12 and final_t = 10000e-12.)

It seems like as soon as a non-zero concentration "touches" the left boundary, then things are different! I thought that maybe I needed to set the velocity at the left-most face (u_face.xvalue(1)) to zero, but I haven't figured out how to do that (it appears that createFaceVariable(m, u) uses only the first element of a vector u).

Any assistance you could give me would be greatly appreciated.

---

Update after sleeping on it...

Adding

u_face.xvalue(1) = 0

after doing createFaceVariable(m, u) did the trick!

The result for u = 1 (except the left boundary) is now

[simulkade]

Hi Dan,

Please excuse me for my late reply. I was extremely busy during the past few weeks. Thanks for your comments on the FVTool. It would be amazing if you can help with the documentation. Please let me know if you need any help from my side.

I have difficulty understanding the physics of the PDE you have formulated. It seems that you have an advection term that acts against diffusion, and since the domain is small the diffusion term is very important indeed.

To me it makes sense that when the velocity is small, the diffusion "wins" and the concentration in the domain becomes equal to the concentration of the right boundary. And of course when the velocity for the advection term is large, it pushes the diffusive front back and finally the system reaches an steady state with a certain concentration profile. Therefore, I don't see any reason for not trusting your result nor any need to change your velocity value at the boundary.

Edit: when I run your script, I get that steady state for velocities larger than 8.0, but then you have to increase the simulation time.