Subsonic Inflow Boundary Condition Problem

[mafiami]

Hi! I'm working on a Ma > 0.3 subsonic jet simulation,but the numerical result of applying the same inflow boundary condition in EB-ConvergingNozzle is not so satisfying, the pressure at the boundary decreases as time advance, so I want to use Characteristic boundary condition in my simulation.
I tried to apply the characteristic boundary condition but it diverged, here's the code I used.

               amrex::Real pb, rhob, ub, vb, wb;
               amrex::Real rm = -uin - 5.0 * a; 
               amrex::Real cb = -0.2 * rm;
               amrex::Real c02 = a*a + 0.2 * (uin * uin + vin * vin + win * win);
               amrex::Real Tb = prob_parm.T0 * (cb * cb / c02);
               pb = prob_parm.pa * prob_parm.scalepa * std::pow((Tb/prob_parm.T0),3.5);
               eos.PYT2R(pb, massfrac, Tb, rhob);
               amrex::Real cp;
               eos.TY2Cp(Tb, massfrac, cp);
               ub = std::sqrt(2.0 * cp * (prob_parm.T0 - Tb));
               const amrex::Real* prob_lo = geomdata.ProbLo();
               const amrex::Real* dx = geomdata.CellSize();
               const amrex::Real xd = prob_lo[0] + 0.5 * dx[0];
               const amrex::Real xb = prob_lo[0];
               vb = 0.0;
               wb = 0.0;
 
               rho = (rhob - rhoin) / (xb - xd) * (x[0] - xd) + rhoin;
               u = (ub - uin) / (xb - xd) * (x[0] - xd) + uin;
               v = (vb - vin) / (xb - xd) * (x[0] - xd) + vin;
               w = (wb - win) / (xb - xd) * (x[0] - xd) + win;
               p = (pb - pin) / (xb - xd) * (x[0] - xd) + pin;           
               eos.RYP2E(rho, massfrac, p, eint);
               eos.REY2T(rho, eint, massfrac, T);

In EB-ConvergingNozzle case, it shows this kind of inlet boundary condition is invalid. Does it mean this kind of boundary condition is not practical in PeleC?
Anything else I can do to make it match the nozzle static pressure?

[baperry2]

I created the EB-ConvergingNozzle case. It was a while ago and I don't remember why I didn't implement the characteristic extrapolation inflow - perhaps it was just that the fixed pressure inlet worked well enough for that case so I didn't bother going further. Looking at it again now, I'm not seeing a physical reason why this wouldn't work at least with GammaLaw EoS, which you appear to be assuming based on your implementation.

Specifying the correct boundary conditions depends a lot on your case and what you're interested in. If your case is a jet into an open domain, you may have success using the inlet BC that takes pressure from the domain interior coupled with the characteristic extrapolation outlets that specify outlet pressure (the default setup for EB-ConvergingNozzle).

In case it's useful to you, I also have an implementation that uses fixed stagnation conditions for the inlet and pulls the Mach number from the domain interior. The code for that is:

      // stuff below assumes ideal gas
      if (eos.identifier() != "Fuego") {
        amrex::Abort("Must use Fuego EOS for this bc");
      }

      amrex::Real massfrac_int[NUM_SPECIES] = {0.0};
      amrex::Real rho_inv = 1.0/s_int[URHO];
      for (int n = 0; n < NUM_SPECIES; n++) {
        massfrac_int[n] = s_int[UFS + n] * rho_inv;
      }
      amrex::Real cs_int;
      eos.RTY2Cs(s_int[URHO], s_int[UTEMP], massfrac_int, cs_int);
      amrex::Real M_int = s_int[UMX] * rho_inv / cs_int;
      amrex::Real rho0, cs0, h0, entrop0;
      CKSBMS(prob_parm.p0, prob_parm.T0, massfrac_ext, entrop0);
      eos.PYT2R(prob_parm.p0, massfrac_int, prob_parm.T0, rho0);
      eos.RTY2H(rho0, prob_parm.T0, massfrac_ext, h0);
      eos.RTY2Cs(rho0, prob_parm.T0, massfrac_ext, cs0);

      amrex::Real cs_ext = cs0;
      amrex::Real M_ext = M_int;
      rho_ext = rho0;
      u_ext = M_ext * cs_ext;
      amrex::Real dT = 1e6;
      constexpr amrex::Real tol = 1e-8;
      constexpr int max_steps = 1000;
      int step = 0;

      while (dT > tol && step < max_steps) {
        step += 1;

        amrex::Real h_ext = h0 - 0.5 * u_ext*u_ext;
        amrex::Real Told = T_ext;
        eos.RHY2T(rho_ext, h_ext, massfrac_ext, T_ext); // adiabatic T
        dT = std::abs(T_ext - Told);
        eos.RTY2P(rho_ext, T_ext, massfrac_ext, p_ext);
        amrex::Real entrop_ext;
        CKSBMS(p_ext, T_ext, massfrac_ext, entrop_ext);

        amrex::Real cp, cv, Rinv;
        eos.RTY2Cp(rho_ext, T_ext, massfrac_ext, cp);
        eos.RTY2Cv(rho_ext, T_ext, massfrac_ext, cv);
        Rinv = 1.0 /(cp - cv);

        amrex::Real dens_rat = std::exp((entrop_ext - entrop0) * Rinv);
        rho_ext *= dens_rat; // density consistent with T_ext and isentropic
        eos.RTY2Cs(rho_ext, T_ext, massfrac_ext, cs_ext);
        u_ext = M_ext * cs_ext;
      }

      if (step >= max_steps) {
        amrex::Abort("Iteration to set inflow condition based on Mach did not converge");
      }

      eos.RTY2P(rho_ext, T_ext, massfrac_ext, p_ext);
      eos.RTY2E(rho_ext, T_ext, massfrac_ext, eint_ext);