Implementing user-defined Permeability Material in PorousFlow

[thewhitewalker10]

Hello,

I’m following the layout of Kozeny-Carman permeability material in PorousFlow to develop my own permeability material. I’m kinda stuck on how the Jacobian/derivative of the K-C material was derived w.r.t to the Variables (i.e., _dpermeability_qp_dvar[_qp][v]). Can someone pls walk me through it? Btw, is there a way to add these derivatives automatically or you need to manually compute and implement these derivatives? My permeability model is not dependent on the porosity, but on the strain as follows:
k=k_m*I + b/a [(b^2/12 - k_m)] (I-M), where b=b_0 + /Deltab and /Deltab=a 〈/epsilon_n -/epsilon_0〉.
I already derived the derivative of k with respect to the strain (/epsilon) but I’m not sure if I need to compute and implement the derivative of k with respect to porosity. Thanks.

[thewhitewalker10]

Hi, @cpgr and @WilkAndy

I’m trying to implement the permeability model described on page 6 of Zill et al. (https://doi.org/10.1016/j.ijrmms.2021.104879). I’m following the layout of PorousFlowPermeabilityKozenyCarman. Is this the right approach to implement a user-defined permeability in PorousFlow? or do I have to follow the layout of PorousFlowPorosityLinear since my permeability does not depend on porosity? I got stuck with computing the derivatives while following PorousFlowPermeabilityKozenyCarman. If I have to go with the PorousFlowPermeabilityKozenyCarman, do I need to find the derivative wrt porosity or the strain? Thanks.

[cpgr]

If permeability doesn't depend on porosity, then the simpler option of following PorousFlowPorosityLinear should be ok (ie you get the strain and it's derivatives as material properties, and use them when calculating your permeability.

PorousFlow isn't set up to make full use of AD yet, so this would need to be by hand.

[thewhitewalker10]

Thanks, @cpgr.

Just a couple of questions: What material provides the derivatives of the porosity (i.e.,_dporosity_qp_dvar and _dporosity_qp_dgradvar ) on lines 66 and 67 of PorousFlowPermeabilityKozenyCarman? I’m guessing that I can use the following derivatives for the strain in my implementation: _dvol_strain_qp_dvar and _dvol_strain_qp_dgradvar. Correct?
Why were _dpermeability_qp_dvar[_qp][v] and _dpermeability_qp_dgradvar[_qp][i][v] calculated on lines 100 and 106 in the PorousFlowPermeabilityKozenyCarman permeability model? I’m trying to implement _dpermeability_qp_dgradvar[_qp][i][v] in my own model, but I can’t find _dvol_strain_qp_dgradvar in any of the source or header files in PorousFlow. How do I obtain _dvol_strain_qp_dgradvar? Thanks.

[WilkAndy]

This "answer" just collects a variety of my thoughts - hopefully one might be useful.

• Your formula for permeability (Zill et al) means that it depends on strain in the normal direction to the fracture. This is a little different than the other perm formulations, which assume that permeability depends on volumetric strain. When creating PorousFlow initially, we considered the idea of allowing Material properties to be functions of arbitrary derivatives of variables, but rejected the idea as being too costly, both in memory and compute (almost every model would repeatedly just compute a bunch of zeroes for the derivatives).
• This is possibly the most obscure and complicated part of PorousFlow, sorry, because the derivatives are a bit mind-bending
• You do not need to calculate the derivatives for PorousFlow to work. If you don't, they will be assumed to be zero and you will get an inexact Jacobian and poorer convergence, BUT maybe that doesn't really matter in your case.
• If you look at

  moose/modules/porous_flow/include/materials/PorousFlowVolumetricStrain.h

  Lines 57 to 62 in ecb6d9c

  	/**
  	* The derivative of the total volumetric strain with respect to the porous flow variables.
  	* Since the total volumetric strain depends on derivatives of the displacement variables,
  	* this should be multiplied by _grad_phi in kernels
  	*/
  	MaterialProperty<std::vector<RealGradient>> & _dvol_total_strain_qp_dvar;

  and

  moose/modules/porous_flow/src/materials/PorousFlowVolumetricStrain.C

  Lines 47 to 48 in ecb6d9c

  	_dvol_total_strain_qp_dvar(
  	declareProperty<std::vector<RealGradient>>("dPorousFlow_total_volumetric_strain_qp_dvar"))

  and the code for calculating this object

  moose/modules/porous_flow/src/materials/PorousFlowVolumetricStrain.C

  Lines 71 to 79 in ecb6d9c

  	_dvol_total_strain_qp_dvar[_qp].resize(_num_var, RealGradient());
  	for (unsigned i = 0; i < _ndisp; ++i)
  	if (_dictator.isPorousFlowVariable(_disp_var_num[i]))
  	{
  	// the i_th displacement is a PorousFlow variable
  	const unsigned int pvar = _dictator.porousFlowVariableNum(_disp_var_num[i]);
  	_dvol_strain_rate_qp_dvar[_qp][pvar](i) = 1.0 / _dt;
  	_dvol_total_strain_qp_dvar[_qp][pvar](i) = 1.0;
  	}

  you'll see that the Material Property called dPorousFlow_total_volumetric_strain_qp_dvar , is "d(volumetric strain) / d(displacement)". I write this in quotes because it's not strictly notationally correct (and that's why the object is a std::vector<RealGradient> not just a Real) To understand further, and understand why this dPorousFlow_total_volumetric_strain_qp_dvar is needed, imagine a Kernel that is exactly residual = volumetric_strain. Then, a Jacobian entry is d(residual) / d(displacement_j) = 1 * grad_j (shape function) . If you don't understand this, please ask, as this is the key point. The "1" in the previous formula is a real number (1) and is exactly a component of dPorousFlow_total_volumetric_strain_qp_dvar. More generally, imagine a Kernel that is residual = f(volumetric_strain). Then, a Jacobian entry is d(residual) / d(displacement_j) = sum_over_i f'(volumetric_strain) * dPorousFlow_total_volumetric_strain_qp_dvar(j, i) grad_i (shape function). In MOOSE language, the "j" is the std::vector part, and the "i" is the RealGradient, and grad_i (shape function) is grad_phi.
• Probably the above will explain why there is no _dvol_strain_qp_dgradvar. The Material called dPorousFlow_total_volumetric_strain_qp_dvar` is closest to what you want.
• If you look at

  moose/modules/porous_flow/include/materials/PorousFlowPermeabilityBase.h

  Lines 28 to 32 in ecb6d9c

  	/// d(quadpoint permeability)/d(PorousFlow variable)
  	MaterialProperty<std::vector<RealTensorValue>> & _dpermeability_qp_dvar;
  	/// d(quadpoint permeability)/d(grad(PorousFlow variable))
  	MaterialProperty<std::vector<std::vector<RealTensorValue>>> & _dpermeability_qp_dgradvar;

  you'll see the definition d(perm)/d(var) and d(perm)/d(gradvar). Consider the case that kernel = perm, so Jacobian = d(perm)/d(var) * phi + d(perm)/d(gradvar) * grad_phi. Hopefully this makes sense: the key is, we need derivatives of Materials with respect to variables, and the gradients of variables, because we need to know whether the Jacobian entries should be multiplied by phi or grad_phi.

[thewhitewalker10]

Thanks, @WilkAndy

From your feedback, I don’t need the derivatives after all. I’m going to go ahead and implement the permeability following PorousFlowPermeabilityKozenyCarman w/o them. Please, keep this ticket open for now, as I need to sit down on your feedback to grasp it fully. I’m going to be chipping in some questions from time to time. Thanks a lot guys @WilkAndy @cpgr !