Modeling Heavy Gas (Radon) Settling Under Gravity Using NavierStokes Physics (or others) Module

[Traiwit]

Hi guys,

I'm working on modeling radon gas transport in underground mine stopes, where the heavy radon gas (density ~9.73 kg/m³) settles under gravity in air (density ~1.225 kg/m³). The simulation needs to handle two scenarios:

Static conditions: No ventilation, radon settles purely due to gravity (buoyancy-driven flow)
Dynamic conditions: With mine ventilation, combining forced convection with gravitational settling

The physical setup involves two connected chambers where radon concentration starts high in one chamber and should settle to the bottom while diffusing to the adjacent chamber over time.

lets focus on Static conditions first, where there is no flow

[Physics]
  [NavierStokes]
    [Flow]
      [flow]
        compressibility = 'incompressible'
        porous_medium_treatment = false
        density = 'rho'  # Should be variable based on radon concentration
        dynamic_viscosity = 'mu'
        initial_velocity = '0 0 0'
        initial_pressure = 0 
        wall_boundaries = 'wall'
        momentum_wall_types = 'noslip'
        
        # For future ventilation cases
        # inlet_boundaries = 'inlet'
        # momentum_inlet_types = 'fixed-velocity'
        # outlet_boundaries = 'outlet'
        # momentum_outlet_types = 'fixed-pressure'
        
        block = '0 1'
      []
    []
  []
  
  [ScalarTransport]
    [Radon]
      coupled_flow_physics = 'flow'
      passive_scalar_names = 'radon_concentration'
      passive_scalar_diffusivity = 2.3e-6  # Radon diffusivity in air
      passive_scalar_coupled_source = 'radon_concentration'
      passive_scalar_coupled_source_coeff = -2.1e-6  # Rn-222 decay
      initial_scalar_variables = '0'
      block = '0 1'
    []
  []
[]

Expected Behavior

Initially: High radon concentration at the top of chamber 1
Short term: Radon settles downward due to higher density, creating convection currents
Long term: Radon accumulates at the bottom of both chambers with concentration gradient

i'm not sure if MOOSE's NS module can model this kinda setup, is there any suggestion?

do I need to do variable density based on concentration within the element & mass/volume fraction, or that's way too complicated?

I tried to turn gravity on [Flow], it's much more difficult to converge and concentration doesn't seem to move down as I expect

Thank you,
T

[GiudGiud]

Hello

Short term: Radon settles downward due to higher density, creating convection currents
Long term: Radon accumulates at the bottom of both chambers with concentration gradient

Your input wont do that. Because you dont have any buoyancy term from the concentration of Radon.
You need to define a custom rho that depends on the local fraction of the radon scalar.

Technically this will be abusing the ScalarTransport Physics which is for passive scalar that do not influence the fluid properties.
It might work anyway and would be cheaper to run

The appropriate Physics is the TwoPhaseMixture physics. This can accomodate change in the fluid properties from the presence of the gas.

[GiudGiud]

it kinda works but it extremely difficult to converge:

the case is ill-posed. You dont have a pressure outlet and you use incompressible flow. So the pressure can be shifted by any amount and still be a valid solution.

You need to pass pressure_pin parameters to the Flow Physics. I find the velocity profile to be odd after that but the scalar concentration believable. The runtime seems fine to me on an M4

        pin_pressure = true
        pinned_pressure_type = 'average'
        pinned_pressure_value = 0

[Traiwit]

Hi @GiudGiud,

Thanks for testing it

I tried that and I still got a few of these Solve Did NOT Converge! at the start, not sure if this an expected behaviour in NS, coming from PorousFlow seeing these make me feel uncomfortable :D

and i agree the flow path is very weird, scalar seems to move along the walls rather than the body of the geometry

