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