Line data    Source code

       1             : //* This file is part of the MOOSE framework
       2             : //* https://www.mooseframework.org
       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 "DomainIntegralQFunction.h"
      11             : 
      12             : registerMooseObject("TensorMechanicsApp", DomainIntegralQFunction);
      13             : 
      14             : InputParameters
      15         703 : DomainIntegralQFunction::validParams()
      16             : {
      17         703 :   InputParameters params = AuxKernel::validParams();
      18         703 :   params.addClassDescription("Computes the q-function for a segment along the crack front, used in "
      19             :                              "the calculation of the J-integral");
      20        1406 :   params.addRequiredParam<Real>("j_integral_radius_inner", "Radius for J-Integral calculation");
      21        1406 :   params.addRequiredParam<Real>("j_integral_radius_outer", "Radius for J-Integral calculation");
      22        1406 :   params.addRequiredParam<UserObjectName>("crack_front_definition",
      23             :                                           "The CrackFrontDefinition user object name");
      24        1406 :   params.addParam<unsigned int>(
      25             :       "crack_front_point_index",
      26             :       "The index of the point on the crack front corresponding to this q function");
      27         703 :   params.set<bool>("use_displaced_mesh") = false;
      28         703 :   return params;
      29           0 : }
      30             : 
      31         624 : DomainIntegralQFunction::DomainIntegralQFunction(const InputParameters & parameters)
      32             :   : AuxKernel(parameters),
      33         624 :     _j_integral_radius_inner(getParam<Real>("j_integral_radius_inner")),
      34        1248 :     _j_integral_radius_outer(getParam<Real>("j_integral_radius_outer")),
      35         624 :     _crack_front_definition(&getUserObject<CrackFrontDefinition>("crack_front_definition")),
      36        1248 :     _has_crack_front_point_index(isParamValid("crack_front_point_index")),
      37         624 :     _crack_front_point_index(
      38        1041 :         _has_crack_front_point_index ? getParam<unsigned int>("crack_front_point_index") : 0),
      39         624 :     _treat_as_2d(false),
      40         624 :     _is_point_on_intersecting_boundary(false)
      41             : {
      42         624 : }
      43             : 
      44             : void
      45         624 : DomainIntegralQFunction::initialSetup()
      46             : {
      47         624 :   _treat_as_2d = _crack_front_definition->treatAs2D();
      48             :   bool using_mesh_cutter = _crack_front_definition->usingMeshCutter();
      49         624 :   if (_treat_as_2d && using_mesh_cutter == false)
      50             :   {
      51         207 :     if (_has_crack_front_point_index)
      52             :     {
      53           0 :       mooseWarning(
      54             :           "crack_front_point_index ignored because CrackFrontDefinition is set to treat as 2D");
      55             :     }
      56             :   }
      57             :   else
      58             :   {
      59         417 :     if (!_has_crack_front_point_index)
      60             :     {
      61           0 :       mooseError("crack_front_point_index must be specified in DomainIntegralQFunction");
      62             :     }
      63             :   }
      64         624 :   _is_point_on_intersecting_boundary =
      65         624 :       _crack_front_definition->isPointWithIndexOnIntersectingBoundary(_crack_front_point_index);
      66         624 : }
      67             : 
      68             : Real
      69      426612 : DomainIntegralQFunction::computeValue()
      70             : {
      71             :   Real dist_to_crack_front;
      72             :   Real dist_along_tangent;
      73      426612 :   projectToFrontAtPoint(dist_to_crack_front, dist_along_tangent);
      74             : 
      75             :   Real q = 1.0;
      76      426612 :   if (dist_to_crack_front > _j_integral_radius_inner &&
      77      419284 :       dist_to_crack_front < _j_integral_radius_outer)
      78        2112 :     q = (_j_integral_radius_outer - dist_to_crack_front) /
      79        2112 :         (_j_integral_radius_outer - _j_integral_radius_inner);
      80      424500 :   else if (dist_to_crack_front >= _j_integral_radius_outer)
      81             :     q = 0.0;
      82             : 
      83      426612 :   if (q > 0.0)
      84             :   {
      85        9432 :     Real tangent_multiplier = 1.0;
      86        9432 :     if (!_treat_as_2d)
      87             :     {
      88             :       const Real forward_segment_length =
      89        5140 :           _crack_front_definition->getCrackFrontForwardSegmentLength(_crack_front_point_index);
      90             :       const Real backward_segment_length =
      91        5140 :           _crack_front_definition->getCrackFrontBackwardSegmentLength(_crack_front_point_index);
      92             : 
      93        5140 :       if (dist_along_tangent >= 0.0)
      94             :       {
      95        3328 :         if (forward_segment_length > 0.0)
      96        3112 :           tangent_multiplier = 1.0 - dist_along_tangent / forward_segment_length;
      97             :       }
      98             :       else
      99             :       {
     100        1812 :         if (backward_segment_length > 0.0)
     101        1784 :           tangent_multiplier = 1.0 + dist_along_tangent / backward_segment_length;
     102             :       }
     103             :     }
     104             : 
     105        9432 :     tangent_multiplier = std::max(tangent_multiplier, 0.0);
     106        9432 :     tangent_multiplier = std::min(tangent_multiplier, 1.0);
     107             : 
     108             :     // Set to zero if a node is on a designated free surface and its crack front node is not.
     109        9432 :     if (_crack_front_definition->isNodeOnIntersectingBoundary(_current_node) &&
     110         948 :         !_is_point_on_intersecting_boundary)
     111         708 :       tangent_multiplier = 0.0;
     112             : 
     113        9432 :     q *= tangent_multiplier;
     114             :   }
     115             : 
     116      426612 :   return q;
     117             : }
     118             : 
     119             : void
     120      426612 : DomainIntegralQFunction::projectToFrontAtPoint(Real & dist_to_front, Real & dist_along_tangent)
     121             : {
     122             :   const Point * crack_front_point =
     123      426612 :       _crack_front_definition->getCrackFrontPoint(_crack_front_point_index);
     124             : 
     125      426612 :   Point p = *_current_node;
     126             :   const RealVectorValue & crack_front_tangent =
     127      426612 :       _crack_front_definition->getCrackFrontTangent(_crack_front_point_index);
     128             : 
     129             :   RealVectorValue crack_node_to_current_node = p - *crack_front_point;
     130      426612 :   dist_along_tangent = crack_node_to_current_node * crack_front_tangent;
     131             :   RealVectorValue projection_point = *crack_front_point + dist_along_tangent * crack_front_tangent;
     132             :   RealVectorValue axis_to_current_node = p - projection_point;
     133      426612 :   dist_to_front = axis_to_current_node.norm();
     134      426612 : }