2D channel flow with periodic BCs: assignment of constant pressure gradient

[alessandr...@gmail.com]

Dear all,

with the aim to simulate a three-dimensional fully developed turbulent channel flow, I am starting with a simple laminar 2D case.

The computational domain is a Lx x Ly rectangle (Ly = 2H, Lx = 14H, where H is the channel half-width) that I am simulating with the "multigrid.h" library.

1) As a preliminary check, I set up a simulation with the following boundary conditions (BCs ):

/**
The top and bottom walls are no-slip. */

u.n[top] = dirichlet(0.);
u.t[top] = dirichlet(0.);
u.n[bottom] = dirichlet(0.);
u.t[bottom] = dirichlet(0.);

/**
 Velocity inlet on left side*/
u.n[left] = y < my_Ly ? dirichlet(-my_U_0*sq(y-0.5*my_Ly)/sq(0.5*my_Ly)+my_U_0) : dirichlet(0);
u.t[left] = dirichlet(0.);

/**
 Pressure outlet on right side*/
p[right]=dirichlet(0.);
u.n[right] = u.n[] > 0 ? neumann(0.): 0;
u.t[right] = neumann(0.);

RESULT:

Results of the simulation are attached as x-velocity component (u_inlet.mp4) and pressure (p_inlet.mp4) videos.

As I am assigning a Poiseuille velocity inflow as left BC, the simulation converges to the (exact?) Poiseuille laminar solution in the whole channel.

2) Now I want to assign periodic boundary conditions on the left side, and insert a constant pressure gradient as a source (acceleration) term in the momentum equation. Therefore, I assign:

/**
The top and bottom walls are no-slip. */

u.n[top] = dirichlet(0.);
u.t[top] = dirichlet(0.);
u.n[bottom] = dirichlet(0.);
u.t[bottom] = dirichlet(0.);

int main() {
  dimensions (nx = 7);  
  periodic(left);
  init_grid (1 << my_initlev);
  size (my_Lx);
  mu = muv;

  TOLERANCE = 1e-4;
  run();
}

/* The pressure gradient required to drive the flow is introduced as a
   source term */

event acceleration (i++) {
  face vector av = a;
  foreach_face()
  av.x[] += my_dpdx;
}

RESULT:

The simulation fails before the first iteration with the following error message attached (out_fail).

What I found strange is that, if I remove the acceleration event, the simulation runs fine and converges to the solution attached as u_periodic.mp4 and p_periodic.mp4. What happens in this case is that the mass flow rate goes to zero, I think because there is no pressure gradient able to sustain the flow and BCs are periodic.

Therefore, my question is: is there any "physical" and/or "numerical" incompatibility in assigning a constant acceleration along a spatial direction where period BCs are imposed?

I attach also two versions of my code:

1) working.c --> periodic BCs with no pressure gradient

2) not_working.c --> same as working.c but with a constant pressure gradient assigned.

Thank you in advance for your help.

Cheers,

Alessandro

[Nicolò Scapin]

Hi Alessandro,

I had quick check but I notice something strange here

event acceleration (i++) {
  face vector av = a;
  foreach_face()
  av.x[] += my_dpdx;
}

You should set the pressure gradient only in the x-direction. So

event acceleration (i++) {
  face vector av = a;

  foreach_face(x)
  av.x[] += my_dpdx;
}

Can you try if this fix your issue?

Best,

Nicolo

[Alessandro Della Pia]

Dear Nicolo,

Thank you very much for your reply.

You are right.

Anyway, I changed the acceleration event as follows

event acceleration (i++) {
  face vector av = a;
  foreach_face(x)
  av.x[] += my_dpdx;
}

but I get the same failure as before.

The curious fact is that the simulation fails also if I set the source term equal to zero as follows

event acceleration (i++) {
  face vector av = a;
  foreach_face(x)

  av.x[] += 0;
}

Moreover, the same happens also if periodic BCs are replaced with velocity_inlet-pressure_outlet PCs described above in the discussion.

Is there any incompatibility between the acceleration event that I made and the "navier-stokes/centered.h” library?

Best,

Alessandro

-- 
You received this message because you are subscribed to a topic in the Google Groups "basilisk-fr" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/basilisk-fr/oFCHO6ph-D4/unsubscribe.
To unsubscribe from this group and all its topics, send an email to basilisk-fr...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/basilisk-fr/5d2a37ad-4f93-4665-8ccc-8650c23f4140n%40googlegroups.com.

[Yi Dai]

Hi, 

Could you try with 

face vector av[] instead of face vector av

Dai

[Yi Dai]

sorry forget to mention this: 

make it a global variable

[Shyam Sunder Yadav]

Dear Alessandro

Do you have to add something like shown below?

  if (is_constant (a.x))
    a = new face vector;

Please check this.

Thanks

Regards

--
You received this message because you are subscribed to the Google Groups "basilisk-fr" group.
To unsubscribe from this group and stop receiving emails from it, send an email to basilisk-fr...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/basilisk-fr/c7fa87a0-47e2-42b3-a889-0f9a622feaa7n%40googlegroups.com.

--

Dr. Shyam Sunder Yadav
Associate Professor
Mechanical Engineering
BITS Pilani
09902346342
http://www.bits-pilani.ac.in/pilani/ssyadav/Profile

The information contained in this electronic communication is intended solely for the individual(s) or entity to which it is addressed. It may contain proprietary, confidential and/or legally privileged information. Any review, retransmission, dissemination, printing, copying or other use of, or taking any action in reliance on the contents of this information by person(s) or entities other than the intended recipient is strictly prohibited and may be unlawful. If you have received this communication in error, please notify us by responding to this email or telephone and immediately and permanently delete all copies of this message and any attachments from your system(s). The contents of this message do not necessarily represent the views or policies of BITS Pilani.

[Shyam Sunder Yadav]

Alessandro

You have the following in centered.h

(const) face vector mu = zerof, a = zerof, alpha = unityf;

"a" is constant face vector, therefore add the following lines in your event init (t = 0){...}

if (is_constant (a.x))
    a = new face vector;

Its working now in my system with this.

Please check.

[Alessandro Della Pia]

Dear Yi, Shyam,

Thank you really much for your precious suggestions.

Based on your comments and on the poiseuille.c file that I found on BASILISK website,

I was able to solve the problem.

I attach the code that is now working.

Thanks again,

Alessandro

[Alessandro Della Pia]

Dear all,

I am trying to simulate a 3D channel flow.

The geometry is Lx x Ly x Lz, where:

Lx=14 dimension along the axis of the channel);

Ly=Lz=2  (normal-to-flow and span wise dimensions).

I am using the multigrid.h library with uniform grid (the code is attached).

Right now, I am able to simulate this geometry with multiple processors only by running the simulation with number of processors np=7.

From my understanding from the website, np=7 means nx=np, ny=nz=1, and so the desired channel is obtained.

But it is not completely clear to me how can I increase the number of processors without changing the geometry of the channel.

For example, can I double the number of processors (from np=7 to np=14) but still simulating the same channel (Lx=14, Ly=Lz=2)?

Thanks in advance for your help,

Alessandro.

[Wouter Mostert]

Dear Alessandro,

I have never tried this, but one way around the funny domain decomposition problem might be to use an octree grid, refined globally to whatever resolution you want on initialization, and without any adapt_wavelet commands. Use embedded boundaries to specify channel geometry. Philosophically it is a waste of amazing AMR capability, but it is functionally multigrid but without the odd MPI-domain-decomposition property.

Wouter