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 "InclusionProperties.h"
      11             : #include "libmesh/utility.h"
      12             : 
      13             : registerMooseObject("TensorMechanicsApp", InclusionProperties);
      14             : 
      15             : InputParameters
      16           0 : InclusionProperties::validParams()
      17             : {
      18           0 :   InputParameters params = Material::validParams();
      19           0 :   params.addRequiredParam<Real>("a", "Ellipse semiaxis");
      20           0 :   params.addRequiredParam<Real>("b", "Ellipse semiaxis");
      21           0 :   params.addRequiredParam<Real>("lambda", "Lame's first parameter");
      22           0 :   params.addRequiredParam<Real>("mu", "Shear modulus (aka Lame's second parameter)");
      23           0 :   params.addRequiredParam<std::vector<Real>>("misfit_strains",
      24             :                                              "Vector of misfit strains in order eps_11, eps_22");
      25           0 :   params.addRequiredParam<MaterialPropertyName>(
      26             :       "stress_name", "Name of the material property where analytical stresses will be stored");
      27           0 :   params.addRequiredParam<MaterialPropertyName>(
      28             :       "strain_name", "Name of the material property where analytical total strains will be stored");
      29           0 :   params.addRequiredParam<MaterialPropertyName>(
      30             :       "energy_name",
      31             :       "Name of the material property where analytical elastic energies will be stored");
      32             : 
      33           0 :   return params;
      34           0 : }
      35             : 
      36           0 : InclusionProperties::InclusionProperties(const InputParameters & parameters)
      37             :   : Material(parameters),
      38           0 :     _a(getParam<Real>("a")),
      39           0 :     _b(getParam<Real>("b")),
      40           0 :     _lambda(getParam<Real>("lambda")),
      41           0 :     _mu(getParam<Real>("mu")),
      42           0 :     _misfit(getParam<std::vector<Real>>("misfit_strains")),
      43           0 :     _stress(declareProperty<RankTwoTensor>(getParam<MaterialPropertyName>("stress_name"))),
      44           0 :     _strain(declareProperty<RankTwoTensor>(getParam<MaterialPropertyName>("strain_name"))),
      45           0 :     _elastic_energy(declareProperty<Real>(getParam<MaterialPropertyName>("energy_name")))
      46             : {
      47           0 :   if (_misfit.size() != 2)
      48           0 :     mooseError("Supply 2 misfit_strains in order eps_11, eps_22 in InclusionProperties.");
      49             : 
      50           0 :   _nu = _lambda / 2.0 / (_lambda + _mu);
      51           0 :   _kappa = 3 - 4 * _nu;
      52           0 :   precomputeInteriorProperties();
      53           0 : }
      54             : 
      55             : void
      56           0 : InclusionProperties::computeQpProperties()
      57             : {
      58           0 :   Real x = _q_point[_qp](0);
      59           0 :   Real y = _q_point[_qp](1);
      60           0 :   if (x * x / _a / _a + y * y / _b / _b < 1)
      61             :   {
      62             :     // Inside the inclusion
      63           0 :     _stress[_qp] = _stress_int;
      64           0 :     _strain[_qp] = _total_strain_int;
      65           0 :     _elastic_energy[_qp] = _elastic_energy_int;
      66             :   }
      67             :   else
      68             :   {
      69             :     // Outside the inclusion
      70           0 :     Real l = 0.5 * (x * x + y * y - _a * _a - _b * _b // Parameter l called lambda in the paper
      71           0 :                     + std::sqrt(Utility::pow<2>((x * x + y * y - _a * _a + _b * _b)) +
      72           0 :                                 4 * (_a * _a - _b * _b) * y * y));
      73           0 :     Real rho_a = _a / sqrt(_a * _a + l);
      74           0 :     Real rho_b = _b / sqrt(_b * _b + l);
      75           0 :     Real m_x = x / (_a * _a + l);
      76           0 :     Real m_y = y / (_b * _b + l);
      77           0 :     Real n_x = m_x / std::sqrt(m_x * m_x + m_y * m_y);
      78           0 :     Real n_y = m_y / std::sqrt(m_x * m_x + m_y * m_y);
      79           0 :     Real T_6 = rho_a * rho_a + rho_b * rho_b - 4 * rho_a * rho_a * n_x * n_x -
      80           0 :                4 * rho_b * rho_b * n_y * n_y - 4;
      81             : 
      82           0 :     Real H11 = rho_a * _b *
      83           0 :                    (_a * rho_b + _b * rho_a + 2 * _a * rho_a * rho_a * rho_b +
      84           0 :                     _b * Utility::pow<3>(rho_a)) /
      85             :                    Utility::pow<2>((_a * rho_b + _b * rho_a)) +
      86           0 :                n_x * n_x * (2 - 6 * rho_a * rho_a + (8 * rho_a * rho_a + T_6) * n_x * n_x);
      87             : 
      88           0 :     Real H22 = rho_b * _a *
      89           0 :                    (_a * rho_b + _b * rho_a + 2 * _b * rho_a * rho_b * rho_b +
      90           0 :                     _a * Utility::pow<3>(rho_b)) /
      91             :                    Utility::pow<2>((_a * rho_b + _b * rho_a)) +
      92           0 :                n_y * n_y * (2 - 6 * rho_b * rho_b + (8 * rho_b * rho_b + T_6) * n_y * n_y);
      93             : 
      94           0 :     Real H12 = (_a * _a * rho_a * rho_a * rho_b * rho_b + _b * _b * rho_a * rho_a +
      95           0 :                 _a * _b * rho_a * rho_b) /
      96           0 :                    Utility::pow<2>((_a * rho_b + _b * rho_a)) -
      97           0 :                rho_b * rho_b * n_x * n_x - rho_a * rho_a * n_y * n_y +
      98           0 :                (4 * rho_a * rho_a + 4 * rho_b * rho_b + T_6) * n_x * n_x * n_y * n_y;
      99             : 
     100           0 :     Real H31 = 2 * (_b * rho_a / (_a * rho_b + _b * rho_a) - n_x * n_x);
     101           0 :     Real H32 = 2 * (_a * rho_b / (_a * rho_b + _b * rho_a) - n_y * n_y);
     102             : 
     103           0 :     Real H41 = n_x * n_y *
     104           0 :                (1 - 3 * rho_a * rho_a + (6 * rho_a * rho_a + 2 * rho_b * rho_b + T_6) * n_x * n_x);
     105           0 :     Real H42 = n_x * n_y *
     106           0 :                (1 - 3 * rho_b * rho_b + (2 * rho_a * rho_a + 6 * rho_b * rho_b + T_6) * n_y * n_y);
     107             : 
     108           0 :     _stress[_qp](0, 0) =
     109           0 :         4 * _mu * rho_a * rho_b / (_kappa + 1) * (H11 * _misfit[0] + H12 * _misfit[1]);
     110           0 :     _stress[_qp](1, 1) =
     111           0 :         4 * _mu * rho_a * rho_b / (_kappa + 1) * (H12 * _misfit[0] + H22 * _misfit[1]);
     112           0 :     _stress[_qp](2, 2) =
     113           0 :         4 * _mu * rho_a * rho_b / (_kappa + 1) * _nu * (H31 * _misfit[0] + H32 * _misfit[1]);
     114           0 :     _stress[_qp](0, 1) = _stress[_qp](1, 0) =
     115           0 :         4 * _mu * rho_a * rho_b / (_kappa + 1) * (H41 * _misfit[0] + H42 * _misfit[1]);
     116             : 
     117           0 :     Real J1 = rho_a * rho_a * rho_b * _b / (_a * rho_b + _b * rho_a);
     118           0 :     Real J11 = Utility::pow<4>(rho_a) * rho_b * _b / (3 * _a * _a) * (2 * _a * rho_b + _b * rho_a) /
     119           0 :                Utility::pow<2>((_a * rho_b + _b * rho_a));
     120           0 :     Real J12 = Utility::pow<3>(rho_a) * Utility::pow<3>(rho_b) /
     121           0 :                Utility::pow<2>((_a * rho_b + _b * rho_a));
     122           0 :     Real J2 = rho_b * rho_b * rho_a * _a / (_a * rho_b + _b * rho_a);
     123           0 :     Real J22 = Utility::pow<4>(rho_b) * rho_a * _a / (3 * _b * _b) * (2 * _b * rho_a + _a * rho_b) /
     124           0 :                Utility::pow<2>((_a * rho_b + _b * rho_a));
     125             : 
     126           0 :     Real G1111 = ((1 - 2 * _nu) * J1 + 3 * _a * _a * J11) / (2 * (1 - _nu)) +
     127           0 :                  rho_a * rho_b * n_x * n_x / (2 * (1 - _nu)) *
     128           0 :                      (2 + 2 * _nu - 6 * rho_a * rho_a + (8 * rho_a * rho_a + T_6) * n_x * n_x);
     129           0 :     Real G1122 = ((2 * _nu - 1) * J1 + _b * _b * J12) / (2 * (1 - _nu)) +
     130           0 :                  rho_a * rho_b / (2 * (1 - _nu)) *
     131           0 :                      ((1 - rho_a * rho_a) * n_y * n_y + (1 - 2 * _nu - rho_b * rho_b) * n_x * n_x +
     132             :                       (4 * rho_a * rho_a + 4 * rho_b * rho_b + T_6) * n_x * n_x * n_y * n_y);
     133           0 :     Real G2211 = ((2 * _nu - 1) * J2 + _a * _a * J12) / (2 * (1 - _nu)) +
     134           0 :                  rho_a * rho_b / (2 * (1 - _nu)) *
     135           0 :                      ((1 - rho_b * rho_b) * n_x * n_x + (1 - 2 * _nu - rho_a * rho_a) * n_y * n_y +
     136             :                       (4 * rho_a * rho_a + 4 * rho_b * rho_b + T_6) * n_x * n_x * n_y * n_y);
     137           0 :     Real G2222 = ((1 - 2 * _nu) * J2 + 3 * _b * _b * J22) / (2 * (1 - _nu)) +
     138           0 :                  rho_a * rho_b * n_y * n_y / (2 * (1 - _nu)) *
     139           0 :                      (2 + 2 * _nu - 6 * rho_b * rho_b + (8 * rho_a * rho_a + T_6) * n_y * n_y);
     140           0 :     Real G1211 =
     141           0 :         rho_a * rho_b * n_x * n_y / (2 * (1 - _nu)) *
     142             :         (1 - 3 * rho_a * rho_a + (6 * rho_a * rho_a + 2 * rho_b * rho_b + T_6) * n_x * n_x);
     143           0 :     Real G1222 =
     144             :         rho_a * rho_b * n_x * n_y / (2 * (1 - _nu)) *
     145             :         (1 - 3 * rho_b * rho_b + (2 * rho_a * rho_a + 6 * rho_b * rho_b + T_6) * n_y * n_y);
     146             : 
     147             :     // Outside the inclusion, total strain = elastic strain so we only need to
     148             :     // calculate one strain tensor
     149           0 :     _strain[_qp](0, 0) = G1111 * _misfit[0] + G1122 * _misfit[1];
     150           0 :     _strain[_qp](1, 1) = G2211 * _misfit[0] + G2222 * _misfit[1];
     151           0 :     _strain[_qp](0, 1) = _strain[_qp](1, 0) = G1211 * _misfit[0] + G1222 * _misfit[1];
     152             : 
     153           0 :     _elastic_energy[_qp] = 0.5 * _stress[_qp].doubleContraction(_strain[_qp]);
     154             :   }
     155           0 : }
     156             : 
     157             : void
     158           0 : InclusionProperties::precomputeInteriorProperties()
     159             : {
     160           0 :   Real t = _b / _a;
     161           0 :   Real T11 = 1 + 2 / t;
     162             :   Real T12 = 1;
     163             :   Real T21 = 1;
     164           0 :   Real T22 = 1 + 2 * t;
     165           0 :   Real T31 = (3 - _kappa) / 2 * (1 + 1 / t);
     166           0 :   Real T32 = (3 - _kappa) / 2 * (1 + t);
     167             : 
     168           0 :   _stress_int(0, 0) =
     169           0 :       -4 * _mu * t / (_kappa + 1) / (1 + t) / (1 + t) * (T11 * _misfit[0] + T12 * _misfit[1]);
     170           0 :   _stress_int(1, 1) =
     171           0 :       -4 * _mu * t / (_kappa + 1) / (1 + t) / (1 + t) * (T21 * _misfit[0] + T22 * _misfit[1]);
     172           0 :   _stress_int(2, 2) =
     173           0 :       -4 * _mu * t / (_kappa + 1) / (1 + t) / (1 + t) * (T31 * _misfit[0] + T32 * _misfit[1]);
     174             : 
     175           0 :   Real S11 = t * (t + 3 + _kappa * (1 + t));
     176           0 :   Real S12 = t * (1 + 3 * t - _kappa * (1 + t));
     177           0 :   Real S21 = t + 3 - _kappa * (1 + t);
     178           0 :   Real S22 = 1 + 3 * t + _kappa * (1 + t);
     179             : 
     180           0 :   _total_strain_int(0, 0) =
     181           0 :       1 / (_kappa + 1) / (1 + t) / (1 + t) * (S11 * _misfit[0] + S12 * _misfit[1]);
     182           0 :   _total_strain_int(1, 1) =
     183           0 :       1 / (_kappa + 1) / (1 + t) / (1 + t) * (S21 * _misfit[0] + S22 * _misfit[1]);
     184             :   // Inside the inclusion, elastic strain = total strain - eigenstrain
     185           0 :   _elastic_strain_int(0, 0) = _total_strain_int(0, 0) - _misfit[0];
     186           0 :   _elastic_strain_int(1, 1) = _total_strain_int(1, 1) - _misfit[1];
     187             : 
     188           0 :   _elastic_energy_int = 0.5 * _stress_int.doubleContraction(_elastic_strain_int);
     189           0 : }