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 "DomainIntegralQFunction.h"
      11             : 
      12             : registerMooseObject("SolidMechanicsApp", DomainIntegralQFunction);
      13             : 
      14             : InputParameters
      15        1406 : DomainIntegralQFunction::validParams()
      16             : {
      17        1406 :   InputParameters params = AuxKernel::validParams();
      18        1406 :   params.addClassDescription("Computes the q-function for a segment along the crack front, used in "
      19             :                              "the calculation of the J-integral");
      20        2812 :   params.addRequiredParam<Real>("j_integral_radius_inner", "Radius for J-Integral calculation");
      21        2812 :   params.addRequiredParam<Real>("j_integral_radius_outer", "Radius for J-Integral calculation");
      22        2812 :   params.addRequiredParam<UserObjectName>("crack_front_definition",
      23             :                                           "The CrackFrontDefinition user object name");
      24        2812 :   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        1406 :   params.set<bool>("use_displaced_mesh") = false;
      28        1406 :   return params;
      29           0 : }
      30             : 
      31        1248 : DomainIntegralQFunction::DomainIntegralQFunction(const InputParameters & parameters)
      32             :   : AuxKernel(parameters),
      33        1248 :     _j_integral_radius_inner(getParam<Real>("j_integral_radius_inner")),
      34        2496 :     _j_integral_radius_outer(getParam<Real>("j_integral_radius_outer")),
      35        1248 :     _crack_front_definition(&getUserObject<CrackFrontDefinition>("crack_front_definition")),
      36        2496 :     _has_crack_front_point_index(isParamValid("crack_front_point_index")),
      37        1248 :     _crack_front_point_index(
      38        2082 :         _has_crack_front_point_index ? getParam<unsigned int>("crack_front_point_index") : 0),
      39        1248 :     _treat_as_2d(false),
      40        1248 :     _is_point_on_intersecting_boundary(false)
      41             : {
      42        1248 : }
      43             : 
      44             : void
      45        1248 : DomainIntegralQFunction::initialSetup()
      46             : {
      47        1248 :   _treat_as_2d = _crack_front_definition->treatAs2D();
      48             :   bool using_mesh_cutter = _crack_front_definition->usingMeshCutter();
      49        1248 :   if (_treat_as_2d && using_mesh_cutter == false)
      50             :   {
      51         414 :     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         834 :     if (!_has_crack_front_point_index)
      60             :     {
      61           0 :       mooseError("crack_front_point_index must be specified in DomainIntegralQFunction");
      62             :     }
      63             :   }
      64        1248 :   _is_point_on_intersecting_boundary =
      65        1248 :       _crack_front_definition->isPointWithIndexOnIntersectingBoundary(_crack_front_point_index);
      66        1248 : }
      67             : 
      68             : Real
      69      853224 : DomainIntegralQFunction::computeValue()
      70             : {
      71             :   Real dist_to_crack_front;
      72             :   Real dist_along_tangent;
      73      853224 :   projectToFrontAtPoint(dist_to_crack_front, dist_along_tangent);
      74             : 
      75             :   Real q = 1.0;
      76      853224 :   if (dist_to_crack_front > _j_integral_radius_inner &&
      77      838568 :       dist_to_crack_front < _j_integral_radius_outer)
      78        4224 :     q = (_j_integral_radius_outer - dist_to_crack_front) /
      79        4224 :         (_j_integral_radius_outer - _j_integral_radius_inner);
      80      849000 :   else if (dist_to_crack_front >= _j_integral_radius_outer)
      81             :     q = 0.0;
      82             : 
      83      853224 :   if (q > 0.0)
      84             :   {
      85       18864 :     Real tangent_multiplier = 1.0;
      86       18864 :     if (!_treat_as_2d)
      87             :     {
      88             :       const Real forward_segment_length =
      89       10280 :           _crack_front_definition->getCrackFrontForwardSegmentLength(_crack_front_point_index);
      90             :       const Real backward_segment_length =
      91       10280 :           _crack_front_definition->getCrackFrontBackwardSegmentLength(_crack_front_point_index);
      92             : 
      93       10280 :       if (dist_along_tangent >= 0.0)
      94             :       {
      95        6656 :         if (forward_segment_length > 0.0)
      96        6224 :           tangent_multiplier = 1.0 - dist_along_tangent / forward_segment_length;
      97             :       }
      98             :       else
      99             :       {
     100        3624 :         if (backward_segment_length > 0.0)
     101        3568 :           tangent_multiplier = 1.0 + dist_along_tangent / backward_segment_length;
     102             :       }
     103             :     }
     104             : 
     105       18864 :     tangent_multiplier = std::max(tangent_multiplier, 0.0);
     106       18864 :     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       18864 :     if (_crack_front_definition->isNodeOnIntersectingBoundary(_current_node) &&
     110        1896 :         !_is_point_on_intersecting_boundary)
     111        1416 :       tangent_multiplier = 0.0;
     112             : 
     113       18864 :     q *= tangent_multiplier;
     114             :   }
     115             : 
     116      853224 :   return q;
     117             : }
     118             : 
     119             : void
     120      853224 : DomainIntegralQFunction::projectToFrontAtPoint(Real & dist_to_front, Real & dist_along_tangent)
     121             : {
     122             :   const Point * crack_front_point =
     123      853224 :       _crack_front_definition->getCrackFrontPoint(_crack_front_point_index);
     124             : 
     125      853224 :   Point p = *_current_node;
     126             :   const RealVectorValue & crack_front_tangent =
     127      853224 :       _crack_front_definition->getCrackFrontTangent(_crack_front_point_index);
     128             : 
     129             :   RealVectorValue crack_node_to_current_node = p - *crack_front_point;
     130      853224 :   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      853224 :   dist_to_front = axis_to_current_node.norm();
     134      853224 : }