Plane Poiseuille Flow Vector form BC problem

[leo mangalath]

Hi all,

I am trying to write my first Dedalus script for a 2D pressure driven channel flow.  I implemented the pressure difference by adding a constant to the pressure gauge condition.

problem.add_equation("integ(p)=dp_dx")

I have initialized the flow field with a velocity of 0, so I was expecting to see it rising to the analytical velocity distribution (Re=177). However, the velocity does not rise and is practically zero throughout the timesteps. Furthermore the pressure is constant throughout the domain, which leads me to believe that I implemented it wrong. I tried looking for similar cases but most cases are from dedalus v2 or not in vector form. I'm sure its a silly mistake but I am very grateful for every help. I have attached my python script below.

Thanks in advance!

Best Leo,

[Ben Brown]

Leo,

     You've added a constant offset to the pressure, everywhere, rather than adding a constant pressure gradient forcing term in the momentum equation.

The best way to do this is to probably prescribe:

dp_dx = <appropriate # value>*ex

where ex is the unit vector pointing in the x direction

and then

dt(u) + grad(p) + ... <other linear terms> = - dp_dx +... <other not-linear terms>

Note that as a not-linear term (a constant), the pressure forcing shows up as an explicit (RHS) term in the equations.

I've done this before for poiseuille flow and similar, and my memory is that it all generally works great.  There may be some details or approaches to doing this better, but this should at least get you started.

--Ben

--
You received this message because you are subscribed to the Google Groups "Dedalus Users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to dedalus-user...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/c0bbf2dc-a0e7-4aa1-8a56-304bb5ad18f9n%40googlegroups.com.

[leo mangalath]

Hi Ben,

I think I made a mistake in sending the last reply. Here is my reply again:

I implemented your suggestion and it did work, thank you! Now I want to replace the constant pressure gradient term with a forcing which is dependent on the flow rate. Furthermore due to the definition of the flow rate, the bulk velocity will have to change accordingly. My approach is to add an additional equation regarding the flow rate definition, where the bulk velocity is calculated by the first coefficient of the fourier series (ubulk=u['c][0][0][0], xbasis,1st node, 1 coefficient).  However, when comparing the a0 coefficient with a manual averaging of the velocity in x direction at node 1, I get different results. I am unsure where my error comes from.

Best,

Leo

[Ben Brown]

Leo [on bcc] and dedalus-users [in case anyone else is grappling with similar in the future],

     I would suggest calculating averages via integration using our integration operators.  That's much more robust, as you otherwise need to exercise real care when dealing with e.g., the chebshev polynomial grid.

Here's an example for calculating these sorts of things in a 2-D cartesian domain:

vol = Lx*Lz

integ = lambda A: de.Integrate(de.Integrate(A, 'x'), 'z') 

avg = lambda A: integ(A)/vol 

x_avg = lambda A: de.Integrate(A, 'x')/(Lx)

"avg" is the volume average, while "x_avg" calculates the average in the horizontal (here Fourier) direction.

What you extracted previously via the coefficients is the zero mode in fourier, but also the zero mode in chebyshev (assuming chebyshev vertical grid), and the quadrature weighting on chebyshev coefficients is different and difficult to turn into other relevant quantities like averages.  The integral operators make this robust and straightforward.

--Ben

On Thu, Sep 28, 2023 at 9:39 AM leo mangalath <leo.ma...@hotmail.de> wrote:

Hi Ben,

Thank you for your quick reply. I have implemented your suggestion in my script. Furthermore, I think it is more convenient for me if I replace the constant pressure gradient with a forcing induced by a constant flow rate. I am planning to do a flow rate which is dependent on the bulk velocity. I tried to calculate the bulk velocity by simply extracting the first Fourier coefficient (by a0=u['c'][0][0][0] (xbasis, 1st node, 1 coefficient). However, when manually calculating the average velocity the results differ. I am still unsure where this error comes from.

Best,

Leo

[Ben Brown]

Leo & all,

     I would suggest calculating averages via integration using our integration operators.  That's much more robust, as you otherwise need to exercise real care when dealing with e.g., the chebshev polynomial grid.

Here's an example for calculating these sorts of things in a 2-D cartesian domain:

vol = Lx*Lz

integ = lambda A: de.Integrate(de.Integrate(A, 'x'), 'z') 

avg = lambda A: integ(A)/vol 

x_avg = lambda A: de.Integrate(A, 'x')/(Lx)

"avg" is the volume average, while "x_avg" calculates the average in the horizontal (here Fourier) direction.

What you extracted previously via the coefficients is the zero mode in fourier, but also the zero mode in chebyshev (assuming chebyshev vertical grid), and the quadrature weighting on chebyshev coefficients is different and difficult to turn into other relevant quantities like averages.  You can see these terms in e.g., Burns et al 2020, Section IV.C (the dedalus instrument paper).  The integral operators make taking an average robust and straightforward.

Best,

--Ben

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/ce0f1065-5ab5-4128-999f-b57eaa371067n%40googlegroups.com.

[leo mangalath]

Hi Ben,
I have finally managed to solve the poiseuille flow case and have verified it with the analytical solution. However, I noticed that my simulations take very long around 15 min for a N=2x64 2D case. I have tried reducing the resolution (to 8x8) and keeping the problem parameters to a minimum but to no avail. I am sure it is just a small misunderstanding of one sim config parameter but I cannot seem to find it. Any help would be greatly appreciated.

Best,
Leo

[leo mangalath]

I just realized that I mistakenly defined a very small channel geometry, which resulted into a very high velocity, which slowed the simulation down. After changing the dimensions to get a more reasonable velocity the simulation just took a few seconds.

[Ben Brown]

Leo,

     Brilliant.  Thanks for the update; I was suspecting something like that.

--Ben

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/c2f5d364-82ee-4098-a0ba-70352100701dn%40googlegroups.com.