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 "SalehaniIrani3DCTraction.h"
      11             : 
      12             : registerMooseObject("SolidMechanicsApp", SalehaniIrani3DCTraction);
      13             : 
      14             : InputParameters
      15         100 : SalehaniIrani3DCTraction::validParams()
      16             : {
      17         100 :   InputParameters params = CZMComputeLocalTractionTotalBase::validParams();
      18         100 :   params.addClassDescription("3D Coupled (3DC) cohesive law of Salehani and Irani with no damage");
      19         200 :   params.addRequiredParam<Real>(
      20             :       "normal_gap_at_maximum_normal_traction",
      21             :       "The value of normal gap at which maximum normal traction is achieved");
      22         200 :   params.addRequiredParam<Real>(
      23             :       "tangential_gap_at_maximum_shear_traction",
      24             :       "The value of tangential gap at which maximum shear traction is achieved");
      25         200 :   params.addRequiredParam<Real>("maximum_normal_traction",
      26             :                                 "The maximum normal traction the interface can sustain");
      27         200 :   params.addRequiredParam<Real>("maximum_shear_traction",
      28             :                                 "The maximum shear traction the interface can sustain");
      29         100 :   return params;
      30           0 : }
      31             : 
      32          48 : SalehaniIrani3DCTraction::SalehaniIrani3DCTraction(const InputParameters & parameters)
      33             :   : CZMComputeLocalTractionTotalBase(parameters),
      34          96 :     _delta_u0({getParam<Real>("normal_gap_at_maximum_normal_traction"),
      35          96 :                std::sqrt(2) * getParam<Real>("tangential_gap_at_maximum_shear_traction"),
      36          96 :                std::sqrt(2) * getParam<Real>("tangential_gap_at_maximum_shear_traction")}),
      37         192 :     _max_allowable_traction({getParam<Real>("maximum_normal_traction"),
      38             :                              getParam<Real>("maximum_shear_traction"),
      39          48 :                              getParam<Real>("maximum_shear_traction")})
      40             : {
      41          48 : }
      42             : 
      43             : void
      44       13344 : SalehaniIrani3DCTraction::computeInterfaceTractionAndDerivatives()
      45             : {
      46       13344 :   _interface_traction[_qp] = computeTraction();
      47       13344 :   _dinterface_traction_djump[_qp] = computeTractionDerivatives();
      48       13344 : }
      49             : 
      50             : RealVectorValue
      51       13344 : 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       53376 :   for (i = 0; i < 3; i++)
      65             :   {
      66       40032 :     aa = _interface_displacement_jump[_qp](i) / _delta_u0(i);
      67       40032 :     if (i > 0)
      68       26688 :       aa *= aa; // square for shear component
      69             : 
      70       40032 :     x += aa;
      71             :   }
      72             : 
      73       13344 :   exp_x = std::exp(-x);
      74             : 
      75       53376 :   for (i = 0; i < 3; i++)
      76             :   {
      77       40032 :     if (i == 0)
      78             :       aa = std::exp(1);
      79             :     else
      80             :       aa = std::sqrt(2 * std::exp(1));
      81             : 
      82       40032 :     a_i = _max_allowable_traction(i) * aa;
      83       40032 :     b_i = _interface_displacement_jump[_qp](i) / _delta_u0(i);
      84       40032 :     traction_local(i) = a_i * b_i * exp_x;
      85             :   }
      86             : 
      87       13344 :   return traction_local;
      88             : }
      89             : 
      90             : RankTwoTensor
      91       13344 : SalehaniIrani3DCTraction::computeTractionDerivatives()
      92             : {
      93       13344 :   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       53376 :   for (i = 0; i < 3; i++)
     116             :   {
     117       40032 :     aa = _interface_displacement_jump[_qp](i) / _delta_u0(i);
     118       40032 :     if (i > 0)
     119       26688 :       aa *= aa;
     120       40032 :     x += aa;
     121             :   }
     122       13344 :   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       53376 :   for (i = 0; i < 3; i++)
     132             :   {
     133       40032 :     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       40032 :     a_i *= _max_allowable_traction(i);
     139       40032 :     b_i = _interface_displacement_jump[_qp](i) / _delta_u0(i);
     140             : 
     141      160128 :     for (j = 0; j < 3; j++)
     142             :     {
     143             : 
     144             :       dbi_dui = 0;
     145      120096 :       if (i == j)
     146       40032 :         dbi_dui = 1 / _delta_u0(j);
     147             : 
     148      120096 :       if (j == 0) // alpha = 1
     149       40032 :         dx_duj = 1. / _delta_u0(j);
     150             :       else // alpha = 2
     151       80064 :         dx_duj = 2. * _interface_displacement_jump[_qp](j) / (_delta_u0(j) * _delta_u0(j));
     152             : 
     153      120096 :       traction_jump_derivatives_local(i, j) =
     154      120096 :           a_i * exp_x * (dbi_dui - b_i * dx_duj); // the minus sign is due to exp(-x)
     155             :     }
     156             :   }
     157             : 
     158       13344 :   return traction_jump_derivatives_local;
     159             : }