2 questions: Flow: Porous Medium + Open Space Setup AND Mixing-Length Turbulence Error

[Traiwit]

Check these boxes if you have followed the posting rules.

• Q&A Navier Stokes is the most appropriate category for my question
• I have consulted the posting Guidelines on the Discussions front page
• I have searched the Discussions forum and my question has not been asked before
• I have searched the MOOSE website and the documentation does not answer my question
• I have formatted my post following the posting guidelines (avoid screenshots if possible, triple back quotes before/after pasted text, etc)

Question

Hi guys,

I have a few questions

1. for Physics/NavierStokes/Flow

if I want to do the flow in open space as well as porous media in the same simulation, can I simply use porous_medium_treatment = true with porosity of porous media = 1 for open space and porosity < 1 for porous media ?

I tried to break this into 2 [flow] blocks, but MOOSE doesn't let me do that

[Physics]
  [NavierStokes]
    [Flow]
      [flow]
        compressibility = 'incompressible'
        
        porous_medium_treatment = true
        porosity = porosity

        density = 'rho'
        dynamic_viscosity = 'mu'

        initial_velocity = '0 0 0'
        initial_pressure = 5

        inlet_boundaries = 'left '
        momentum_inlet_types = 'fixed-velocity'
        momentum_inlet_function = '0 ${u_inlet} 0'

        wall_boundaries = 'active'
        momentum_wall_types = 'noslip'

        outlet_boundaries = 'right right_porous'
        momentum_outlet_types = 'fixed-pressure-zero-gradient fixed-pressure-zero-gradient'
        pressure_function = '${p_outlet} 5'

        # mass_advection_interpolation = 'average'
        # momentum_advection_interpolation = 'average'

        transient = false
         block = '12 2'
      []
    []

  []
[]

[FunctorMaterials]
  [porosity]
    type = GenericFunctorMaterial
    prop_names = 'porosity'
    prop_values = porosity_fx
  []
[]


[Functions]
    [./porosity_fx]
      type = ParsedFunction
      expression = if(z>=-45,0.2,1)
    [../]
[]

2. when I set turbulence_handling = 'mixing-length' it gives this error, not sure what is the cause of this:

*** ERROR ***
We caught a MetaPhysicL error in while performing element or face loops. This is potentially due to AD not having a sufficiently large derivative container size. To increase the AD container size, you can run configure in the MOOSE root directory with the '--with-derivative-size=<n>' option and then recompile. Other causes of MetaPhysicL logic errors include evaluating functions where they are not defined or differentiable like sqrt (which gets called for vector norm functions) or log with arguments <= 0

below is my physics block:

[Physics]
  [NavierStokes]
    [Flow]
      [flow]
        compressibility = 'incompressible'
        
        porous_medium_treatment = true
        porosity = porosity

        density = 'rho'
        dynamic_viscosity = 'mu'

        initial_velocity = '0 0 0'
        initial_pressure = 5

        inlet_boundaries = 'left '
        momentum_inlet_types = 'fixed-velocity'
        momentum_inlet_function = '0 ${u_inlet} 0'

        wall_boundaries = 'active'
        momentum_wall_types = 'noslip'

        outlet_boundaries = 'right right_porous'
        momentum_outlet_types = 'fixed-pressure-zero-gradient fixed-pressure-zero-gradient'
        pressure_function = '${p_outlet} 5'

        # mass_advection_interpolation = 'average'
        # momentum_advection_interpolation = 'average'

        transient = false
         block = '12 2'
      []
    []

  []
[]

[Physics/NavierStokes]
  [FluidHeatTransfer]
    [heat]
      coupled_turbulence_physics = 'turbulence'
      thermal_conductivity = 'k'
      specific_heat = 'cp'

      fluid_temperature_variable = 'T_fluid'
      initial_temperature = 30 # '${T_inlet}'
      
      energy_inlet_types = 'heatflux'
      energy_inlet_functors = '${fparse u_inlet * rho * cp * T_inlet}'

      energy_wall_types = 'heatflux'
      energy_wall_functors = '0'

      ambient_convection_alpha = 'h_cv'
      ambient_temperature = 'T_solid'

      # energy_advection_interpolation = 'average'

       block = '12 2'
    []
  []
