Wall boundary conditions (e.g. NoSlipWallIsothermal) do lead to wrong behaviour using the soret effect.

[stfsschneider]

When the soret effect is turned on (simple, fuego), we can set up a test case, where species are "produced" at the wall due to zero gradient boundary condition. For soret, a zero gradient boundary conidition for the species is actually not correct, since the fluxes due to thermal diffusion have to be balanced by fluxes in the species profiles, leading to non zero gradient species profiles at the walls. I created a testcase (pipe with two walls on the top and the bottom with perodic bcs on the left and right) and set a temperature gradient from 300K to 700K from the bottom of the domain to the top (top wall temperature = 700K, bottom wall temperatur = 300K) and initialized the species profiles with a mixture of air and H2 (distributed equally over the domain). Over time the mass fraction of H2 increases and the mass fluxes of O2 and N2 decrease, which is physically not correct.
If you have any questions, I can also share the test case and we can discuss on this!

Thanks

[baperry2]

I think I follow what you're saying - basically at isothermal walls the species BC should be no flux rather than zero gradient, which is usually equivalent but ends up being subtly different when Soret effects are active. Depending on approach, implementation of a solution could be a bit tricky as the diffusion solve is done implicitly using AMReX's linear solvers, and setting complex boundary conditions for these solves is often not straightforward. It might easiest to just turn off the Soret fluxes at the walls, which should conserve mass but may also incur some error. @ThomasHowarth, do you have any thoughts on this from when you did the Soret implementation?

@stfsschneider I'd be happy to discuss proposed fixes and guide you through the PR process if you are interested in implementing one. Sharing the test case would be a good first step.

[stfsschneider]

Thank you for you quick and helpful response!

Yes exactly, that sums it up quiet nicely! In OpenFOAM we were also facing a similar issue and solved it by rearranging the following equation to Yi (with a Finite volume discretization of the gradYi at the boundary/interface):
YV = - rho * Dmixi * gradYi - rho * Dmixi / W * Yi * gradMMean - Dtherm / T * gradT= 0

Would something similar be possible (accessing values of the field in the boundary conditions)? Otherwise a RobinBC which seems to be implemented in AmRex might be useful.

I would also be very happy to implement the changes, at the moment I still try to figure out at which point i would have to make changes to the code (as you also pointed out this is not straightforward)

[ThomasHowarth]

I didn't consider this when implementing Soret effects. Preferential diffusion effects near hot walls have been seen in literature, but obviously we don't want leakage through the boundary. As you've both said, the two options here would be:

1. ('Lazy' option) Neglect the Soret term at the boundary interface. This would result in zero species flux at the boundary, but isn't technically correct.
2. Impose a gradient in the light species that causes a flux that acts opposite and equal to the Soret flux. I don't think this would require Robin boundary conditions, but would require a non-zero, and variable, Neumann boundary condition. An additional complication here is that the boundary conditions for the two species would be coupled via the gradWbar term.

We also need to think about how this applies for isothermal EB. I'm even less familiar with how EB boundaries are imposed than non-EB ones, but the same problem will apply there.

[stfsschneider]

Yes, you are right, a Robin BC would only be needed if we want a specified mass flux (unequal to zero). Is there any similar BC you know of in Pele, that uses a variable Neumann BC?

[stfsschneider]

Can you allow me to push to a new branch? I can then push the test case to a new branch so you can also have a look into it

[ThomasHowarth]

AMReX offers support for inhomogeneous Neumann conditions (see Boundary conditions), but as far as I can tell this particular boundary type has not been implemented anywhere in Pele. I imagine we would have to change the section where it sets the DiffusionLinOpBC around this function and associated calls for the light species, and we'd then compute what the RHS needs to be and feed it a MultiFAB with these values and the backend linear solvers will do all the work for us. @baperry2 / @marchdf / @esclapez or others can maybe confirm this before I go breaking things? Also, I think the simplest way to see a case would be to have your own repository and copy a link to the case here.

[stfsschneider]

Here is a link to the testcase: https://github.com/stfsschneider/PeleLMeX/tree/development/Exec/RegTests/TestZeroFlux

[baperry2]

@stfsschneider Thanks for linking to the test case. We've been trying to keep the main repo cleaner by having folks create branches in their own forks as you have now done.

I think the approach Thomas proposes should generally work, although I'm not too familiar with setting BCs for the linear solvers either. The radiation solver in Pele uses variable Robin conditions so you may be able to get an idea of how to set things up for the variable Neumann conditions from that: https://github.com/AMReX-Combustion/PelePhysics/blob/19ccf8557a69d18cec8354ddde4af8e97d384cef/Source/Radiation/POneMulti.H#L117

@ThomasHowarth agreed there would be the same issue for isothermal EB. Although if I'm following the Soret code correctly, it looks like Soret fluxes aren't computed through the EB (correct for adiabatic EB), so for isothermal EB effectively the current status would be the "lazy" solution.

[drummerdoc]

@ThomasHowarth Yes, I believe you are correct in your suggestion. Somewhat nontrivial to implement though. Would welcome attempts to do so by anyone on the line here :)

