Line data    Source code

       1             : //* This file is part of the MOOSE framework
       2             : //* https://mooseframework.inl.gov
       3             : //*
       4             : //* All rights reserved, see COPYRIGHT for full restrictions
       5             : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
       6             : //*
       7             : //* Licensed under LGPL 2.1, please see LICENSE for details
       8             : //* https://www.gnu.org/licenses/lgpl-2.1.html
       9             : 
      10             : #include "PorousFlowDarcyVelocityComponentLowerDimensional.h"
      11             : #include "MooseMesh.h"
      12             : #include "Assembly.h"
      13             : 
      14             : registerMooseObject("PorousFlowApp", PorousFlowDarcyVelocityComponentLowerDimensional);
      15             : 
      16             : InputParameters
      17         456 : PorousFlowDarcyVelocityComponentLowerDimensional::validParams()
      18             : {
      19         456 :   InputParameters params = PorousFlowDarcyVelocityComponent::validParams();
      20         912 :   params.addCoupledVar("aperture", 1.0, "Aperture of the fracture");
      21         456 :   params.addClassDescription(
      22             :       "Darcy velocity on a lower-dimensional element embedded in a higher-dimensional mesh.  Units "
      23             :       "m^3.s^-1.m^-2, or m.s^-1.  Darcy velocity =  -(k_ij * krel /(mu * a) (nabla_j P - w_j)), "
      24             :       "where k_ij is the permeability tensor, krel is the relative permeability, mu is the fluid "
      25             :       "viscosity, P is the fluid pressure, a is the fracture aperture and w_j is the fluid weight. "
      26             :       " The difference between this AuxKernel and PorousFlowDarcyVelocity is that this one "
      27             :       "projects gravity along the element's tangent direction.  NOTE!  For a meaningful answer, "
      28             :       "your permeability tensor must NOT contain terms that rotate tangential vectors to "
      29             :       "non-tangential vectors.");
      30         456 :   return params;
      31           0 : }
      32             : 
      33         245 : PorousFlowDarcyVelocityComponentLowerDimensional::PorousFlowDarcyVelocityComponentLowerDimensional(
      34         245 :     const InputParameters & parameters)
      35         245 :   : PorousFlowDarcyVelocityComponent(parameters), _aperture(coupledValue("aperture"))
      36             : {
      37         245 :   if (isNodal())
      38           0 :     paramError("variable", "This AuxKernel only supports Elemental fields");
      39         245 : }
      40             : 
      41             : Real
      42       62043 : PorousFlowDarcyVelocityComponentLowerDimensional::computeValue()
      43             : {
      44       62043 :   const unsigned elem_dim = _current_elem->dim();
      45       62043 :   if (elem_dim == _mesh.dimension())
      46           3 :     mooseError("The variable ",
      47           3 :                _var.name(),
      48             :                " must must be defined on lower-dimensional elements "
      49             :                "only since it employs "
      50             :                "PorousFlowDarcyVelocityComponentLowerDimensional\n");
      51             : 
      52             :   // MOOSE will automatically make grad(P) lie along the element's tangent direction
      53             :   // but we need to project gravity in that direction too
      54             : 
      55             :   // tang_xi is the element's tangent vector in xi direction
      56       62040 :   const std::vector<RealGradient> & tang_xi = _assembly.getFE(FEType(), elem_dim)->get_dxyzdxi();
      57             :   RealVectorValue tangential_gravity =
      58       62040 :       (_gravity * tang_xi[_qp] / tang_xi[_qp].norm_sq()) * tang_xi[_qp];
      59       62040 :   if (elem_dim == 2)
      60             :   {
      61             :     // tang_eta is the element's tangent vector in eta direction
      62             :     const std::vector<RealGradient> & tang_eta =
      63        1920 :         _assembly.getFE(FEType(), elem_dim)->get_dxyzdeta();
      64             :     const RealGradient normal_to_xi =
      65        1920 :         tang_eta[_qp] - (tang_eta[_qp] * tang_xi[_qp] / tang_xi[_qp].norm_sq()) * tang_xi[_qp];
      66        1920 :     tangential_gravity += (_gravity * normal_to_xi / normal_to_xi.norm_sq()) * normal_to_xi;
      67             :   }
      68             : 
      69       62040 :   return -(_permeability[_qp] *
      70       62040 :            (_grad_p[_qp][_ph] - _fluid_density_qp[_qp][_ph] * tangential_gravity) *
      71       62040 :            _relative_permeability[_qp][_ph] / _fluid_viscosity[_qp][_ph])(_component) /
      72       62040 :          _aperture[_qp];
      73             : }