[]

[Physics/NavierStokes]
  [ScalarTransport]
    [Radon]
      coupled_turbulence_physics = 'turbulence'
      passive_scalar_names = 'scalar'
      coupled_flow_physics = 'flow'


      passive_scalar_diffusivity = 0.05
      passive_scalar_source = 'radon_source'

      initial_scalar_variables = 'ini_scalar_fx'

      passive_scalar_inlet_types = 'fixed-value'
      passive_scalar_inlet_function = '0'

      # passive_scalar_advection_interpolation = 'average'

      block = '12 2'
    []
  []
[]

[Physics/NavierStokes]
[Turbulence]
  [turbulence]
    # turbulence_handling = 'mixing-length'

    coupled_flow_physics = flow
    fluid_heat_transfer_physics = heat
    scalar_transport_physics = Radon

    passive_scalar_schmidt_number = 1.0

    # # von_karman_const = ${von_karman_const}
    # mixing_length_delta = 1e9
    # mixing_length_walls = 'active'
    # mixing_length_aux_execute_on = 'initial'

    block = '12 2'
  []
[]
[]

Thank you!
Traiwit

Additional information

Mesh size and type:
Reynolds number: not sure
Discretization (finite element CG/DG, finite volume, etc):
Models (turbulence, porous media, etc): turbulence, porous media, flow, heat
Solver method (fully coupled, segregated, multiapps, etc): fully coupled
Base input you started from:

[GiudGiud]

Hello

1. yes. Porosity of 1 works fine and gets the exact same result as non-porous flow

2. we caught a MetaPhysicL error in while performing element or face loops. This is potentially due to AD not having a sufficiently large derivative container size. To increase the AD container size, you can run configure in the MOOSE root directory with the '--with-derivative-size=<n>' option and then recompile.

follow what the message says please. it s a known issue because of the gradient of the velocities on the faces in that model, which makes for a large stencil and a large AD container being used as a result

[Traiwit]

Hi @GiudGiud, thank you!

this's gonna sound really dumb, but what is the 'derivative container size' ? is this going to be different for each problem?

[GiudGiud]

it s an array that contains the derivatives, attached to each AD number.

So an ADReal is a Real + this array
the array stores the value and index of the derivatives of this number with regards to the nonlinear varibales

is this going to be different for each problem?

it will be different based on which kernels you use, and how far the stencils go in finite volume. But only the max size is set at configure time. A size of 140 should do it for mixing-length + 3D + nonlinear FV

[Traiwit]

Hi @GiudGiud,

--with-derivative-size=140 still give the same error, I also tried 1,500,000 + removed heat transfer from physics also still give the same error

when turbulence_handling = 'none' it works fine though, this only happen with turbulence_handling = 'mixing-length'

all the kernels/actions I use are as posted above, below is the size of my problem

Parallelism:
  Num Processors:          32
  Num Threads:             2

Mesh: 
  Parallel Type:           replicated
  Mesh Dimension:          3
  Spatial Dimension:       3
  Nodes:                   
    Total:                 100293
    Local:                 3535
    Min/Max/Avg:           2417/3649/3134
  Elems:                   
    Total:                 93400
    Local:                 2837
    Min/Max/Avg:           2837/2967/2918
  Num Subdomains:          3
  Num Partitions:          32
  Partitioner:             metis

Nonlinear System:
  Num DOFs:                201600
  Num Local DOFs:          0
  Variables:               { "superficial_vel_x" "superficial_vel_y" "superficial_vel_z" "pressure" "T_fluid" 
                             "scalar" } 
  Finite Element Types:    "MONOMIAL" 
  Approximation Orders:    "CONSTANT" 