[ThomasHowarth]

Feel free to check my working, but we do need a Robin boundary condition, not Neumann as I initially thought, and it needs to be applied to every species:

[stfsschneider]

Yes, I checked your derivation and get the exact same equation. One could also use a term containing molar masses instead of the chi_i and the temperature on the left hand side, but this is probably easier (since less dependencies on other variables/fields).

It, however, is only valid, if no correction velocity is assumed for the species. Otherwise, the correction velocity would also be required to take into account in the derivation.

[ThomasHowarth]

Since the diffusive flux is zero for all species at the boundary by the way we have constructed this, then the sum of the fluxes is zero, and there will be no correction velocity

[stfsschneider]

Yes, you are correct, I didn't think about that.

[baperry2]

One thing to consider - I haven't fully this out but I think it will change things - is that to discretely enforce species conservation, the fluxes must balance as computed numerically. This is relevant because the Wbar and Soret terms are computed on a lagged basis relative to the gradYk fluxes that are computed implicitly on each SDC iteration. Taking into account SDC iteration indices may result in a constraint that looks like a variable inhomogeneous Neumann BC rather than a Robin BC. It's also worth noting that a zero gradient BC is applied in computing the wbar terms and that also wouldn't be physically accurate. That may be harder to correct but also less impactful.

[ThomasHowarth]

Do you think I could therefore just impose an inhomogenous Neumann condition with the lagged Wbar and Soret fluxes? i.e. something like
$$\rho D_{mix} \nabla Y_{k}^{(j)} = W_{flux}^{(j-1)} + S_{flux}^{(j-1)}$$ for SDC iteration j
This, in theory, is much simpler than coupling everything together.

[drummerdoc]

Losing some of the thread here, but Thomas, your suggestion seems consistent with the overall design, provided you still use a correction velocity adjustment after the solve giving the flux listed above, since that alone doesn't guarantee the fluxes summed over the species is zero. right?

[baperry2]

Yeah that's pretty much what I was thinking. And Wflux would be zero at the boundary unless an improved BC is applied when computing gradW.

As Thomas mentioned above, no need to deal with the correction velocity because the species are individually forced to have no flux, so the total flux must be zero.

[esclapez]

I like the lagged version, it is a lot simpler and consistent with the way LMeX uses SDC iterations.
The relatively fast species diffusion solve, especially on GPU, relies on the fact that all the species are solved together in a single MLMG solve, which requires all the species to have the same BCs. So from an implementation point of view either all the species have to be treated with non-homogeneous Neumann on isothermal walls, or we keep the homogeneous Neumann but add a source term on the RHS of the solve with the added flux near the isoT walls for the light species.

[stfsschneider]

As I currently understand, for an isothermal wall, a dirichlet boundary condition is currently used for the species, in which the values on the faces are set as values from the centers. One would therefore have to change the boundary condition and catch the Soret effect and the isothermal wall in this function.
Moreover, one would need to make changes in the diffuse_scalar function and set the soretFluxes and wbarFluxes divided by rho and Dmix in the setLevelBC function. Do we need to make changes to any other methods of the DiffusionOp?

[ThomasHowarth]

So I think I've found where I need to make these changes, however, I've noticed another difficulty, in that the boundary conditions are being imposed for $\nabla (\rho Y_{k})$, rather than $\nabla Y_{k}$. Since we are creating gradients in the mass fractions, we will also be creating a density gradient - any suggestions for this?

[drummerdoc]

@ThomasHowarth I'm confused by this. It very definitely should be working with Grad(y), not Grad(rho.Y). We copy rhoY from the state and divide it by rho prior to working with it in the operator. Rho is explicitly multiplied back in when forming the rhs and the CC part of the operator.

@stfsschneider By copying out the cell-centered values "to the faces", that is NOT setting a Dirichlet value, that is enforcing a Neumann condition on mass fraction. Directly after that operation, the gradient is computed in the applyOp that will get a zero gradient identically, and thus a zero value for rho.D.Grad(Y).

[stfsschneider]

@ThomasHowarth Would a gradient in density be wrong? Since the composition at the wall and at the cell-center of the wall nearest cell differ, the thermochemical state and also the density at the wall do differ, or am I missing something? But as @drummerdoc pointed out, the way I understand the code is, that rhoY is divided by rho before the solve.

@drummerdoc Yes, I understand this, I maybe did not formulate it well enough. But what I meant is, aren't we using a LinOpBCType::Dirichlet since, a amrex::BCType::ext_dir is set for the species in case of IsothermalWalls and therefore enforcing a Neumann BC implicitly instead of explicitly using a Neumann BC? Would this also work (Neumann BC and zero for the second argument in the setLevelBC function levelbcdata)?
If this is the case one could change it to an inhomogeneous NeumannBC and use a amrec::BCType::foextrap for the species, but I am not sure, if I understood it correctly in all details.

[ThomasHowarth]

@drummerdoc you start second guessing yourself at every stage in this sort of thing...

Just opened a PR - hopefully what I've done makes sense. I'm getting zero flux on the walls in the (slightly modified) test case @stfsschneider provided with what I've implemented.