Skin friction model replacement

[nlaing613]

The current Karman-Nikuradse skin friction model does not account for compressibility effects (density and temperature variations in a turbulent compressible boundary layer). The consequence is that skin friction is currently over-predicted (depending on Mach at fixed Re, as much as 50%, but in the cases of interest, roughly 20%).

To fix this, the analytical, closed-form DeChant skin friction model will be implemented. This provides an empirical correlation between skin friction, Mach number, local to wall temperature ratio, and Reynolds number. The current model constants implement a recovery factor for air, and as tabulated, assumes γ = 1.4. The result is a three-dimensional interpolated table-lookup. For reacting cases, this interpolation may have to be replaced with an in-place root-finding for skin friction due to variable Prandtl number and γ in high temperatures and differing mass fractions.

[nlaing613]

In the recent commit, the boundary layer class now must receive a user-selected SkinFriction model and HeatFlux model (optional). These will be instantiated in the runfile as the tuple "wall_models=(SkinFriction, HeatFlux)". The user is permitted to add a wall temperature without specifying the heat flux model, as this is needed for the compressible skin friction model.

If no model is selected for SkinFriction, then the boundary layer will not be instatiated. This is roughly equivalent to the boolean that used to be in its place: "include_boundary_layer".

[nlaing613]

@TimothyEDawson

I'm trying to ensure the functionality of the debugging file I made first and I'm getting the error:
` Traceback:
File "/examples/hifire2/constant_area_duct.py", line 106, in
ss.advance_simulation(tFinal)

File "/src/stanshock/components/combustor.py", line 288, in advance_simulation
self.t, state_array, gamma_star, e0_star = self.time_integrator.advance(
dt=dt, ...<3 lines> ... e0_star=e0_star)

File "/src/stanshock/numerics/time_integration.py", line 329, in advance
_, state_array_stage, gamma_star_temp, e0_star_temp = advance(
dt=cdt * dt,...<3 lines>...e0_star=e0_star_stage)

File "/src/stanshock/numerics/time_integration.py", line 138, in advance
dydt = self.rhs.source(time=t, ...<2 lines>...e0_star=e0_star_local )

File "/src/stanshock/system/base.py", line 204, in precompute_for_source
face_states = self.face_extrapolator(state)

File "/src/stanshock/numerics/face_extrapolation.py", line 476, in call
face_states_array: Array = weno5(
~~~~~^
r=state.density,
^^^^^^^^^^^^^^^^
...<5 lines>...
n_scalars_rho_sum=self.n_scalars_rho_sum,
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
)

ZeroDivisionError: division by zero`

I believe this must occur because of the double flux being disabled (which I selected in the constant_area_duct.py), but I am a little unsure. I will continue to poke around, but if you have any ideas, it'd be a big help. Thanks!

[TimothyEDawson]

@nlaing613 You're right, that's a bug I've encountered before as well, and you're correct that it comes from the approach used to disable the double flux. I simply added the following after line 373 in face_extrapolation.py as a maybe-temporary workaround:

rLR = max(rLR, 1e-30)

[nlaing613]

@TimothyEDawson One more thing, I think I've caught a problem with the strategy used to ensure WENO only uses the interior points in its stencil when double flux is disabled. This occurs in InviscidFlux.source_implementation(). The ghost cells are assigned (as you had mentioned) very large values in "DeactivateWenoCells(GhostCellPrimitives). This eventually gets passed to the face_states via boundary_conditions.update_face_states(...).

When it comes time to run the Riemann solver, (line 316) in InvsicidFlux.source_implementation(), the call to gamma = self.physics.get_gamma(face_states) fails because the resulting temperature values inputted to get_cp() include values that are far out of bounds for the table (I get a minimum temperature in the ghost cells of T_min = 5.39e-09).

A temporary fix I'll try is simply clipping the resulting temperature to Tmin and Tmax from the table, but I thought I'd just let you know what I found.

 File "/StanShock/src/stanshock/physics/thermotable.py", line 72, in get_cp_compiled
    raise ValueError(msg)
