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 "SalehaniIrani3DCTraction.h"
      11             : 
      12             : registerMooseObject("TensorMechanicsApp", SalehaniIrani3DCTraction);
      13             : 
      14             : InputParameters
      15          50 : SalehaniIrani3DCTraction::validParams()
      16             : {
      17          50 :   InputParameters params = CZMComputeLocalTractionTotalBase::validParams();
      18          50 :   params.addClassDescription("3D Coupled (3DC) cohesive law of Salehani and Irani with no damage");
      19         100 :   params.addRequiredParam<Real>(
      20             :       "normal_gap_at_maximum_normal_traction",
      21             :       "The value of normal gap at which maximum normal traction is achieved");
      22         100 :   params.addRequiredParam<Real>(
      23             :       "tangential_gap_at_maximum_shear_traction",
      24             :       "The value of tangential gap at which maximum shear traction is achieved");
      25         100 :   params.addRequiredParam<Real>("maximum_normal_traction",
      26             :                                 "The maximum normal traction the interface can sustain");
      27         100 :   params.addRequiredParam<Real>("maximum_shear_traction",
      28             :                                 "The maximum shear traction the interface can sustain");
      29          50 :   return params;
      30           0 : }
      31             : 
      32          24 : SalehaniIrani3DCTraction::SalehaniIrani3DCTraction(const InputParameters & parameters)
      33             :   : CZMComputeLocalTractionTotalBase(parameters),
      34          48 :     _delta_u0({getParam<Real>("normal_gap_at_maximum_normal_traction"),
      35          48 :                std::sqrt(2) * getParam<Real>("tangential_gap_at_maximum_shear_traction"),
      36          48 :                std::sqrt(2) * getParam<Real>("tangential_gap_at_maximum_shear_traction")}),
      37          96 :     _max_allowable_traction({getParam<Real>("maximum_normal_traction"),
      38             :                              getParam<Real>("maximum_shear_traction"),
      39          24 :                              getParam<Real>("maximum_shear_traction")})
      40             : {
      41          24 : }
      42             : 
      43             : void
      44        7480 : SalehaniIrani3DCTraction::computeInterfaceTractionAndDerivatives()
      45             : {
      46        7480 :   _interface_traction[_qp] = computeTraction();
      47        7480 :   _dinterface_traction_djump[_qp] = computeTractionDerivatives();
      48        7480 : }
      49             : 
      50             : RealVectorValue
      51        7480 : SalehaniIrani3DCTraction::computeTraction()
      52             : {
      53             : 
      54             :   // The convention for ordering the traction is N, T, S, where N is the normal direction, and T and
      55             :   // S are two arbitrary tangential directions.
      56             :   RealVectorValue traction_local;
      57             : 
      58             :   // temporary containers for auxiliary calculations
      59             :   Real aa, x, exp_x, a_i, b_i;
      60             :   // index variable to avoid multiple redefinitions
      61             :   unsigned int i;
      62             : 
      63             :   x = 0;
      64       29920 :   for (i = 0; i < 3; i++)
      65             :   {
      66       22440 :     aa = _interface_displacement_jump[_qp](i) / _delta_u0(i);
      67       22440 :     if (i > 0)
      68       14960 :       aa *= aa; // square for shear component
      69             : 
      70       22440 :     x += aa;
      71             :   }
      72             : 
      73        7480 :   exp_x = std::exp(-x);
      74             : 
      75       29920 :   for (i = 0; i < 3; i++)
      76             :   {
      77       22440 :     if (i == 0)
      78             :       aa = std::exp(1);
      79             :     else
      80             :       aa = std::sqrt(2 * std::exp(1));
      81             : 
      82       22440 :     a_i = _max_allowable_traction(i) * aa;
      83       22440 :     b_i = _interface_displacement_jump[_qp](i) / _delta_u0(i);
      84       22440 :     traction_local(i) = a_i * b_i * exp_x;
      85             :   }
      86             : 
      87        7480 :   return traction_local;
      88             : }
      89             : 
      90             : RankTwoTensor
      91        7480 : SalehaniIrani3DCTraction::computeTractionDerivatives()
      92             : {
      93        7480 :   RankTwoTensor traction_jump_derivatives_local;
      94             : 
      95             :   // this function computes the partial derivates of Tn[0][:], Tt[1][:], Ts[2][:]
      96             :   // wrt dun, dut, dus
      97             :   // T_i = a_i*b_i*exp(-x) with:
      98             :   // a_i = \sigma_i,max * (\alpha_i*e)^{1/\alpha_i} with \alpha_i = 1 for i==n
      99             :   // \alpha_i = 2 for i!=n
     100             :   // b_i = \delta_u,i / \delta_0,i
     101             :   // x = sum_i=1^3{(\delta_u,i / \delta_0,i)^\alpha_i}  with \alpha_i = 1 for i==n
     102             :   // \alpha_i = 2 for i!=n
     103             : 
     104             :   // dTi_duj = a_i * ( dBi_duj * exp(-x) + b_i * exp(-x) * dx_duj  )
     105             :   //         = a_i * ( exp(-x) * (dBi_duj + b_i * dx_duj ) )
     106             : 
     107             :   // temporary containers for auxiliary calculations
     108             :   Real aa, exp_x, x;
     109             :   // index variables to avoid multiple redefinitions
     110             :   unsigned int i, j;
     111             : 
     112             :   // compute x and the exponential term
     113             :   aa = 0;
     114             :   x = 0;
     115       29920 :   for (i = 0; i < 3; i++)
     116             :   {
     117       22440 :     aa = _interface_displacement_jump[_qp](i) / _delta_u0(i);
     118       22440 :     if (i > 0)
     119       14960 :       aa *= aa;
     120       22440 :     x += aa;
     121             :   }
     122        7480 :   exp_x = std::exp(-x);
     123             : 
     124             :   // compute partial derivatives in local coordiante wrt the displacement jump
     125             :   //            | dTn/dun dTn/dut dTn/dus |
     126             :   // dTi_duj  = | dTt/dun dTt/dut dTt/dus | = _traction_derivative[i][j]
     127             :   //            | dTs/dun dTs/dut dTs/dus |
     128             :   Real a_i, b_i;
     129             :   Real dbi_dui, dx_duj;
     130             : 
     131       29920 :   for (i = 0; i < 3; i++)
     132             :   {
     133       22440 :     if (i == 0) // alpha = 1
     134             :       a_i = std::exp(1);
     135             :     else // alpha = 2
     136             :       a_i = std::sqrt(2 * std::exp(1));
     137             : 
     138       22440 :     a_i *= _max_allowable_traction(i);
     139       22440 :     b_i = _interface_displacement_jump[_qp](i) / _delta_u0(i);
     140             : 
     141       89760 :     for (j = 0; j < 3; j++)
     142             :     {
     143             : 
     144             :       dbi_dui = 0;
     145       67320 :       if (i == j)
     146       22440 :         dbi_dui = 1 / _delta_u0(j);
     147             : 
     148       67320 :       if (j == 0) // alpha = 1
     149       22440 :         dx_duj = 1. / _delta_u0(j);
     150             :       else // alpha = 2
     151       44880 :         dx_duj = 2. * _interface_displacement_jump[_qp](j) / (_delta_u0(j) * _delta_u0(j));
     152             : 
     153       67320 :       traction_jump_derivatives_local(i, j) =
     154       67320 :           a_i * exp_x * (dbi_dui - b_i * dx_duj); // the minus sign is due to exp(-x)
     155             :     }
     156             :   }
     157             : 
     158        7480 :   return traction_jump_derivatives_local;
     159             : }