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 "InterfaceOrientationMultiphaseMaterial.h"
      11             : #include "MooseMesh.h"
      12             : #include "MathUtils.h"
      13             : 
      14             : registerMooseObject("PhaseFieldApp", InterfaceOrientationMultiphaseMaterial);
      15             : 
      16             : InputParameters
      17         188 : InterfaceOrientationMultiphaseMaterial::validParams()
      18             : {
      19         188 :   InputParameters params = Material::validParams();
      20         188 :   params.addClassDescription(
      21             :       "This Material accounts for the the orientation dependence "
      22             :       "of interfacial energy for multi-phase multi-order parameter phase-field model.");
      23         376 :   params.addRequiredParam<MaterialPropertyName>("kappa_name",
      24             :                                                 "Name of the kappa for the given phase");
      25         376 :   params.addRequiredParam<MaterialPropertyName>(
      26             :       "dkappadgrad_etaa_name", "Name of the derivative of kappa w.r.t. the gradient of eta");
      27         376 :   params.addRequiredParam<MaterialPropertyName>(
      28             :       "d2kappadgrad_etaa_name",
      29             :       "Name of the second derivative of kappa w.r.t. the gradient of eta");
      30         376 :   params.addParam<Real>(
      31         376 :       "anisotropy_strength", 0.04, "Strength of the anisotropy (typically < 0.05)");
      32         376 :   params.addParam<unsigned int>("mode_number", 4, "Mode number for anisotropy");
      33         376 :   params.addParam<Real>(
      34         376 :       "reference_angle", 90, "Reference angle for defining anisotropy in degrees");
      35         376 :   params.addParam<Real>("kappa_bar", 0.1125, "Average value of the interface parameter kappa");
      36         376 :   params.addRequiredCoupledVar("etaa", "Order parameter for the current phase alpha");
      37         376 :   params.addRequiredCoupledVar("etab", "Order parameter for the neighboring phase beta");
      38         188 :   return params;
      39           0 : }
      40             : 
      41         144 : InterfaceOrientationMultiphaseMaterial::InterfaceOrientationMultiphaseMaterial(
      42         144 :     const InputParameters & parameters)
      43             :   : Material(parameters),
      44         288 :     _kappa_name(getParam<MaterialPropertyName>("kappa_name")),
      45         144 :     _dkappadgrad_etaa_name(getParam<MaterialPropertyName>("dkappadgrad_etaa_name")),
      46         144 :     _d2kappadgrad_etaa_name(getParam<MaterialPropertyName>("d2kappadgrad_etaa_name")),
      47         288 :     _delta(getParam<Real>("anisotropy_strength")),
      48         288 :     _j(getParam<unsigned int>("mode_number")),
      49         288 :     _theta0(getParam<Real>("reference_angle")),
      50         288 :     _kappa_bar(getParam<Real>("kappa_bar")),
      51         144 :     _kappa(declareProperty<Real>(_kappa_name)),
      52         144 :     _dkappadgrad_etaa(declareProperty<RealGradient>(_dkappadgrad_etaa_name)),
      53         144 :     _d2kappadgrad_etaa(declareProperty<RealTensorValue>(_d2kappadgrad_etaa_name)),
      54         144 :     _etaa(coupledValue("etaa")),
      55         144 :     _grad_etaa(coupledGradient("etaa")),
      56         144 :     _etab(coupledValue("etab")),
      57         288 :     _grad_etab(coupledGradient("etab"))
      58             : {
      59             :   // this currently only works in 2D simulations
      60         144 :   if (_mesh.dimension() != 2)
      61           0 :     mooseError("InterfaceOrientationMultiphaseMaterial requires a two-dimensional mesh.");
      62         144 : }
      63             : 
      64             : void
      65     1081800 : InterfaceOrientationMultiphaseMaterial::computeQpProperties()
      66             : {
      67             :   const Real tol = 1e-9;
      68             :   const Real cutoff = 1.0 - tol;
      69             : 
      70             :   // Normal direction of the interface
      71             :   Real n = 0.0;
      72     1081800 :   const RealGradient nd = _grad_etaa[_qp] - _grad_etab[_qp];
      73             :   const Real nx = nd(0);
      74             :   const Real ny = nd(1);
      75     1081800 :   const Real n2x = nd(0) * nd(0);
      76     1081800 :   const Real n2y = nd(1) * nd(1);
      77             :   const Real nsq = nd.norm_sq();
      78     1081800 :   const Real n2sq = nsq * nsq;
      79     1081800 :   if (nsq > tol)
      80      873132 :     n = nx / std::sqrt(nsq);
      81             : 
      82      873132 :   if (n > cutoff)
      83             :     n = cutoff;
      84             : 
      85     1050640 :   if (n < -cutoff)
      86             :     n = -cutoff;
      87             : 
      88             :   // Calculate the orientation angle
      89     1081800 :   const Real angle = std::acos(n) * MathUtils::sign(ny);
      90             : 
      91             :   // Compute derivatives of the angle wrt n
      92     1081800 :   const Real dangledn = -MathUtils::sign(ny) / std::sqrt(1.0 - n * n);
      93     1081800 :   const Real d2angledn2 = -MathUtils::sign(ny) * n / (1.0 - n * n) / std::sqrt(1.0 - n * n);
      94             : 
      95             :   // Compute derivative of n wrt grad_eta
      96             :   RealGradient dndgrad_etaa;
      97     1081800 :   if (nsq > tol)
      98             :   {
      99      873132 :     dndgrad_etaa(0) = ny * ny;
     100      873132 :     dndgrad_etaa(1) = -nx * ny;
     101      873132 :     dndgrad_etaa /= nsq * std::sqrt(nsq);
     102             :   }
     103             : 
     104             :   // Compute the square of dndgrad_etaa
     105             :   RealTensorValue dndgrad_etaa_sq;
     106     1081800 :   if (nsq > tol)
     107             :   {
     108      873132 :     dndgrad_etaa_sq(0, 0) = n2y * n2y;
     109      873132 :     dndgrad_etaa_sq(0, 1) = -nx * ny * n2y;
     110      873132 :     dndgrad_etaa_sq(1, 0) = -nx * ny * n2y;
     111      873132 :     dndgrad_etaa_sq(1, 1) = n2x * n2y;
     112      873132 :     dndgrad_etaa_sq /= n2sq * nsq;
     113             :   }
     114             : 
     115             :   // Compute the second derivative of n wrt grad_eta
     116             :   RealTensorValue d2ndgrad_etaa2;
     117     1081800 :   if (nsq > tol)
     118             :   {
     119      873132 :     d2ndgrad_etaa2(0, 0) = -3.0 * nx * n2y;
     120      873132 :     d2ndgrad_etaa2(0, 1) = -ny * n2y + 2.0 * ny * n2x;
     121      873132 :     d2ndgrad_etaa2(1, 0) = -ny * n2y + 2.0 * ny * n2x;
     122      873132 :     d2ndgrad_etaa2(1, 1) = -nx * n2x + 2.0 * nx * n2y;
     123      873132 :     d2ndgrad_etaa2 /= n2sq * std::sqrt(nsq);
     124             :   }
     125             : 
     126             :   // Calculate interfacial coefficient kappa and its derivatives wrt the angle
     127     1081800 :   Real anglediff = _j * (angle - _theta0 * libMesh::pi / 180.0);
     128     1081800 :   _kappa[_qp] =
     129     1081800 :       _kappa_bar * (1.0 + _delta * std::cos(anglediff)) * (1.0 + _delta * std::cos(anglediff));
     130     1081800 :   Real dkappadangle =
     131     1081800 :       -2.0 * _kappa_bar * _delta * _j * (1.0 + _delta * std::cos(anglediff)) * std::sin(anglediff);
     132     1081800 :   Real d2kappadangle =
     133     1081800 :       2.0 * _kappa_bar * _delta * _delta * _j * _j * std::sin(anglediff) * std::sin(anglediff) -
     134     1081800 :       2.0 * _kappa_bar * _delta * _j * _j * (1.0 + _delta * std::cos(anglediff)) *
     135             :           std::cos(anglediff);
     136             : 
     137             :   // Compute derivatives of kappa wrt grad_eta
     138     1081800 :   _dkappadgrad_etaa[_qp] = dkappadangle * dangledn * dndgrad_etaa;
     139     1081800 :   _d2kappadgrad_etaa[_qp] = d2kappadangle * dangledn * dangledn * dndgrad_etaa_sq +
     140     1081800 :                             dkappadangle * d2angledn2 * dndgrad_etaa_sq +
     141     1081800 :                             dkappadangle * dangledn * d2ndgrad_etaa2;
     142     1081800 : }