Auxiliary System:
  Num DOFs:                93400
  Num Local DOFs:          2837
  Variables:               "porosity_output" 
  Finite Element Types:    "MONOMIAL" 
  Approximation Orders:    "CONSTANT" 

[GiudGiud]

try 300?

and did you re-compile after?

[Traiwit]

Hi @GiudGiud

300 works! yeh i didnt re-compile, my bad :)

just wondering is there any down side when i set this to --with-derivative-size=300 instead of the default value? will it make my simulation slower or taking more memory?

[GiudGiud]

taking more memory

anything that caches AD numbers (notably gradient caching in INSFV) will scale linearly from this number.

[Traiwit]

Hi @GiudGiud

I have another question re: friction_coeffs inside [Physics]/[NavierStokes]/[Flow]

I'm not able to assign value to this friction_coeffs at all (actually dont know exactly what this is)

first I tried a single value like

  [NavierStokes]
    [Flow]
      [flow]
        compressibility = 'incompressible'
        porous_medium_treatment = true
        porosity = porosity
        ...
        ...
        ...
        friction_blocks = '2'
        friction_types = 'DARCY'
        friction_coeffs = '1e-7'

        # transient = false
         block = '12 2'
      []
    []

  []

i got this error

No functor ever provided with name 'const_Fdrags_mat', which was requested by 'flow_momentum_friction_2_x,flow_momentum_friction_2_y,flow_momentum_friction_2_z,flow_pins_momentum_friction_correction_2_x,flow_pins_momentum_friction_correction_2_y,flow_pins_momentum_friction_correction_2_z'.

then i provided

        [./constant_vector]
          type = ParsedVectorFunction
          expression_x = '1e-7'
          expression_y = '1e-7'
          expression_z = '1e-7'
        [../]

       friction_coeffs = 'constant_vector'

then i got this error:

A call to SubProblem::getFunctor requested a functor named 'constant_vector' that returns the type: 'libMesh::VectorValue<MetaPhysicL::DualNumber<double, MetaPhysicL::SemiDynamicSparseNumberArray<double, unsigned long, MetaPhysicL::NWrapper<150ul> >, true> >'. However, that functor already exists and returns a different type: 'MetaPhysicL::DualNumber<double, MetaPhysicL::SemiDynamicSparseNumberArray<double, unsigned long, MetaPhysicL::NWrapper<150ul> >, true>'

I also tried

  [const_Fdrags_mat]
    type = FunctorErgunDragCoefficients
    porosity = porosity
    pebble_diameter = ${pebble_diameter}
    Darcy_multiplier = 1
  []

and still failed

I'm not sure what to do with this anymore, my ultimate objective is to control frictional force via Darcy coefficient explained in

https://mooseframework.inl.gov/source/fvkernels/PINSFVMomentumFriction.html

best would be by related to permeability rather than the pebble_diameter

please help when you have time

Thank you!
Traiwit

[GiudGiud]

Hello

Use the ParsedVectorFunctorMaterial to define this vector functor

Functions are Real valued functors at this time, even this vector function you used. We would need a new VectorFunction class to make it a vector-valued functor

Guillaume

[Traiwit]

thank you @GiudGiud

so i tried quite a few different value of Darcy_coefficient range from 1e-5 to 1e-15

the results seem to remain the same, like this value doesn't do anything

        friction_blocks = '2'
        friction_types = 'DARCY'
        friction_coeffs = 'Darcy_coefficient'

        # transient = false
         block = '12 2'

    [darcy]
      type = ADGenericVectorFunctorMaterial
      prop_names = 'Darcy_coefficient'
      prop_values = '1e-15 1e-15 1e-15'
    []

I just wonder if the friction_coeffs actually goes into PINSFVMomentumFriction

in the [Physics]/[NavierStokes]/[Flow] (WCNSFVFlowPhysics)