Time Step 17, time = 10.23, dt = 2.56
 0 Nonlinear |R| = 7.668309e+02
 1 Nonlinear |R| = 2.024912e+03
 2 Nonlinear |R| = 4.748294e+02
    Linear solve did not converge due to DIVERGED_ITS iterations 100
 3 Nonlinear |R| = 1.337957e+04
    Linear solve did not converge due to DIVERGED_ITS iterations 100
 4 Nonlinear |R| = 8.893192e+03
    Linear solve did not converge due to DIVERGED_ITS iterations 100
 5 Nonlinear |R| = 4.716313e+05
    Linear solve did not converge due to DIVERGED_ITS iterations 100
 6 Nonlinear |R| = 5.258530e+07
    Linear solve did not converge due to DIVERGED_ITS iterations 100
 7 Nonlinear |R| = 3.603765e+11
    Linear solve did not converge due to DIVERGED_ITS iterations 100
 8 Nonlinear |R| = 6.655545e+14
  Nonlinear solve did not converge due to DIVERGED_DTOL iterations 8
 Solve Did NOT Converge!
  Finished Solving                                                                       [  5.35 s] [  287 MB]
Aborting as solve did not converge

Solve failed, cutting timestep.

Time Step 17, time = 8.95, dt = 1.28
 0 Nonlinear |R| = 7.593435e+02
 1 Nonlinear |R| = 3.756791e+02
 2 Nonlinear |R| = 6.369841e+01
 3 Nonlinear |R| = 1.156062e+00
 4 Nonlinear |R| = 8.997860e-03
 Solve Converged!

[GiudGiud]

not sure if this an expected behaviour in NS

yeah it's business as usual in NS. We cannot just take any time step in NS, the Newton method does not always find a solution even though the implicit time stepping should be stable for any size of time step.

The longer time steps at the end of the transient are a positive sign.

Also you should be using the weakly-compressible approximation since you have a variable rho.
And you should define the drho_dt functor with the dependence on the time derivative of the scalar. This might require a custom material though you might be able to place functor materials for time derivatives and parsed expressions. never tried

[Traiwit]

got it, tbh i'm not sure if what i'm doing is physically correct, but the results look okay (for now)

class below by Claude Code:

#include "ScalarDependentDensityTimeDerivativeFunctor.h"

registerMooseObject("moose_behApp", ScalarDependentDensityTimeDerivativeFunctor);

InputParameters
ScalarDependentDensityTimeDerivativeFunctor::validParams()
{
  InputParameters params = FunctorMaterial::validParams();
  params.addClassDescription("Functor material that computes the time derivative of density when density depends on a scalar variable");
  
  params.addRequiredParam<std::vector<MooseFunctorName>>("scalar_functor", "The scalar functor that density depends on");
  params.addRequiredParam<Real>("density_coefficient", "Coefficient relating scalar to density (drho/dscalar)");
  params.addParam<MooseFunctorName>("drho_dt_name", "drho_dt", "Name of the time derivative of density property");
  
  return params;
}

ScalarDependentDensityTimeDerivativeFunctor::ScalarDependentDensityTimeDerivativeFunctor(const InputParameters & parameters)
  : FunctorMaterial(parameters),
    _density_coeff(getParam<Real>("density_coefficient"))
{
  const auto & scalar_names = getParam<std::vector<MooseFunctorName>>("scalar_functor");
  const auto & drho_dt_name = getParam<MooseFunctorName>("drho_dt_name");
  
  addFunctorProperty<ADReal>(
    drho_dt_name,
    [this, scalar_names](const auto & r, const auto & t) -> ADReal
    [{](`url`)
      const auto & scalar_functor = getFunctor<ADReal>(scalar_names[0]);
      return _density_coeff * scalar_functor.dot(r, t);
    });
}

  # Time derivative of density for weakly-compressible formulation
  # Custom functor material that computes drho_dt = (1/500) * d(scalar)/dt
  [density_time_derivative]
    type = ScalarDependentDensityTimeDerivativeFunctor
    scalar_functor = 'scalar'
    density_coefficient = 0.002  # 1/500 coefficient from rho = rho_base + scalar/500
    drho_dt_name = drho_dt
  []

[GiudGiud]

[{](`url`)

in the functor looks weird and wrong

      const auto & scalar_functor = getFunctor<ADReal>(scalar_names[0]);

call this before and pass it to the lambda function instead of "scalar_names". Way more efficient to avoid string comparison on every evaluation

Otherwise looks believable but my job still safe for now haha