diff --git a/include/tensor_computes/FFTMechanics.h b/include/tensor_computes/FFTMechanics.h index 7c93b4c9..f7097d46 100644 --- a/include/tensor_computes/FFTMechanics.h +++ b/include/tensor_computes/FFTMechanics.h @@ -65,6 +65,7 @@ class FFTMechanics : public TensorOperator<> const torch::Tensor * const _applied_macroscopic_strain; const bool _verbose; + const bool _accept_nonconverged; using TensorOperatorBase::_dim; using TensorOperatorBase::_domain; diff --git a/include/tensor_computes/HyperElasticIsotropic.h b/include/tensor_computes/HyperElasticIsotropic.h index aae36a65..2923798c 100644 --- a/include/tensor_computes/HyperElasticIsotropic.h +++ b/include/tensor_computes/HyperElasticIsotropic.h @@ -31,7 +31,11 @@ class HyperElasticIsotropic : public TensorOperator<> const torch::Tensor _tI4s; const torch::Tensor _tII; + /// deformation gradient tensor const torch::Tensor & _tF; + /// optional eigenstrain tensor + const torch::Tensor * _pFstar; + const torch::Tensor & _tmu; const torch::Tensor & _tK; torch::Tensor & _tK4; diff --git a/src/tensor_computes/FFTMechanics.C b/src/tensor_computes/FFTMechanics.C index ced49666..eb2f5647 100644 --- a/src/tensor_computes/FFTMechanics.C +++ b/src/tensor_computes/FFTMechanics.C @@ -42,6 +42,7 @@ FFTMechanics::validParams() "Applied macroscopic strain"); params.addParam("F", "F", "Deformation gradient tensor."); params.addParam("verbose", false, "Print non-linear residuals."); + params.addParam("accept_nonconverged", false, "Accept a non-converged solution"); return params; } @@ -69,7 +70,8 @@ FFTMechanics::FFTMechanics(const InputParameters & parameters) _applied_macroscopic_strain(isParamValid("applied_macroscopic_strain") ? &getInputBuffer("applied_macroscopic_strain") : nullptr), - _verbose(getParam("verbose")) + _verbose(getParam("verbose")), + _accept_nonconverged(getParam("accept_nonconverged")) { // Build projection tensor once const auto & q = _domain.getKGrid(); @@ -148,7 +150,7 @@ FFTMechanics::computeBuffer() // print nonlinear residual to the screen if (_verbose) - _console << "|R|=" << anorm << "\t|R/R0|=" << rnorm << '\n'; + _console << "|R|=" << anorm << "\t|R/R0|=" << rnorm << std::endl; // check convergence if ((rnorm < _nl_rel_tol || anorm < _nl_abs_tol) && iiter > 0) @@ -157,7 +159,17 @@ FFTMechanics::computeBuffer() iiter++; if (iiter > _nl_max_its) - paramError("nl_max_its", - "Exceeded the maximum number of nonlinear iterations without converging."); + { + if (_accept_nonconverged) + paramWarning( + "nl_max_its", "Accepting non-converged solution |R|=", anorm, "\t|R/R0|=", rnorm); + else + paramError( + "nl_max_its", + "Exceeded the maximum number of nonlinear iterations without converging with |R|=", + anorm, + "\t|R/R0|=", + rnorm); + } } } diff --git a/src/tensor_computes/HyperElasticIsotropic.C b/src/tensor_computes/HyperElasticIsotropic.C index 47e88e64..0e097c5a 100644 --- a/src/tensor_computes/HyperElasticIsotropic.C +++ b/src/tensor_computes/HyperElasticIsotropic.C @@ -17,6 +17,7 @@ HyperElasticIsotropic::validParams() InputParameters params = TensorOperator<>::validParams(); params.addClassDescription("Hyperelastic isotropic constitutive model."); params.addRequiredParam("F", "Deformation gradient tensor"); + params.addParam("Fstar", "Optional eigenstrain tensor"); params.addRequiredParam("mu", "Deformation gradient tensor"); params.addRequiredParam("K", "Deformation gradient tensor"); params.addParam("tangent_operator", "dstressdstrain", "Stiffness tensor"); @@ -32,6 +33,7 @@ HyperElasticIsotropic::HyperElasticIsotropic(const InputParameters & parameters) _tI4s((_tI4 + _tI4rt) / 2.0), _tII(MooseTensor::dyad22(_tI, _tI)), _tF(getInputBuffer("F")), + _pFstar(isParamValid("Fstar") ? &getInputBuffer("Fstar") : nullptr), _tmu(getInputBuffer("mu")), _tK(getInputBuffer("K")), _tK4(getOutputBuffer("tangent_operator")) @@ -43,10 +45,25 @@ HyperElasticIsotropic::computeBuffer() { using namespace MooseTensor; + const auto F_e = _pFstar ? dot22(_tF, torch::inverse(*_pFstar)) : _tF; + + // Compute Green–Lagrange strain E = 0.5 (F_e^T * F_e - I) + const auto E = 0.5 * (dot22(trans2(F_e), F_e) - _tI); + + // Build material stiffness tensor const auto C4 = _tK.reshape(_domain.getValueShape({1, 1, 1, 1})) * _tII + 2. * _tmu.reshape(_domain.getValueShape({1, 1, 1, 1})) * (_tI4s - 1. / 3. * _tII); - const auto S = ddot42(C4, .5 * (dot22(trans2(_tF), _tF) - _tI)); - _u = dot22(_tF, S); - _tK4 = dot24(S, _tI4) + ddot44(ddot44(_tI4rt, dot42(dot24(_tF, C4), trans2(_tF))), _tI4rt); + // Second Piola-Kirchhoff stress: S = C : E + const auto S = ddot42(C4, E); + + // First Piola-Kirchhoff stress: P = F_e * S + _u = dot22(F_e, S); + + _tK4 = dot24(S, _tI4) + ddot44(ddot44(_tI4rt, dot42(dot24(F_e, C4), trans2(F_e))), _tI4rt); + + // Consistent tangent: K_ijkl = F^e_{im} C_{jlmn} F^e_{kn} + delta_{ik} S_{jl} + // const auto K_geo = torch::einsum("...im,...jlmn,...kn->...ijkl", {F_e, C4, F_e}); + // const auto K_stress = torch::einsum("ik,jl,...jl->...ijkl", {_ti, _ti, S}); + // _tK4 = K_geo + K_stress; }