From 36bfbcadcc67959ed48e88d6ef37ed3d1bcf8a97 Mon Sep 17 00:00:00 2001 From: Mahmud Samdani Date: Mon, 3 Nov 2025 14:38:29 -0600 Subject: [PATCH] Descriptive commit message --- include/base/FerretConfig.h | 2 + .../TPElectrostrictiveCouplingPDerivative.h | 57 +++++++ .../materials/TPComputeFerroelectricStrain.h | 55 +++++++ .../TPElectrostrictiveCouplingPEnergy.h | 58 +++++++ moose | 1 + .../TPElectrostrictiveCouplingPDerivative.C | 153 ++++++++++++++++++ src/materials/TPComputeFerroelectricStrain.C | 78 +++++++++ .../TPElectrostrictiveCouplingPEnergy.C | 80 +++++++++ 8 files changed, 484 insertions(+) create mode 100644 include/kernels/TPElectrostrictiveCouplingPDerivative.h create mode 100644 include/materials/TPComputeFerroelectricStrain.h create mode 100644 include/postprocessors/TPElectrostrictiveCouplingPEnergy.h create mode 160000 moose create mode 100644 src/kernels/TPElectrostrictiveCouplingPDerivative.C create mode 100644 src/materials/TPComputeFerroelectricStrain.C create mode 100644 src/postprocessors/TPElectrostrictiveCouplingPEnergy.C diff --git a/include/base/FerretConfig.h b/include/base/FerretConfig.h index cbedbebda..619c7f38e 100644 --- a/include/base/FerretConfig.h +++ b/include/base/FerretConfig.h @@ -1,4 +1,6 @@ /* include/base/FerretConfig.h. Generated from FerretConfig.h by configure. */ +/* include/base/FerretConfig.h. Generated from FerretConfig.h by configure. */ +/* include/base/FerretConfig.h. Generated from FerretConfig.h by configure. */ /* include/base/FerretConfig.h.in. Generated from configure.ac by autoheader. */ /* Flag indicating that Boost special functions are available */ diff --git a/include/kernels/TPElectrostrictiveCouplingPDerivative.h b/include/kernels/TPElectrostrictiveCouplingPDerivative.h new file mode 100644 index 000000000..1e5459ef1 --- /dev/null +++ b/include/kernels/TPElectrostrictiveCouplingPDerivative.h @@ -0,0 +1,57 @@ +/* + This file is part of FERRET, an add-on module for MOOSE + FERRET is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + You should have received a copy of the GNU General Public License + along with this program. If not, see . + For help with FERRET please contact J. Mangeri + and be sure to track new changes at github.com/mangerij/ferret +**/ + +#ifndef TPELECTROSTRICTIVECOUPLINGPDERIVATIVE_H +#define TPELECTROSTRICTIVECOUPLINGPDERIVATIVE_H + +#include "Kernel.h" + +class TPElectrostrictiveCouplingPDerivative: public Kernel +{ +public: + TPElectrostrictiveCouplingPDerivative(const InputParameters & parameters); + + static InputParameters validParams(); + +protected: + virtual Real computeQpResidual(); + virtual Real computeQpJacobian(); + virtual Real computeQpOffDiagJacobian(unsigned int jvar); + +private: + const unsigned int _component; + const unsigned int _polar_x_var; + const unsigned int _polar_y_var; + const unsigned int _polar_z_var; + const VariableValue & _polar_x; + const VariableValue & _polar_y; + const VariableValue & _polar_z; + const MaterialProperty & _C11; + const MaterialProperty & _C12; + const MaterialProperty & _C13; + const MaterialProperty & _C33; + const MaterialProperty & _C44; + const MaterialProperty & _C66; + const MaterialProperty & _Q11; + const MaterialProperty & _Q12; + const MaterialProperty & _Q13; + const MaterialProperty & _Q33; + const MaterialProperty & _Q44; + const MaterialProperty & _Q66; + const std::string _base_name; + const MaterialProperty & _strain; +}; +#endif //TPELECTROSTRICTIVECOUPLINGPDERIVATIVE_H diff --git a/include/materials/TPComputeFerroelectricStrain.h b/include/materials/TPComputeFerroelectricStrain.h new file mode 100644 index 000000000..80f33d480 --- /dev/null +++ b/include/materials/TPComputeFerroelectricStrain.h @@ -0,0 +1,55 @@ +/* + This file is part of FERRET, an add-on module for MOOSE + + FERRET is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program. If not, see . + + For help with FERRET please contact J. Mangeri + and be sure to track new changes at github.com/mangerij/ferret + +**/ + +#ifndef TPCOMPUTEFERROELECTRICSTRAIN_H +#define TPCOMPUTEFERROELECTRICSTRAIN_H + +#include "Material.h" +#include "RankTwoTensor.h" +#include "TPComputeFerroelectricStrain.h" +#include "ComputeEigenstrainBase.h" + +/** + * TPComputeFerroelectricStrain the base class for computing spontaneous polar strain contributions (Tetragonal) + */ +class TPComputeFerroelectricStrain : public ComputeEigenstrainBase +{ +public: + TPComputeFerroelectricStrain(const InputParameters & parameters); + + static InputParameters validParams(); + void computeQpEigenstrain(); +private: + const VariableValue & _polar_x; + const VariableValue & _polar_y; + const VariableValue & _polar_z; + + const MaterialProperty & _Q11; + const MaterialProperty & _Q12; + const MaterialProperty & _Q13; + const MaterialProperty & _Q33; + const MaterialProperty & _Q44; + const MaterialProperty & _Q66; + std::vector _vals; + RankTwoTensor _polar_strain; +}; + +#endif //COMPUTEFERROELECTRICSTRAIN_H diff --git a/include/postprocessors/TPElectrostrictiveCouplingPEnergy.h b/include/postprocessors/TPElectrostrictiveCouplingPEnergy.h new file mode 100644 index 000000000..61c9c2c62 --- /dev/null +++ b/include/postprocessors/TPElectrostrictiveCouplingPEnergy.h @@ -0,0 +1,58 @@ +/* + This file is part of FERRET, an add-on module for MOOSE + + FERRET is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program. If not, see . + + For help with FERRET please contact J. Mangeri + and be sure to track new changes at github.com/mangerij/ferret + +**/ + +#ifndef TPELECTROSTRICTIVECOUPLINGPENERGY_H +#define TPELECTROSTRICTIVECOUPLINGPENERGY_H + +#include "ElementIntegralPostprocessor.h" + +class TPElectrostrictiveCouplingPEnergy : public ElementIntegralPostprocessor +{ +public: + TPElectrostrictiveCouplingPEnergy(const InputParameters & parameters); + + static InputParameters validParams(); + +protected: + virtual Real computeQpIntegral(); + +private: + const VariableValue & _polar_x; + const VariableValue & _polar_y; + const VariableValue & _polar_z; + const MaterialProperty & _C11; + const MaterialProperty & _C12; + const MaterialProperty & _C13; + const MaterialProperty & _C33; + const MaterialProperty & _C44; + const MaterialProperty & _C66; + const MaterialProperty & _Q11; + const MaterialProperty & _Q12; + const MaterialProperty & _Q13; + const MaterialProperty & _Q33; + const MaterialProperty & _Q44; + const MaterialProperty & _Q66; + const Real _energy_scale; + const std::string _base_name; + const MaterialProperty & _strain; +}; + +#endif diff --git a/moose b/moose new file mode 160000 index 000000000..77e8c30ba --- /dev/null +++ b/moose @@ -0,0 +1 @@ +Subproject commit 77e8c30ba61dcd5caba613bc1c4ecc602925737f diff --git a/src/kernels/TPElectrostrictiveCouplingPDerivative.C b/src/kernels/TPElectrostrictiveCouplingPDerivative.C new file mode 100644 index 000000000..7c09bdb9a --- /dev/null +++ b/src/kernels/TPElectrostrictiveCouplingPDerivative.C @@ -0,0 +1,153 @@ +/* + This file is part of FERRET, an add-on module for MOOSE + + FERRET is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a co_polar_y[_qp] of the GNU General Public License + along with this program. If not, see . + + For help with FERRET please contact J. Mangeri + and be sure to track new changes at github.com/mangerij/ferret + +**/ + +#include "TPElectrostrictiveCouplingPDerivative.h" + +class TPElectrostrictiveCouplingPDerivative; + +registerMooseObject("FerretApp", TPElectrostrictiveCouplingPDerivative); + +InputParameters TPElectrostrictiveCouplingPDerivative::validParams() +{ + InputParameters params = Kernel::validParams(); + params.addClassDescription("Calculates a residual contribution due to the variation w.r.t polarization of the electrostrictive coupling energy. Note: for cubic parent phase only."); + params.addRequiredCoupledVar("polar_x", "The x component of the polarization"); + params.addRequiredCoupledVar("polar_y", "The y component of the polarization"); + params.addCoupledVar("polar_z", 0.0, "The z component of the polarization"); + params.addParam("base_name", "Material property base name"); + params.addRequiredParam("component", "An integer corresponding to the direction in order parameter space this kernel acts in (e.g. for unrotated functionals 0 for q_x, 1 for q_y, 2 for q_z)."); + return params; +} + +TPElectrostrictiveCouplingPDerivative::TPElectrostrictiveCouplingPDerivative(const InputParameters & parameters) + :Kernel(parameters), + _component(getParam("component")), + _polar_x_var(coupled("polar_x")), + _polar_y_var(coupled("polar_y")), + _polar_z_var(coupled("polar_z")), + _polar_x(coupledValue("polar_x")), + _polar_y(coupledValue("polar_y")), + _polar_z(coupledValue("polar_z")), + _C11(getMaterialProperty("C11")), + _C12(getMaterialProperty("C12")), + _C13(getMaterialProperty("C13")), + _C33(getMaterialProperty("C33")), + _C44(getMaterialProperty("C44")), + _C66(getMaterialProperty("C66")), + _Q11(getMaterialProperty("Q11")), + _Q12(getMaterialProperty("Q12")), + _Q13(getMaterialProperty("Q13")), + _Q33(getMaterialProperty("Q33")), + _Q44(getMaterialProperty("Q44")), + _Q66(getMaterialProperty("Q66")), + _strain(getMaterialPropertyByName(_base_name + "total_strain")) +{ +} + +Real +TPElectrostrictiveCouplingPDerivative::computeQpResidual() +{ + if (_component == 0) + { + return (2*_C11[_qp]*_polar_x[_qp]*_Q13[_qp]*(-_strain[_qp](1,1) - _strain[_qp](2,2) + (Utility::pow<2>(_polar_y[_qp]) + Utility::pow<2>(_polar_z[_qp]))*(_Q11[_qp] + _Q12[_qp]) + 2*Utility::pow<2>(_polar_x[_qp])*_Q13[_qp]) + 2*_C13[_qp]*_polar_x[_qp]*(-2*_strain[_qp](0,0)*_Q13[_qp] + 2*(Utility::pow<2>(_polar_y[_qp]) + Utility::pow<2>(_polar_z[_qp]))*Utility::pow<2>(_Q13[_qp]) - (_strain[_qp](1,1) + _strain[_qp](2,2))*_Q13[_qp] + (Utility::pow<2>(_polar_y[_qp]) + Utility::pow<2>(_polar_z[_qp]))*(_Q11[_qp] + _Q12[_qp])*_Q13[_qp] + 4*Utility::pow<2>(_polar_x[_qp])*_Q13[_qp]*_Q13[_qp]) + 2*_polar_x[_qp]*(_C12[_qp]*_Q13[_qp]*(-_strain[_qp](1,1) - _strain[_qp](2,2) + (Utility::pow<2>(_polar_y[_qp]) + Utility::pow<2>(_polar_z[_qp]))*(_Q11[_qp] + _Q12[_qp]) + 2*Utility::pow<2>(_polar_x[_qp])*_Q13[_qp]) + _C33[_qp]*_Q13[_qp]*(-_strain[_qp](0,0) + (Utility::pow<2>(_polar_y[_qp]) + Utility::pow<2>(_polar_z[_qp]))*_Q13[_qp] + Utility::pow<2>(_polar_x[_qp])*_Q13[_qp])) - _C44[_qp]*(_strain[_qp](0,1)*_polar_y[_qp] + _strain[_qp](0,2)*_polar_z[_qp])*_Q44[_qp] + _C44[_qp]*_polar_x[_qp]*(Utility::pow<2>(_polar_y[_qp]) + Utility::pow<2>(_polar_z[_qp]))*Utility::pow<2>(_Q44[_qp]))*_test[_i][_qp]; + } + else if (_component == 1) + { + return (-2*_C12[_qp]*_strain[_qp](2,2)*_polar_y[_qp]*_Q11[_qp] + 2*_C12[_qp]*_polar_y[_qp]*Utility::pow<2>(_polar_z[_qp])*Utility::pow<2>(_Q11[_qp]) - 2*_C12[_qp]*_strain[_qp](1,1)*_polar_y[_qp]*_Q12[_qp] + 4*_C12[_qp]*Utility::pow<3>(_polar_y[_qp])*_Q11[_qp]*_Q12[_qp] + 2*_C12[_qp]*_polar_y[_qp]*Utility::pow<2>(_polar_z[_qp])*Utility::pow<2>(_Q12[_qp]) - 2*_C33[_qp]*_strain[_qp](0,0)*_polar_y[_qp]*_Q13[_qp] + 2*_C12[_qp]*Utility::pow<2>(_polar_x[_qp])*_polar_y[_qp]*_Q11[_qp]*_Q13[_qp] + 2*_C12[_qp]*Utility::pow<2>(_polar_x[_qp])*_polar_y[_qp]*_Q12[_qp]*_Q13[_qp] + 2*_C33[_qp]*Utility::pow<3>(_polar_y[_qp])*Utility::pow<2>(_Q13[_qp]) + 2*_C33[_qp]*_polar_y[_qp]*Utility::pow<2>(_polar_z[_qp])*Utility::pow<2>(_Q13[_qp]) + 2*_C11[_qp]*_polar_y[_qp]*(-(_strain[_qp](1,1)*_Q11[_qp]) - _strain[_qp](2,2)*_Q12[_qp] + 2*Utility::pow<2>(_polar_z[_qp])*_Q11[_qp]*_Q12[_qp] + Utility::pow<2>(_polar_y[_qp])*(Utility::pow<2>(_Q11[_qp]) + Utility::pow<2>(_Q12[_qp])) + Utility::pow<2>(_polar_x[_qp])*(_Q11[_qp] + _Q12[_qp])*_Q13[_qp]) + 2*_C33[_qp]*Utility::pow<2>(_polar_x[_qp])*_polar_y[_qp]*_Q13[_qp]*_Q13[_qp] + 2*_C13[_qp]*_polar_y[_qp]*(-(_strain[_qp](0,0)*(_Q11[_qp] + _Q12[_qp])) - (_strain[_qp](1,1) + _strain[_qp](2,2))*_Q13[_qp] + 2*_Q13[_qp]*((Utility::pow<2>(_polar_y[_qp]) + Utility::pow<2>(_polar_z[_qp]))*(_Q11[_qp] + _Q12[_qp]) + Utility::pow<2>(_polar_x[_qp])*_Q13[_qp]) + Utility::pow<2>(_polar_x[_qp])*(_Q11[_qp] + _Q12[_qp])*_Q13[_qp]) - _C44[_qp]*_strain[_qp](0,1)*_polar_x[_qp]*_Q44[_qp] + _C44[_qp]*Utility::pow<2>(_polar_x[_qp])*_polar_y[_qp]*Utility::pow<2>(_Q44[_qp]) - _C66[_qp]*_strain[_qp](1,2)*_polar_z[_qp]*_Q66[_qp] + _C66[_qp]*_polar_y[_qp]*Utility::pow<2>(_polar_z[_qp])*Utility::pow<2>(_Q66[_qp]))*_test[_i][_qp]; + } + else if (_component == 2) + { + return (-2*_C11[_qp]*_strain[_qp](2,2)*_polar_z[_qp]*_Q11[_qp] + 2*_C11[_qp]*Utility::pow<3>(_polar_z[_qp])*Utility::pow<2>(_Q11[_qp]) - 2*_C11[_qp]*_strain[_qp](1,1)*_polar_z[_qp]*_Q12[_qp] + 4*_C11[_qp]*Utility::pow<2>(_polar_y[_qp])*_polar_z[_qp]*_Q11[_qp]*_Q12[_qp] + 2*_C11[_qp]*Utility::pow<3>(_polar_z[_qp])*Utility::pow<2>(_Q12[_qp]) - 2*_C33[_qp]*_strain[_qp](0,0)*_polar_z[_qp]*_Q13[_qp] + 2*_C11[_qp]*Utility::pow<2>(_polar_x[_qp])*_polar_z[_qp]*_Q11[_qp]*_Q13[_qp] + 2*_C11[_qp]*Utility::pow<2>(_polar_x[_qp])*_polar_z[_qp]*_Q12[_qp]*_Q13[_qp] + 2*_C33[_qp]*Utility::pow<2>(_polar_y[_qp])*_polar_z[_qp]*Utility::pow<2>(_Q13[_qp]) + 2*_C33[_qp]*Utility::pow<3>(_polar_z[_qp])*Utility::pow<2>(_Q13[_qp]) + 2*_C12[_qp]*_polar_z[_qp]*(-(_strain[_qp](1,1)*_Q11[_qp]) - _strain[_qp](2,2)*_Q12[_qp] + 2*Utility::pow<2>(_polar_z[_qp])*_Q11[_qp]*_Q12[_qp] + Utility::pow<2>(_polar_y[_qp])*(Utility::pow<2>(_Q11[_qp]) + Utility::pow<2>(_Q12[_qp])) + Utility::pow<2>(_polar_x[_qp])*(_Q11[_qp] + _Q12[_qp])*_Q13[_qp]) + 2*_C33[_qp]*Utility::pow<2>(_polar_x[_qp])*_polar_z[_qp]*_Q13[_qp]*_Q13[_qp] + 2*_C13[_qp]*_polar_z[_qp]*(-(_strain[_qp](0,0)*(_Q11[_qp] + _Q12[_qp])) - (_strain[_qp](1,1) + _strain[_qp](2,2))*_Q13[_qp] + 2*_Q13[_qp]*((Utility::pow<2>(_polar_y[_qp]) + Utility::pow<2>(_polar_z[_qp]))*(_Q11[_qp] + _Q12[_qp]) + Utility::pow<2>(_polar_x[_qp])*_Q13[_qp]) + Utility::pow<2>(_polar_x[_qp])*(_Q11[_qp] + _Q12[_qp])*_Q13[_qp]) - _C44[_qp]*_strain[_qp](0,2)*_polar_x[_qp]*_Q44[_qp] + _C44[_qp]*Utility::pow<2>(_polar_x[_qp])*_polar_z[_qp]*Utility::pow<2>(_Q44[_qp]) - _C66[_qp]*_strain[_qp](1,2)*_polar_y[_qp]*_Q66[_qp] + _C66[_qp]*Utility::pow<2>(_polar_y[_qp])*_polar_z[_qp]*Utility::pow<2>(_Q66[_qp]))*_test[_i][_qp]; + } + else + return 0.0; +} + +Real +TPElectrostrictiveCouplingPDerivative::computeQpJacobian() +{ + if (_component == 0) + { + return _phi[_j][_qp]*(-2*_C12[_qp]*_strain[_qp](1,1)*_Q13[_qp] - 2*_C12[_qp]*_strain[_qp](2,2)*_Q13[_qp] + 2*_C12[_qp]*Utility::pow<2>(_polar_y[_qp])*_Q11[_qp]*_Q13[_qp] + 2*_C12[_qp]*Utility::pow<2>(_polar_z[_qp])*_Q11[_qp]*_Q13[_qp] + 2*_C12[_qp]*Utility::pow<2>(_polar_y[_qp])*_Q12[_qp]*_Q13[_qp] + 2*_C12[_qp]*Utility::pow<2>(_polar_z[_qp])*_Q12[_qp]*_Q13[_qp] + 12*_C12[_qp]*Utility::pow<2>(_polar_x[_qp])*Utility::pow<2>(_Q13[_qp]) + 2*_C11[_qp]*_Q13[_qp]*(-_strain[_qp](1,1) - _strain[_qp](2,2) + (Utility::pow<2>(_polar_y[_qp]) + Utility::pow<2>(_polar_z[_qp]))*(_Q11[_qp] + _Q12[_qp]) + 6*Utility::pow<2>(_polar_x[_qp])*_Q13[_qp]) - 2*_C33[_qp]*_strain[_qp](0,0)*_Q13[_qp] + 2*_C33[_qp]*Utility::pow<2>(_polar_y[_qp])*_Q13[_qp]*_Q13[_qp] + 2*_C33[_qp]*Utility::pow<2>(_polar_z[_qp])*_Q13[_qp]*_Q13[_qp] + 6*_C33[_qp]*Utility::pow<2>(_polar_x[_qp])*Utility::pow<2>(_Q13[_qp]) + 2*_C13[_qp]*(-2*_strain[_qp](0,0)*_Q13[_qp] + 2*(Utility::pow<2>(_polar_y[_qp]) + Utility::pow<2>(_polar_z[_qp]))*Utility::pow<2>(_Q13[_qp]) - (_strain[_qp](1,1) + _strain[_qp](2,2))*_Q13[_qp] + (Utility::pow<2>(_polar_y[_qp]) + Utility::pow<2>(_polar_z[_qp]))*(_Q11[_qp] + _Q12[_qp])*_Q13[_qp] + 12*Utility::pow<2>(_polar_x[_qp])*_Q13[_qp]*_Q13[_qp]) + _C44[_qp]*Utility::pow<2>(_polar_y[_qp])*Utility::pow<2>(_Q44[_qp]) + _C44[_qp]*Utility::pow<2>(_polar_z[_qp])*Utility::pow<2>(_Q44[_qp]))*_test[_i][_qp]; + } + else if (_component == 1) + { + return _phi[_j][_qp]*(-2*_C12[_qp]*_strain[_qp](2,2)*_Q11[_qp] + 2*_C12[_qp]*Utility::pow<2>(_polar_z[_qp])*Utility::pow<2>(_Q11[_qp]) - 2*_C12[_qp]*_strain[_qp](1,1)*_Q12[_qp] + 12*_C12[_qp]*Utility::pow<2>(_polar_y[_qp])*_Q11[_qp]*_Q12[_qp] + 2*_C12[_qp]*Utility::pow<2>(_polar_z[_qp])*Utility::pow<2>(_Q12[_qp]) - 2*_C33[_qp]*_strain[_qp](0,0)*_Q13[_qp] + 2*_C12[_qp]*Utility::pow<2>(_polar_x[_qp])*_Q11[_qp]*_Q13[_qp] + 2*_C12[_qp]*Utility::pow<2>(_polar_x[_qp])*_Q12[_qp]*_Q13[_qp] + 6*_C33[_qp]*Utility::pow<2>(_polar_y[_qp])*Utility::pow<2>(_Q13[_qp]) + 2*_C33[_qp]*Utility::pow<2>(_polar_z[_qp])*Utility::pow<2>(_Q13[_qp]) + 2*_C11[_qp]*(-(_strain[_qp](1,1)*_Q11[_qp]) - _strain[_qp](2,2)*_Q12[_qp] + 2*Utility::pow<2>(_polar_z[_qp])*_Q11[_qp]*_Q12[_qp] + 3*Utility::pow<2>(_polar_y[_qp])*(Utility::pow<2>(_Q11[_qp]) + Utility::pow<2>(_Q12[_qp])) + Utility::pow<2>(_polar_x[_qp])*(_Q11[_qp] + _Q12[_qp])*_Q13[_qp]) + 2*_C33[_qp]*Utility::pow<2>(_polar_x[_qp])*_Q13[_qp]*_Q13[_qp] + 2*_C13[_qp]*(-(_strain[_qp](0,0)*(_Q11[_qp] + _Q12[_qp])) - (_strain[_qp](1,1) + _strain[_qp](2,2))*_Q13[_qp] + 2*_Q13[_qp]*((3*Utility::pow<2>(_polar_y[_qp]) + Utility::pow<2>(_polar_z[_qp]))*(_Q11[_qp] + _Q12[_qp]) + Utility::pow<2>(_polar_x[_qp])*_Q13[_qp]) + Utility::pow<2>(_polar_x[_qp])*(_Q11[_qp] + _Q12[_qp])*_Q13[_qp]) + _C44[_qp]*Utility::pow<2>(_polar_x[_qp])*Utility::pow<2>(_Q44[_qp]) + _C66[_qp]*Utility::pow<2>(_polar_z[_qp])*Utility::pow<2>(_Q66[_qp]))*_test[_i][_qp]; + } + else if (_component == 2) + { + return _phi[_j][_qp]*(-2*_C12[_qp]*_strain[_qp](2,2)*_Q11[_qp] + 2*_C12[_qp]*Utility::pow<2>(_polar_z[_qp])*Utility::pow<2>(_Q11[_qp]) - 2*_C12[_qp]*_strain[_qp](1,1)*_Q12[_qp] + 12*_C12[_qp]*Utility::pow<2>(_polar_y[_qp])*_Q11[_qp]*_Q12[_qp] + 2*_C12[_qp]*Utility::pow<2>(_polar_z[_qp])*Utility::pow<2>(_Q12[_qp]) - 2*_C33[_qp]*_strain[_qp](0,0)*_Q13[_qp] + 2*_C12[_qp]*Utility::pow<2>(_polar_x[_qp])*_Q11[_qp]*_Q13[_qp] + 2*_C12[_qp]*Utility::pow<2>(_polar_x[_qp])*_Q12[_qp]*_Q13[_qp] + 6*_C33[_qp]*Utility::pow<2>(_polar_y[_qp])*Utility::pow<2>(_Q13[_qp]) + 2*_C33[_qp]*Utility::pow<2>(_polar_z[_qp])*Utility::pow<2>(_Q13[_qp]) + 2*_C11[_qp]*(-(_strain[_qp](1,1)*_Q11[_qp]) - _strain[_qp](2,2)*_Q12[_qp] + 2*Utility::pow<2>(_polar_z[_qp])*_Q11[_qp]*_Q12[_qp] + 3*Utility::pow<2>(_polar_y[_qp])*(Utility::pow<2>(_Q11[_qp]) + Utility::pow<2>(_Q12[_qp])) + Utility::pow<2>(_polar_x[_qp])*(_Q11[_qp] + _Q12[_qp])*_Q13[_qp]) + 2*_C33[_qp]*Utility::pow<2>(_polar_x[_qp])*_Q13[_qp]*_Q13[_qp] + 2*_C13[_qp]*(-(_strain[_qp](0,0)*(_Q11[_qp] + _Q12[_qp])) - (_strain[_qp](1,1) + _strain[_qp](2,2))*_Q13[_qp] + 2*_Q13[_qp]*((3*Utility::pow<2>(_polar_y[_qp]) + Utility::pow<2>(_polar_z[_qp]))*(_Q11[_qp] + _Q12[_qp]) + Utility::pow<2>(_polar_x[_qp])*_Q13[_qp]) + Utility::pow<2>(_polar_x[_qp])*(_Q11[_qp] + _Q12[_qp])*_Q13[_qp]) + _C44[_qp]*Utility::pow<2>(_polar_x[_qp])*Utility::pow<2>(_Q44[_qp]) + _C66[_qp]*Utility::pow<2>(_polar_z[_qp])*Utility::pow<2>(_Q66[_qp]))*_test[_i][_qp]; + } + else + return 0.0; +} + +Real +TPElectrostrictiveCouplingPDerivative::computeQpOffDiagJacobian(unsigned int jvar) +{ + if (_component == 0) + { + if (jvar == _polar_y_var) + { + return _phi[_j][_qp]*(4*_polar_x[_qp]*_polar_y[_qp]*((_C11[_qp] + _C12[_qp])*(_Q11[_qp] + _Q12[_qp])*_Q13[_qp] + 2*_C13[_qp]*Utility::pow<2>(_Q13[_qp]) + _C13[_qp]*(_Q11[_qp] + _Q12[_qp])*_Q13[_qp] + _C33[_qp]*_Q13[_qp]*_Q13[_qp]) - _C44[_qp]*_strain[_qp](0,1)*_Q44[_qp] + 2*_C44[_qp]*_polar_x[_qp]*_polar_y[_qp]*Utility::pow<2>(_Q44[_qp]))*_test[_i][_qp]; + } + else if (jvar == _polar_z_var) + { + return _phi[_j][_qp]*(4*_polar_x[_qp]*_polar_z[_qp]*((_C11[_qp] + _C12[_qp])*(_Q11[_qp] + _Q12[_qp])*_Q13[_qp] + 2*_C13[_qp]*Utility::pow<2>(_Q13[_qp]) + _C13[_qp]*(_Q11[_qp] + _Q12[_qp])*_Q13[_qp] + _C33[_qp]*_Q13[_qp]*_Q13[_qp]) - _C44[_qp]*_strain[_qp](0,2)*_Q44[_qp] + 2*_C44[_qp]*_polar_x[_qp]*_polar_z[_qp]*Utility::pow<2>(_Q44[_qp]))*_test[_i][_qp]; + } + else + { + return 0.0; + } + } + else if (_component == 1) + { + if (jvar == _polar_x_var) + { + return _phi[_j][_qp]*(4*_polar_x[_qp]*_polar_y[_qp]*((_C11[_qp] + _C12[_qp])*(_Q11[_qp] + _Q12[_qp])*_Q13[_qp] + 2*_C13[_qp]*Utility::pow<2>(_Q13[_qp]) + _C13[_qp]*(_Q11[_qp] + _Q12[_qp])*_Q13[_qp] + _C33[_qp]*_Q13[_qp]*_Q13[_qp]) - _C44[_qp]*_strain[_qp](0,1)*_Q44[_qp] + 2*_C44[_qp]*_polar_x[_qp]*_polar_y[_qp]*Utility::pow<2>(_Q44[_qp]))*_test[_i][_qp]; + } + else if (jvar == _polar_z_var) + { + return _phi[_j][_qp]*(-(_C66[_qp]*_strain[_qp](1,2)*_Q66[_qp]) + 2*_polar_y[_qp]*_polar_z[_qp]*(4*_C11[_qp]*_Q11[_qp]*_Q12[_qp] + 2*_C12[_qp]*(Utility::pow<2>(_Q11[_qp]) + Utility::pow<2>(_Q12[_qp])) + 2*_Q13[_qp]*(2*_C13[_qp]*(_Q11[_qp] + _Q12[_qp]) + _C33[_qp]*_Q13[_qp]) + _C66[_qp]*Utility::pow<2>(_Q66[_qp])))*_test[_i][_qp]; + } + else + { + return 0.0; + } + } + else if (_component == 2) + { + if (jvar == _polar_x_var) + { + return _phi[_j][_qp]*(4*_polar_x[_qp]*_polar_z[_qp]*((_C11[_qp] + _C12[_qp])*(_Q11[_qp] + _Q12[_qp])*_Q13[_qp] + 2*_C13[_qp]*Utility::pow<2>(_Q13[_qp]) + _C13[_qp]*(_Q11[_qp] + _Q12[_qp])*_Q13[_qp] + _C33[_qp]*_Q13[_qp]*_Q13[_qp]) - _C44[_qp]*_strain[_qp](0,2)*_Q44[_qp] + 2*_C44[_qp]*_polar_x[_qp]*_polar_z[_qp]*Utility::pow<2>(_Q44[_qp]))*_test[_i][_qp]; + } + else if (jvar == _polar_y_var) + { + return _phi[_j][_qp]*(-(_C66[_qp]*_strain[_qp](1,2)*_Q66[_qp]) + 2*_polar_y[_qp]*_polar_z[_qp]*(4*_C11[_qp]*_Q11[_qp]*_Q12[_qp] + 2*_C12[_qp]*(Utility::pow<2>(_Q11[_qp]) + Utility::pow<2>(_Q12[_qp])) + 2*_Q13[_qp]*(2*_C13[_qp]*(_Q11[_qp] + _Q12[_qp]) + _C33[_qp]*_Q13[_qp]) + _C66[_qp]*Utility::pow<2>(_Q66[_qp])))*_test[_i][_qp]; + } + else + { + return 0.0; + } + } + else + return 0.0; +} diff --git a/src/materials/TPComputeFerroelectricStrain.C b/src/materials/TPComputeFerroelectricStrain.C new file mode 100644 index 000000000..5642a077e --- /dev/null +++ b/src/materials/TPComputeFerroelectricStrain.C @@ -0,0 +1,78 @@ +/* + This file is part of FERRET, an add-on module for MOOSE + + FERRET is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program. If not, see . + + For help with FERRET please contact J. Mangeri + and be sure to track new changes at github.com/mangerij/ferret + +**/ + +#include "TPComputeFerroelectricStrain.h" +#include "RankTwoTensor.h" + +registerMooseObject("FerretApp", TPComputeFerroelectricStrain); + +InputParameters TPComputeFerroelectricStrain::validParams() +{ + InputParameters params = ComputeEigenstrainBase::validParams(); + params.addClassDescription("Compute the spontaneous polar contribution to the strain."); + params.addRequiredCoupledVar("polar_x", "The x component of the polarization"); + params.addCoupledVar("polar_y", 0.0, "The y component of the polarization"); + params.addCoupledVar("polar_z", 0.0, "The z component of the polarization"); + return params; +} + +TPComputeFerroelectricStrain::TPComputeFerroelectricStrain(const InputParameters & parameters) : + ComputeEigenstrainBase(parameters), + _polar_x(coupledValue("polar_x")), + _polar_y(coupledValue("polar_y")), + _polar_z(coupledValue("polar_z")), + _Q11(getMaterialProperty("Q11")), + _Q12(getMaterialProperty("Q12")), + _Q13(getMaterialProperty("Q13")), + _Q33(getMaterialProperty("Q33")), + _Q44(getMaterialProperty("Q44")), + _Q66(getMaterialProperty("Q66")), + _vals(6), + _polar_strain() +{ +} + +void +TPComputeFerroelectricStrain::computeQpEigenstrain() +{ + RealVectorValue w(_polar_x[_qp], _polar_y[_qp], _polar_z[_qp]); + + // _vals[0] = _Q11[_qp]*w(0)*w(0) + _Q12[_qp]*(w(1)*w(1)+w(2)*w(2)); // eps_1 + _vals[0] = _Q33[_qp]*w(0)*w(0) + _Q13[_qp]*(w(1)*w(1)+w(2)*w(2)); // eps_1 = Q33*(Px^2) + Q13*(Py^2) + Q13*(Pz^2) + + // _vals[1] = _Q11[_qp]*w(1)*w(1) + _Q12[_qp]*(w(0)*w(0)+w(2)*w(2)); // eps_2 + _vals[1] = _Q13[_qp]*w(0)*w(0) + _Q11[_qp]*w(1)*w(1) + _Q12[_qp]*w(2)*w(2); // eps_2 = Q13*(Px^2) + Q11(Py^2) + Q12*(Pz^2) + + // _vals[2] = _Q11[_qp]*w(2)*w(2) + _Q12[_qp]*(w(0)*w(0)+w(1)*w(1)); // eps_3 + _vals[2] = _Q13[_qp]*w(0)*w(0) + _Q12[_qp]*w(1)*w(1) + _Q11[_qp]*w(2)*w(2); // eps_3 = Q13*(Px^2) + Q12*(Py^2) + Q11*(Pz^2) + + // _vals[3] = _Q44[_qp]*(w(1)*w(2)); //eps_4 23 + _vals[3] = _Q66[_qp]*(w(1)*w(2)); //eps_4 23 = Q66*Py*Pz + + // _vals[4] = _Q44[_qp]*(w(0)*w(2)); //eps_5 13 + _vals[4] = _Q44[_qp]*(w(0)*w(2)); //eps_5 13 = Q44*Px*Pz + + // _vals[5] = _Q44[_qp]*(w(0)*w(1)); //eps_6 12 + _vals[5] = _Q44[_qp]*(w(0)*w(1)); //eps_6 12 = Q44*Px*Py + + _polar_strain.fillFromInputVector(_vals); + _eigenstrain[_qp] = _polar_strain; +} diff --git a/src/postprocessors/TPElectrostrictiveCouplingPEnergy.C b/src/postprocessors/TPElectrostrictiveCouplingPEnergy.C new file mode 100644 index 000000000..ceeed07df --- /dev/null +++ b/src/postprocessors/TPElectrostrictiveCouplingPEnergy.C @@ -0,0 +1,80 @@ +/* + This file is part of FERRET, an add-on module for MOOSE + + FERRET is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a co_polar_y[_qp] of the GNU General Public License + along with this program. If not, see . + + For help with FERRET please contact J. Mangeri + and be sure to track new changes at github.com/mangerij/ferret + +**/ + +#include "TPElectrostrictiveCouplingPEnergy.h" + +registerMooseObject("FerretApp", TPElectrostrictiveCouplingPEnergy); + +InputParameters TPElectrostrictiveCouplingPEnergy::validParams() +{ + InputParameters params = ElementIntegralPostprocessor::validParams(); + params.addClassDescription("Calculates a volume integral over the electrostrictive coupling energy density."); + params.addRequiredCoupledVar("polar_x", "The x component of the polarization vector"); + params.addRequiredCoupledVar("polar_y", "The y component of the polarization vector"); + params.addCoupledVar("polar_z", 0.0, "The z component of the polarization vector"); + params.addParam("base_name", "Material property base name"); + params.addParam("energy_scale", 1.0, "the energy scale, useful for transition between eV and J"); + return params; +} + +TPElectrostrictiveCouplingPEnergy::TPElectrostrictiveCouplingPEnergy(const InputParameters & parameters) : + ElementIntegralPostprocessor(parameters), + _polar_x(coupledValue("polar_x")), + _polar_y(coupledValue("polar_y")), + _polar_z(coupledValue("polar_z")), + _C11(getMaterialProperty("C11")), + _C12(getMaterialProperty("C12")), + _C13(getMaterialProperty("C13")), + _C33(getMaterialProperty("C33")), + _C44(getMaterialProperty("C44")), + _C66(getMaterialProperty("C66")), + _Q11(getMaterialProperty("Q11")), + _Q12(getMaterialProperty("Q12")), + _Q13(getMaterialProperty("Q13")), + _Q33(getMaterialProperty("Q33")), + _Q44(getMaterialProperty("Q44")), + _Q66(getMaterialProperty("Q66")), + _energy_scale(getParam("energy_scale")), + _strain(getMaterialPropertyByName(_base_name + "total_strain")) +{ +} + +Real +TPElectrostrictiveCouplingPEnergy::computeQpIntegral() +{ + return _energy_scale*(_C11[_qp]*Utility::pow<2>(-_strain[_qp](2,2) + Utility::pow<2>(_polar_z[_qp])*_Q11[_qp] + Utility::pow<2>(_polar_y[_qp])*_Q12[_qp] + Utility::pow<2>(_polar_x[_qp])*_Q13[_qp]) + - 2*_C12[_qp]*(_strain[_qp](2,2) - Utility::pow<2>(_polar_z[_qp])*_Q11[_qp] - Utility::pow<2>(_polar_y[_qp])*_Q12[_qp] - Utility::pow<2>(_polar_x[_qp])*_Q13[_qp])* (-_strain[_qp](1,1) + Utility::pow<2>(_polar_y[_qp])*_Q11[_qp] + Utility::pow<2>(_polar_z[_qp])*_Q12[_qp] + Utility::pow<2>(_polar_x[_qp])*_Q13[_qp]) + + _C11[_qp]*Utility::pow<2>(-_strain[_qp](1,1) + Utility::pow<2>(_polar_y[_qp])*_Q11[_qp] + Utility::pow<2>(_polar_z[_qp])*_Q12[_qp] + Utility::pow<2>(_polar_x[_qp])*_Q13[_qp]) + - 2*_C13[_qp]*(_strain[_qp](2,2) - Utility::pow<2>(_polar_z[_qp])*_Q11[_qp] - Utility::pow<2>(_polar_y[_qp])*_Q12[_qp] - Utility::pow<2>(_polar_x[_qp])*_Q13[_qp])* (-_strain[_qp](0,0) + (Utility::pow<2>(_polar_y[_qp]) + Utility::pow<2>(_polar_z[_qp]))*_Q13[_qp] + Utility::pow<2>(_polar_x[_qp])*_Q33[_qp]) + - 2*_C13[_qp]*(_strain[_qp](1,1) - Utility::pow<2>(_polar_y[_qp])*_Q11[_qp] - Utility::pow<2>(_polar_z[_qp])*_Q12[_qp] - Utility::pow<2>(_polar_x[_qp])*_Q13[_qp])* (-_strain[_qp](0,0) + (Utility::pow<2>(_polar_y[_qp]) + Utility::pow<2>(_polar_z[_qp]))*_Q13[_qp] + Utility::pow<2>(_polar_x[_qp])*_Q33[_qp]) + + _C33[_qp]*Utility::pow<2>(-_strain[_qp](0,0) + (Utility::pow<2>(_polar_y[_qp]) + Utility::pow<2>(_polar_z[_qp]))*_Q13[_qp] + Utility::pow<2>(_polar_x[_qp])*_Q33[_qp]) + + _C44[_qp]*Utility::pow<2>(_strain[_qp](0,1) - _polar_x[_qp]*_polar_y[_qp]*_Q44[_qp]) + + _C44[_qp]*Utility::pow<2>(_strain[_qp](0,2) - _polar_x[_qp]*_polar_z[_qp]*_Q44[_qp]) + + _C66[_qp]*Utility::pow<2>(_strain[_qp](1,2) - _polar_y[_qp]*_polar_z[_qp]*_Q66[_qp]))/2; +} + + + + + + + +