void
WCNSFVFlowPhysics::addMomentumFrictionKernels()
{

  std::cout << "DEBUG: Reached computeFrictionWCoefficient()" << std::endl;
  std::cout.flush();


  unsigned int num_friction_blocks = _friction_blocks.size();
  unsigned int num_used_blocks = num_friction_blocks ? num_friction_blocks : 1;

it prints the DEBUG but it doesn't print whatever inside PINSFVMomentumFriction

ADReal
PINSFVMomentumFriction::computeFrictionWCoefficient(const Moose::ElemArg & elem_arg,
                                                    const Moose::StateArg & state)
{

  std::cout << "DEBUG: PINSFVMomentumFriction::calculating Friction coefficient "  << std::endl;
  std::cout.flush();

  // Forward declaration of the coeffcients to be used and returned by the model
  ADReal coefficient = 0.0;
  ADReal speed = 0.0;
  ADReal rho = 0.0;
  ADReal mu = 0.0;
  /// Speed treatment
  if (_use_Forchheimer_friction_model && _speed)
    speed = (*_speed)(elem_arg, state);
  else if (_use_Forchheimer_friction_model && !_speed)
  {
    ADRealVectorValue superficial_velocity;
    if (_dim >= 1)
    {
      superficial_velocity(0) = (*_u_var)(elem_arg, state);
      if (_dim >= 2)
      {
        superficial_velocity(1) = (*_v_var)(elem_arg, state);
        if (_dim >= 3)
          superficial_velocity(2) = (*_w_var)(elem_arg, state);
      }
    }
    speed = NS::computeSpeed(superficial_velocity);
    if (_is_porous_medium)
    {
      speed *= (1 / _epsilon(elem_arg, state));
    }
  }

  /// Fluid properies
  if (_use_Darcy_friction_model && _standard_friction_formulation)
  {
    mu = (*_mu)(elem_arg, state);
  }

  if (_use_Forchheimer_friction_model && _standard_friction_formulation)
  {
    rho = _rho(elem_arg, state);
  }

  ////////////////////////////////////////////////////////////////////
  ///// Switching across formulation cases
  ////////////////////////////////////////////////////////////////////
  if (_standard_friction_formulation)
  {
    if (_use_Darcy_friction_model)
      coefficient += mu * (*_D)(elem_arg, state)(_index);
    if (_use_Forchheimer_friction_model)
      coefficient += rho / 2 * (*_F)(elem_arg, state)(_index)*speed;
  }
  else
  {
    if (_use_Darcy_friction_model)
      coefficient += (*_D)(elem_arg, state)(_index);
    if (_use_Forchheimer_friction_model)
      coefficient += (*_F)(elem_arg, state)(_index)*speed;
  }

  std::cout << "DEBUG: Friction coefficient = " << coefficient << std::endl;
  std::cout.flush();


  return coefficient;
}

[Traiwit]

okay i tried a higher value, and it makes the difference

so it's f_i (friction factor) not D (Darcy coefficient)

https://mooseframework.inl.gov/source/fvkernels/PINSFVMomentumFriction.html

[GiudGiud]

There are two formulations available. (Standard / false to use_standard). One should let you specify the Darcy coefficient while the other one lets you specify the linear coefficient when writing friction as F times velocity.

[Traiwit]

Hi @GiudGiud I have a couple of questions regarding Navier-Stokes in MOOSE (again):

1. Is this a deliberate choice to avoid over-constraining the system? In my case (mine ventilation/cave flow simulation), the actual control is usually flux-based (m³/s), not pressure. I understand the pressure outlet is more natural for NS, but is there any safe workaround to prescribe a known outlet flux or velocity without causing instabilities?

2. is there any good way to do quasi-steady state with NS, currently i just follow my Porousflow workflow, ie update property (porosity/darcy coeff/source) via PropertyReadFile time based, then set TimeSteppers as below, you see 10000 and 10010, so that i can reduce the DT when properties changes (so it's easier to converge), but just wondering if there is a better way to do this?

  [uobj_read_porosity_fromCsv]
    type = PropertyReadFileTimeBased   // (my own object)
    prop_file_name = 'exodus_meshes_with_nodesets_v5_smaller_stepped/XXXX_porosity1.csv
                      exodus_meshes_with_nodesets_v5_smaller_stepped/XXXX_porosity2.csv
                      exodus_meshes_with_nodesets_v5_smaller_stepped/XXXX_porosity3.csv'
    nprop = 1
    time_per_file = '10000 20000 30000'
    read_type = element
    execute_on = 'TIMESTEP_BEGIN'
  []

   [porosity_M0]
    type = FunctionAux
    variable = porosity_M0
    function = Porosity_csv_read_fx
    execute_on = 'initial timestep_begin'
  []

--------------------------------------

    [fixed_points]
      type = TimeSequenceStepper
      time_sequence = ' 10000 10010 20000 20010'
    []

    [TimeStepper_IterationAdaptiveDT]
      type = IterationAdaptiveDT
      dt = 0.5
    []

[GiudGiud]

Hello

    momentum_inlet_functors = 'inlet_vel_L inlet_vel_R'

should be

    momentum_inlet_functors = 'inlet_vel_L; inlet_vel_R'

[GiudGiud]

Is this a deliberate choice to avoid over-constraining the system? In my case (mine ventilation/cave flow simulation), the actual control is usually flux-based (m³/s), not pressure. I understand the pressure outlet is more natural for NS, but is there any safe workaround to prescribe a known outlet flux or velocity without causing instabilities?

There s a choice to avoid under-constraining the system. For incompressible flow, there is a null space for pressure, where pressure can be offset globally by any value and it s still a solution to the equation.
So if you do not want a pressure outlet, you will need to either pin (with a constraint, there are examples):

• the average pressure vaue
• the pressure value at a single point

The Physics/Action lets you do that (see pressure pin parameters).

[GiudGiud]

is there any good way to do quasi-steady state with NS, currently i just follow my Porousflow workflow, ie update property (porosity/darcy coeff/source) via PropertyReadFile time based, then set TimeSteppers as below, you see 10000 and 10010, so that i can reduce the DT when properties changes (so it's easier to converge), but just wondering if there is a better way to do this?

yes and no.
If it converges, you can do steady state solves by not including the time derivative. If using the Physics syntax, then set transient = false in the Physics.
Then time steps in a transient executioner become load steps.
If steady state solves dont converge, then you will need to use a Transient executioner, then make the PropertyReadFile read another file on every execution. The TimesEnableControl lets you set this up.
Some time steps will be used to converge a load step, some time steps will represent the change in load steps

[Traiwit]

Hi @GiudGiud, thank you for your reply

re: q1, so there is no way to work around setting outlet as flux/velocity at all?

im using flow Physics syntax and the only options i got for momentum_outlet_types are fixed-pressure, zero-gradient, fixed-pressure-zero-gradient

re: q2. yes, I tried transient = false but it wont converge, i guess I can work with this, just need to deal with adaptive timestepper better (lower DT then cave geometry/flow change)

[GiudGiud]

re: q1, so there is no way to work around setting outlet as flux/velocity at all?

there is. You can use the OutflowBC classes, such as
https://mooseframework.inl.gov/source/fvbcs/INSFVMassAdvectionOutflowBC.html
https://mooseframework.inl.gov/source/fvbcs/INSFVMomentumAdvectionOutflowBC.html
https://mooseframework.inl.gov/source/fvbcs/NSFVOutflowTemperatureBC.html

im using flow Physics syntax and the only options i got for momentum_outlet_types are fixed-pressure, zero-gradient, fixed-pressure-zero-gradient

zero-gradient is the one that allows the outflow BCs to be used.
I dont see the outflow in the energy physics, so you'll need to add NSFVOutflowTemperatureBC manually outside the Physics

re: q2. yes, I tried transient = false but it wont converge, i guess I can work with this, just need to deal with adaptive timestepper better (lower DT then cave geometry/flow change)

Steady states from fresh are hard with navier stokes. You could try initializing the first load step with a transient result solution, and see if the next load steps converge without help

otherwise a flexible time stepping and a bunch of Control objects to turn on/off various terms are the way to do this