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 "InterfaceOrientationMaterial.h"
      11             : #include "MooseMesh.h"
      12             : #include "MathUtils.h"
      13             : 
      14             : registerMooseObject("PhaseFieldApp", InterfaceOrientationMaterial);
      15             : 
      16             : InputParameters
      17          82 : InterfaceOrientationMaterial::validParams()
      18             : {
      19          82 :   InputParameters params = Material::validParams();
      20          82 :   params.addClassDescription("2D interfacial anisotropy");
      21         164 :   params.addParam<Real>(
      22         164 :       "anisotropy_strength", 0.04, "Strength of the anisotropy (typically < 0.05)");
      23         164 :   params.addParam<unsigned int>("mode_number", 6, "Mode number for anisotropy");
      24         164 :   params.addParam<Real>(
      25         164 :       "reference_angle", 90, "Reference angle for defining anisotropy in degrees");
      26         164 :   params.addParam<Real>("eps_bar", 0.01, "Average value of the interface parameter epsilon");
      27         164 :   params.addRequiredCoupledVar("op", "Order parameter defining the solid phase");
      28          82 :   return params;
      29           0 : }
      30             : 
      31          63 : InterfaceOrientationMaterial::InterfaceOrientationMaterial(const InputParameters & parameters)
      32             :   : Material(parameters),
      33          63 :     _delta(getParam<Real>("anisotropy_strength")),
      34         126 :     _j(getParam<unsigned int>("mode_number")),
      35         126 :     _theta0(getParam<Real>("reference_angle")),
      36         126 :     _eps_bar(getParam<Real>("eps_bar")),
      37          63 :     _eps(declareProperty<Real>("eps")),
      38          63 :     _deps(declareProperty<Real>("deps")),
      39          63 :     _depsdgrad_op(declareProperty<RealGradient>("depsdgrad_op")),
      40          63 :     _ddepsdgrad_op(declareProperty<RealGradient>("ddepsdgrad_op")),
      41          63 :     _op(coupledValue("op")),
      42         126 :     _grad_op(coupledGradient("op"))
      43             : {
      44             :   // this currently only works in 2D simulations
      45          63 :   if (_mesh.dimension() != 2)
      46           0 :     mooseError("InterfaceOrientationMaterial requires a two-dimensional mesh.");
      47          63 : }
      48             : 
      49             : void
      50     3383296 : InterfaceOrientationMaterial::computeQpProperties()
      51             : {
      52             :   const Real tol = 1e-9;
      53             :   const Real cutoff = 1.0 - tol;
      54             : 
      55             :   // cosine of the gradient orientation angle
      56             :   Real n = 0.0;
      57     3383296 :   const Real nsq = _grad_op[_qp].norm_sq();
      58     3383296 :   if (nsq > tol)
      59     2514944 :     n = _grad_op[_qp](0) / std::sqrt(nsq);
      60             : 
      61     2514944 :   if (n > cutoff)
      62             :     n = cutoff;
      63             : 
      64     3383296 :   if (n < -cutoff)
      65             :     n = -cutoff;
      66             : 
      67     3383296 :   const Real angle = std::acos(n) * MathUtils::sign(_grad_op[_qp](1));
      68             : 
      69             :   // Compute derivative of angle wrt n
      70     3383296 :   const Real dangledn = -MathUtils::sign(_grad_op[_qp](1)) / std::sqrt(1.0 - n * n);
      71             : 
      72             :   // Compute derivative of n with respect to grad_op
      73             :   RealGradient dndgrad_op;
      74     3383296 :   if (nsq > tol)
      75             :   {
      76     2514944 :     dndgrad_op(0) = _grad_op[_qp](1) * _grad_op[_qp](1);
      77     2514944 :     dndgrad_op(1) = -_grad_op[_qp](0) * _grad_op[_qp](1);
      78     2514944 :     dndgrad_op /= (_grad_op[_qp].norm_sq() * _grad_op[_qp].norm());
      79             :   }
      80             : 
      81             :   // Calculate interfacial parameter epsilon and its derivatives
      82     3383296 :   _eps[_qp] = _eps_bar * (_delta * std::cos(_j * (angle - _theta0 * libMesh::pi / 180.0)) + 1.0);
      83     3383296 :   _deps[_qp] = -_eps_bar * _delta * _j * std::sin(_j * (angle - _theta0 * libMesh::pi / 180.0));
      84     3383296 :   Real d2eps =
      85     3383296 :       -_eps_bar * _delta * _j * _j * std::cos(_j * (angle - _theta0 * libMesh::pi / 180.0));
      86             : 
      87             :   // Compute derivatives of epsilon and its derivative wrt grad_op
      88     3383296 :   _depsdgrad_op[_qp] = _deps[_qp] * dangledn * dndgrad_op;
      89     3383296 :   _ddepsdgrad_op[_qp] = d2eps * dangledn * dndgrad_op;
      90     3383296 : }