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 "INSProjection.h"
      11             : #include "MooseMesh.h"
      12             : #include "NS.h"
      13             : 
      14             : registerMooseObject("NavierStokesApp", INSProjection);
      15             : 
      16             : InputParameters
      17           0 : INSProjection::validParams()
      18             : {
      19           0 :   InputParameters params = Kernel::validParams();
      20             : 
      21           0 :   params.addClassDescription("This class computes the 'projection' part of the 'split' method for "
      22             :                              "solving incompressible Navier-Stokes.");
      23             :   // Coupled variables
      24           0 :   params.addRequiredCoupledVar("a1", "x-acceleration");
      25           0 :   params.addCoupledVar("a2", "y-acceleration"); // only required in 2D and 3D
      26           0 :   params.addCoupledVar("a3", "z-acceleration"); // only required in 3D
      27           0 :   params.addRequiredCoupledVar(NS::pressure, "pressure");
      28             : 
      29             :   // Required parameters
      30           0 :   params.addRequiredParam<unsigned>(
      31             :       "component",
      32             :       "0,1,2 depending on if we are solving the x,y,z component of the momentum equation");
      33             : 
      34             :   // Optional parameters
      35           0 :   params.addParam<MaterialPropertyName>("rho_name", "rho", "density name");
      36             : 
      37           0 :   return params;
      38           0 : }
      39             : 
      40           0 : INSProjection::INSProjection(const InputParameters & parameters)
      41             :   : Kernel(parameters),
      42             : 
      43             :     // Coupled variables
      44           0 :     _a1(coupledValue("a1")),
      45           0 :     _a2(_mesh.dimension() >= 2 ? coupledValue("a2") : _zero),
      46           0 :     _a3(_mesh.dimension() == 3 ? coupledValue("a3") : _zero),
      47             : 
      48             :     // Gradients
      49           0 :     _grad_p(coupledGradient(NS::pressure)),
      50             : 
      51             :     // Variable numberings
      52           0 :     _a1_var_number(coupled("a1")),
      53           0 :     _a2_var_number(_mesh.dimension() >= 2 ? coupled("a2") : libMesh::invalid_uint),
      54           0 :     _a3_var_number(_mesh.dimension() == 3 ? coupled("a3") : libMesh::invalid_uint),
      55           0 :     _p_var_number(coupled(NS::pressure)),
      56             : 
      57             :     // Required parameters
      58           0 :     _component(getParam<unsigned>("component")),
      59             : 
      60             :     // Material properties
      61           0 :     _rho(getMaterialProperty<Real>("rho_name"))
      62             : {
      63           0 : }
      64             : 
      65             : Real
      66           0 : INSProjection::computeQpResidual()
      67             : {
      68             :   // Vector object for a
      69           0 :   RealVectorValue a(_a1[_qp], _a2[_qp], _a3[_qp]);
      70             : 
      71             :   // Vector object for test function (only the component'th entry is non-zero)
      72             :   RealVectorValue test;
      73           0 :   test(_component) = _test[_i][_qp];
      74             : 
      75             :   // "Symmetric" part, -a.test
      76           0 :   Real symmetric_part = -a(_component) * _test[_i][_qp];
      77             : 
      78             :   // The pressure part, (1/_rho[_qp]) * (grad(p).v)
      79           0 :   Real pressure_part = (1. / _rho[_qp]) * (_grad_p[_qp] * test);
      80             : 
      81             :   // Return the result
      82           0 :   return symmetric_part + pressure_part;
      83             : }
      84             : 
      85             : Real
      86           0 : INSProjection::computeQpJacobian()
      87             : {
      88             :   // There will be a diagonal component from the time derivative term...
      89           0 :   return 0.;
      90             : }
      91             : 
      92             : Real
      93           0 : INSProjection::computeQpOffDiagJacobian(unsigned jvar)
      94             : {
      95           0 :   if (((jvar == _a1_var_number) && (_component == 0)) ||
      96           0 :       ((jvar == _a2_var_number) && (_component == 1)) ||
      97           0 :       ((jvar == _a3_var_number) && (_component == 2)))
      98             :   {
      99             :     // The symmetric term's Jacobian is only non-zero when the
     100             :     // component of 'a' being differentiated is the same as _component.
     101           0 :     return -_phi[_j][_qp] * _test[_i][_qp];
     102             :   }
     103             : 
     104           0 :   else if (jvar == _p_var_number)
     105             :   {
     106           0 :     return (1. / _rho[_qp]) * (_grad_phi[_j][_qp](_component) * _test[_i][_qp]);
     107             :   }
     108             : 
     109             :   else
     110             :     return 0;
     111             : }