EB Isothermal BC for a supersonic channel case (@Mach 2.0)

[RSuryaNarayan]

Hello everyone, I am trying to simulate supersonic flow inside a channel using PeleC. To do so I have taken the EB-Bluffbody case in the Exec/RegTests folder and made the following changes to the inputs.3d file

1. Made changes to the domain

   geometry.prob_hi     =  15.0  6.0  1.5
   amr.n_cell           =  160   64   16

2. Modified BCs using:

   pelec.lo_bc       =  "Hard"     "NoSlipWall" "Interior"
   pelec.hi_bc       =  "FOExtrap" "NoSlipWall" "Interior"

3. Made changes to the EB to cover the entire domain

   eb2.use_eb2 = 1
   eb2.geom_type = box
   eb2.box_lo = 0.0 0.0 0.0
   eb2.box_hi = 15.0 6.0 1.5 
   eb2.box_has_fluid_inside = 1
   ebd.boundary_grad_stencil_type = 0

4. Changed the problem parameters

   prob.p = 1013250.0
   prob.rho = 0.00116
   prob.vx_in =  70000.0
   prob.vy_in =  0.0
   prob.Re_L = 625.0
   prob.Pr = 0.7

   Also since I want the walls to be isothermal at 300K I left the

   pelec.eb_boundary_T =300
   pelec.eb_isothermal =1

   as is. My question is, is this methodology right for the problem at hand or I need to modify something more to run this case as intended? The results obtained from the above inputs file look something like so: Here are the contours of the Mach number:

There is a nice shock wave and its reflections, however, the walls don't seem to maintain the 300K isothermal boundary condition, rather the temperatures vary along the length of the channel. Hence the whole question on the methodology behind the implementation. Thanks in advance.

[marchdf]

Hi,

Thanks for reaching out! Since this is basically just a channel with no complex EB. I would have started by not using EB at all. And then if you need to impose a specific BC for temperature, I would have used "Hard" BC in the input file for the wall BCs and then edited bcnormal function to impose a no-slip condition and a specified wall temperature. Does that make sense? Removing the EB also means you can use other hydro methods (e.g. PPM) which are less diffusive than MOL. Please let me know if you have any other questions.

[RSuryaNarayan]

Thanks for your reply @marchdf . I was indeed contemplating on those lines as using EB wasn't really necessary for this case. I was stuck with the isothermal BC. The NoSlipWall simply takes care of the no-slip at the channel walls. I have the following questions:

1. How do I distinguish between top-bottom and inflow-outflow faces in bcnormal
2. Do I just set the state vector variable (s_ext[UTEMP] in case of temperature snd s_ext[UMX] for x-momentum) to the required value in the bcnormal? The documentation suggests I redefine and write bcnormal again for all of this.

[marchdf]

So here are some pointers:

1. You can use idir and sgn that are passed in as arguments to bcnormal to determine what face you are on. See for example https://github.com/AMReX-Combustion/PeleC/blob/development/Exec/RegTests/EB-C10/prob.H#L72
2. You basically set all the s_ext components (density, momentum, energies, temp). For a no slip wall you would just reflect velocities s_ext[UMX] = - s_int[UMX], same for y/z. And then set the internal energy, total energy, and temp to what it should be according to the wall BC temp.

You can post the bcnormal code here if you want and I can take a look once you've implemented it.

[RSuryaNarayan]

So based on my understanding so far, here's what I've tried implementing: I set the no-slip and isothermal only for the walls of the top and bottom face using the idir and sgn. The other variables are basically the same as that for this EB-BluffBody case.