ValueError: Temperature out of bounds: <object type:float64> not in range [<object type:float64>, <object type:float64>]

Turns out I also have to do this in ViscousFlux, and get_temperature, too. I must be doing something wrong, so any ideas would be a big help. I did implement what you recommended, to be sure, it should be after the for-loop in "face_extrapolation.py"?

rLR = 0.0
for kSc in range(n_scalars_rho_sum):
    rLR += U[mn + kSc]
rLR = max(rLR, 1e-30)

After implementing the temperature clipping in ViscousFlux, I'm now also getting the issue of:

File "/src/stanshock/system/base.py", line 232, in source
    self.precompute_for_source(time, state_array_local, gamma_star, e0_star)

  File "/src/stanshock/system/base.py", line 189, in precompute_for_source
    state = self.physics.conservative_to_primitive(
        state_array_local, gamma_star, e0_star
    )
  File "/src/stanshock/physics/fluid_base.py", line 332, in conservative_to_primitive
    state.temperature = self.get_temperature(state)

  File "/src/stanshock/physics/thermotable.py", line 216, in get_temperature
    state = self.set_state(state)
  File "/src/stanshock/physics/cantera_interface.py", line 48, in set_state
    self.sol.UVY = (

  File "envs/shock/lib/python3.13/site-packages/cantera/composite.py", line 704, in __setattr__
    super().__setattr__(name, value)

  File "envs/shock/lib/python3.13/site-packages/cantera/composite.py", line 1383, in setter
    setattr(self._phase, name, (A[index], B[index], C[index]))

  File "thermo.pyx", line 1481, in cantera.thermo.ThermoPhase.UVY.__set__
cantera._utils.CanteraError: 
CanteraError thrown by Phase::setDensity:
density must be positive. density = 0

[TimothyEDawson]

@nlaing613 Hmmmm, so the part I don't understand is why the face states are invalid. What
should happen is that the face state at the boundary is simply extrapolated from the interior points, which should be valid.

Your code matches my fix, yes. I wonder if there was anything else I had done to workaround such problems, or if maybe there's something else going on in your specific case.

In general I am just not very happy with the DeactivateWenoCells as currently implemented, and would like to implement a better way to handle that in the future.

[TimothyEDawson]

@nlaing613 I think there may be a bug in the viscous flux, possibly related to this being a single-species simulation as I've only tested it with multiple species which include mass diffusion. For now, maybe try with include_diffusion=False and see how far you get with that. Your results may be a bit off since it won't include the viscous momentum flux and viscous heating which are also handled within ViscousFlux.

Meanwhile I'll try and make sure the constant area duct case works with viscous flux, I may open a separate pull request for that.

[TimothyEDawson]

@nlaing613 I've made some updates and simplifications to what you had. Some of the simplifications include turning a lot of optional values in the WallState dataclass into mandatory values (since currently they are computed regardless), and turning a lot of Array | float values into simply Array values.

I noticed that there are some changes from the original mode of operation which we may want to consider reverting/reimplementing. Namely:

• Several properties are now always computed even if they won't be required for the given model(s) (e.g. cp, k, mu). Previously those were only computed if a heat flux model required it, though I do see the CompressibleReactingSkinFriction requires the Prandtl number so it's no longer that straightforward.
• The original get_nusselt_number routine applied different equations for low and high Reynold's numbers, and this seems to have been removed entirely. It would be good to retain this functionality as an option somewhere in the wall models.

I also caught a bug in the boundary conditions for the constant_area_duct example case which comes down to the location="left"/"right" needing to be passed in at initialization of the SpecifiedFace objects. It defaults to left, so what was supposed to be the right face was applying p2 to the inflow (left) boundary pressure instead! Woops! I almost missed that, and I'm not liking the duplicated information. I'll have to think about whether there's a good way to avoid this in the future.