bcnormal(
  const amrex::Real* /*x[AMREX_SPACEDIM]*/,
  const amrex::Real* /*s_int[NVAR]*/,
  amrex::Real s_ext[NVAR],
  const int /*idir*/,
  const int /*sgn*/,
  const amrex::Real /*time*/,
  amrex::GeometryData const& /*geomdata*/,
  ProbParmDevice const& prob_parm)
{
  //considering both X and Y are set to "Hard"
  //impose the BC on the top walls:
  if ((idir==1 && sgn==1) || (idir==1 && sgn==-1) //not sure of the numbering scheme adopted here! assumed 1 for y-normal (0 for x and 2 for z) 
  //and sgn toggles between 1 and -1 for positive and negative normals hence (1,1) for Y normal and top (1,-1) for Y normal and bottom
  {
    //set the no-slip BC    
    s_ext[UMX] = -s_int[UMX];
    s_ext[UMY] = -s_int[UMY];
    s_ext[UMZ] = -s_int[UMZ];
    //set the isothermal wall
    s_ext[UTEMP] = 300.0; //for the isothermal wall at 300 K
    //for other variables
    s_ext[URHO] = prob_parm.rho;
    s_ext[UEINT] = prob_parm.rho * prob_parm.eint;
    s_ext[UEDEN] = prob_parm.rho *
                 (prob_parm.eint + 0.5 * (prob_parm.vx_in * prob_parm.vx_in +
                                          prob_parm.vy_in * prob_parm.vy_in));
    for (int n = 0; n < NUM_SPECIES; n++)
      s_ext[UFS + n] = prob_parm.rho * prob_parm.massfrac[n];
  }
  else
  {
  s_ext[URHO] = prob_parm.rho;
  s_ext[UMX] = prob_parm.rho * prob_parm.vx_in;
  s_ext[UMY] = prob_parm.rho * prob_parm.vy_in;
  s_ext[UMZ] = 0.0;
  s_ext[UEINT] = prob_parm.rho * prob_parm.eint;
  s_ext[UEDEN] = prob_parm.rho *
                 (prob_parm.eint + 0.5 * (prob_parm.vx_in * prob_parm.vx_in +
                                          prob_parm.vy_in * prob_parm.vy_in));
  s_ext[UTEMP] = prob_parm.T;
  for (int n = 0; n < NUM_SPECIES; n++)
    s_ext[UFS + n] = prob_parm.rho * prob_parm.massfrac[n];
    }
}

[RSuryaNarayan]

The compilation runs fine except for the part on the s_ext[UMX] = -s_int[UMX] because they are of different types

[hsitaram]

 s_ext[UEINT] = prob_parm.rho * prob_parm.eint;
s_ext[UEDEN] = prob_parm.rho *
                 (prob_parm.eint + 0.5 * (prob_parm.vx_in * prob_parm.vx_in +
                                          prob_parm.vy_in * prob_parm.vy_in));
  for (int n = 0; n < NUM_SPECIES; n++)
      s_ext[UFS + n] = prob_parm.rho * prob_parm.massfrac[n];```

These lines for your isothermal wall may not be correct.
Remember you would still want to impose zerogradient for species and pressure. So i would s_int for each of these variables too.

[RSuryaNarayan]

@hsitaram Thanks for that. so finally would that work out to: s_ext[UEINT] = s_int[UEINT]?

[hsitaram]

s_ext[URHO,UEDEN,UFS+sp] should be set to s_int. UEINT is not really used in the conserved to primitive conversion. s_ext[UEINT]=s_int[UEINT] should work ok. Not sure what the compile issue is: bcnormal takes s_ext and s_int as regular Real arrays of size NVAR. you can drop the const amrex::Real* for s_int in the argument list. just do amrex::Real s_int[NVAR]

[RSuryaNarayan]

Yes that dropping of the * worked.

[hsitaram]

that is correct.

[emotheau]

In s_ext[UMX] = -s_int[UMX], does s_int is taking values from the stencil inside the domain? Or is it just taking the last interior cell of the domain and making a copy to all the exterior ghost-cells?

[RSuryaNarayan]

I have a similar question. Intuitively one would just set the wall value to zero for a no-slip. So I am guessing that the s_ext represents the ghost cell and the s_int represents the last column/row of cells?

[hsitaram]

@emotheau - the second, its taking value from the last interior cell and doing the "reflect_odd" for the exterior ghost-cells.

[RSuryaNarayan]

@hsitaram just one last clarification from your side (thanks for being so patient so far!) The ghost cells here in the pic are represented by the blue dots and the interior by the red dots. The physical boundary is also indicated. So does reflect odd simply set the blue dots to negative of the red dots?

[drummerdoc]

Right, it's to reflect the internal state across the physical boundary and across zero...so the row below the orange circles would be negated and used to fill the row above the blue circles if grow>1, etc.

[RSuryaNarayan]

Thanks a lot @hsitaram @marchdf @drummerdoc !! that clarified a lot of stuff we theorized about Pele so far!

[emotheau]

@emotheau - the second, its taking value from the last interior cell and doing the "reflect_odd" for the exterior ghost-cells.

So if you have more than 1 exterior ghost-cell (for PPM for example), it is not actually doing properly a mirror relation to enforce the BC, right?

[RSuryaNarayan]

So @drummerdoc this reflect odd procedure will be implemented across the blue-orange cells in addition to the cells tick-marked ryt?

[drummerdoc]

That is what the filcc routine in AMReX does when it is called to fill cells at a reflect odd boundary.

[RSuryaNarayan]

Just a follow up - the bcnormal function I finally used was like so:

AMREX_GPU_DEVICE
AMREX_FORCE_INLINE
void
bcnormal(
  const amrex::Real* /*x[AMREX_SPACEDIM]*/,
  const amrex::Real s_int[NVAR],
  amrex::Real s_ext[NVAR],
  const int idir,
  const int sgn,
  const amrex::Real /*time*/,
  amrex::GeometryData const& /*geomdata*/,
  ProbParmDevice const& prob_parm)
{
  //considering both X and Y are set to "Hard"
  //impose the BC on the top walls:
  if ((idir==1 && sgn==1) || (idir==1 && sgn==-1)) 
  {
    //set the no-slip BC    
    s_ext[UMX] = -s_int[UMX];
    s_ext[UMY] = -s_int[UMY];
    s_ext[UMZ] = -s_int[UMZ];
    //set the isothermal wall
    s_ext[UTEMP] = 300.0; //for the isothermal wall at 300 K
    //for other variables set zero grad
    s_ext[URHO] = s_int[URHO];
    s_ext[UEINT] = s_int[UEINT];
    s_ext[UEDEN] = s_int[UEDEN];
    for (int n = 0; n < NUM_SPECIES; n++)
      s_ext[UFS + n] = s_int[UFS+n];
  }
  else
  {
  s_ext[URHO] = prob_parm.rho;
  s_ext[UMX] = prob_parm.rho * prob_parm.vx_in;
  s_ext[UMY] = prob_parm.rho * prob_parm.vy_in;
  s_ext[UMZ] = 0.0;
  s_ext[UEINT] = prob_parm.rho * prob_parm.eint;
  s_ext[UEDEN] = prob_parm.rho *
                 (prob_parm.eint + 0.5 * (prob_parm.vx_in * prob_parm.vx_in +
                                          prob_parm.vy_in * prob_parm.vy_in));
  s_ext[UTEMP] = prob_parm.T;
  for (int n = 0; n < NUM_SPECIES; n++)
    s_ext[UFS + n] = prob_parm.rho * prob_parm.massfrac[n];
    }
}

This seemed to produce the exact same result from the EB = TRUE case! I am not really sure whether the isothermal wall should show up in the plot considering we are imposing the hard temperature BC on the ghost cell (which don't get plotted) using s_ext[UTEMP]=300. So what I see in the plot at the wall are cells one row below the ghost cells whose temperature is allowed to vary?

[drummerdoc]

Right. Specifically, Dirichlet boundary data is actually enforced on the bounding WALL. It can be confusing, but the AMReX codes use a cell-centered container to pass this information, both for regular grids and for EB. In both cases the container has the correct shape (one slot for each value). Internal to the code, at the level of the stencils applied for advection and diffusion, the code knows that the Dirichlet data is to be applied on the face.

[RSuryaNarayan]

@drummerdoc thanks for those prompt replies. Made it even more clear. I still have one question lingering - what exactly does the Hard keyword trigger in the inputs file and how does it influence the way bcnormal parses arguments? Because it seems to me we aren't using that part anywhere in the bcnormal function (It seems to me I could also set NoSlipWall and then define bcnormal the way I did for the Hard case).

[drummerdoc]

It's related to BCs specifically for compressible flows, and refers to whether the external boundary condition is held fixed, independent of the internal state. If so, perturbations traveling to that boundary from inside the domain will typically reflect as if they hit a hard wall. Alternatively, there is an entire "NSCBC" infrastructure to dial in a "soft" condition that absorbs the outgoing waves with varying degrees of reflection. I hesitate to bring that up because the next question will be how to control the NSCBC response, and I afraid I'm not completely up on the status of the porting of that infrastructure to the current PeleC framework (which underwent considerable development associated with porting to GPUs).

[RSuryaNarayan]

Okay but I am specifically coding in the response of the wall (using the state vectors) in bcnormal. So what I am unable to understand is the role the keywords in the inputs file do given I have to manually implement them in the bcnormal function. Maybe I am wrong in assuming that ALL of the keywords in the inputs file need bcnormal definitions?

[drummerdoc]

Sorry :) for that I have to appeal to the folks that implemented the most recent interface.

[RSuryaNarayan]

No worries. I think I can find an answer to that in the PeleC.cpp file that parses arguments from the inputs to the the main program using ParmParse. Also @drummerdoc if I would want to impose a turbulent inflow BC for the inlet in the channel, would bcnormal still work?

[drummerdoc]

Sure, subject to similar characteristics issues. You would enforce the Dirichlet values on the upstream side of the inflow face. The algorithm extrapolates the internal state to back to the other side of that face and a Riemann solver resolves the wave structure to determine the ultimate state on the face. If the inflow is subsonic, you have the issue that acoustic signals can propagate back against the flow. Physically, a sound wave will approach the inlet face and the pressure will not be identical on both sides of the face. The Riemann solver would compute a reflection of that acoustic wave back into the domain, and your inlet flow would contain the flow plus the acoustic reflection. Again, there is the possibility of "softening" this with an NSCBC approach that has parameters for how absorbing the boundary is.

Of course, you could just ignore all this and sit back and admire the pretty acoustic reflections that result. they'll likely be tiny waves rattling around your domain but they will be there.

[emotheau]

Back in the days, the Hard keyword was put to trigger the bcnornal routine that will fill the ghost-cells at the domain BCs. Unless you use NSCBC ( https://ccse.lbl.gov/people/motheau/Manuscripts_website/2017_AIAA_CFD_Motheau.pdf ), it's just leading to unphysical problems. Regarding your questions above, my understanding is that you are imposing incorrectly a wall BC because the ghost-cells are not properly filled. I went to exactly the same problems many years ago...

[RSuryaNarayan]

Thanks for that reference. @emotheau would that mean the other keywords don't trigger the bcnormal function? Also the wall BC for the momentum has simply been implemented both using NoSlipWall and Hard+bcnormal (both seem to fetch the same results). For the isothermal part I just set the ghost cell to the required wall temperature. Could you please explain on what part would be wrong? thanks

[hsitaram]

@emotheau. Ppm will just drop to first order at the wall when the outer ghost cells are filled with the same value as the inner ghost cells. I don't think @Rsuryanarayan 's way is incorrect. Can you also set

S_ext[temp] = 2.0*300 - S_int[temp]

This will be more accurate

On Wed, Mar 31, 2021, 00:57 RSuryaNarayan @.***> wrote:

Thanks for that reference. @emotheau https://github.com/emotheau would that mean the other keywords don't trigger the bcnormal function? Also the wall BC for the momentum has simply been implemented both using NoSlipWall and Hard+bcnormal (both seem to fetch the same results). For the isothermal part I just set the ghost cell to the required wall temperature. Could you please explain on what part would be wrong? thanks

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/AMReX-Combustion/PeleC/issues/263#issuecomment-810823607, or unsubscribe https://github.com/notifications/unsubscribe-auth/ABYOQM5PMT3GYIW57EHL36TTGLBWXANCNFSM42BXXYFA .

[RSuryaNarayan]

@hsitaram we did try out the cell averaging technique that you just mentioned (600- s_int[UTEMP]) the results were exactly the same. Also these conformed to the EB =TRUE case (where the box covered the domain) However we are using MOL right now.

[RSuryaNarayan]

Here's the MOL result for 5000 timesteps

[hsitaram]

Looks ok to me. I see a layer of higher temperature in the boundary layer, which is typical of compressible boundary layers with isothermal walls. Have you run with a more refined mesh?

On Wed, Mar 31, 2021, 08:18 RSuryaNarayan @.***> wrote:

Here's the MOL result for 5000 timesteps [image: Screenshot from 2021-03-31 19-46-07] https://user-images.githubusercontent.com/52670594/113159170-05fb4f80-925a-11eb-8032-6fe67a80b6fe.png

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/AMReX-Combustion/PeleC/issues/263#issuecomment-811104880, or unsubscribe https://github.com/notifications/unsubscribe-auth/ABYOQM7CNFYRV5QMS4MI7Z3TGMVMZANCNFSM42BXXYFA .

[RSuryaNarayan]

@hsitaram Thanks for that heads up with the boundary layer. I am trying to validate these with the NASA dataset for supersonic channel flows but I am yet to figure out a way to vizualize the dataset. Regarding running with a finer mesh, not yet. I'll have to connect to my lab cluster to get more computing power. Currently running on my laptop. Will post the result when I get :-). Also I am trying out pelec.do_ppm = 1 to see if we can get sharper shocks

[RSuryaNarayan]

Doubled the mesh used previously from (192, 64, 16) to (384, 128, 32) for a (18,6,1.5) domain. Here are some results: . mach number: You may find the tarball of the plot files along with the movies for the various variables here. I think not much has changed after doubling the mesh.

[RSuryaNarayan]

I have a question on space-dependent BCs using bcnormal. Consider a turbulent inflow boundary condition in this case, where some perturbation is required on the inlet velocity profile that depends on both space and time. For instance, the inlet velocity function could be: u_net = u_avg + u_pert (y,t) Where the u_pert represents the perturbation component that varies with the y-direction and time. There has been a nice provision to include time in the bcnormal function using Real time, however, I am not sure of the space component, is it represented by the first argument? Am I right in assuming x[0], x[1] and x[1] represent the X, Y, and Z coordinate? In that case the turbulent inflow for a simple sinusoidal perturbation:

bcnormal(
  const amrex::Real x[AMREX_SPACEDIM],
  const amrex::Real s_int[NVAR],
  amrex::Real s_ext[NVAR],
  const int idir,
  const int sgn,
  const amrex::Real time,
  amrex::GeometryData const& /*geomdata*/,
  ProbParmDevice const& prob_parm)
{
 amrex::Real u_turb, u_mean{70000};
         if ((idir==0 && sgn==1) )
          {
               u_turb = u_mean + 0.05 *u_mean*sin(2*pi*x[1]/1000)*sin(2*pi*time/1000); // the 1000s in the denos are simply some placeholders for some characteristic wave-number and time scale which we can adjust later on
               s_ext[UMX] = u_turb-s_int[UMX];
           }
}

Any corrections to the above function? thanks in advance

[jrood-nrel]

I just want to make a comment that I have discussions enabled for this project, so feel free to use that for more discussion related inquiries like this issue.