diff --git a/src/model/solid_mechanics/materials/material_plastic/material_drucker_prager.hh b/src/model/solid_mechanics/materials/material_plastic/material_drucker_prager.hh index 8b2eedfa4..12bc3fe63 100644 --- a/src/model/solid_mechanics/materials/material_plastic/material_drucker_prager.hh +++ b/src/model/solid_mechanics/materials/material_plastic/material_drucker_prager.hh @@ -1,140 +1,142 @@ /** * @file material_drucker_pruger.hh * * @author Mohit Pundir * * @date creation: Wed Sep 09 2020 * @date last modification: Wed Sep 09 2020 * * @brief Specialization of the material class for isotropic * plasticity with Drucker-Pruger yield criterion * * * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu 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 Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "aka_voigthelper.hh" #include "material_plastic.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MATERIAL_DRUCKER_PRAGER_HH__ #define __AKANTU_MATERIAL_DRUCKER_PRAGER_HH__ namespace akantu { /** * Material plastic with a Drucker-pruger yield criterion */ template class MaterialDruckerPrager : public MaterialPlastic { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: MaterialDruckerPrager(SolidMechanicsModel & model, const ID & id = ""); MaterialDruckerPrager(SolidMechanicsModel & model, UInt dim, const Mesh & mesh, FEEngine & fe_engine, const ID & id = ""); protected: using voigt_h = VoigtHelper; void initialize(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// constitutive law for all element of a type void computeStress(ElementType el_type, GhostType ghost_type = _not_ghost) override; /// compute the tangent stiffness matrix for an element type void computeTangentModuli(const ElementType & el_type, Array & tangent_matrix, GhostType ghost_type = _not_ghost) override; protected: /// Infinitesimal deformations inline void computeStressOnQuad( const Matrix & grad_u, const Matrix & previous_grad_u, Matrix & sigma, const Matrix & previous_sigma, Matrix & inelas_strain, const Matrix & previous_inelas_strain, const Real & sigma_th, const Real & previous_sigma_th); /// Finite deformations inline void computeStressOnQuad( const Matrix & grad_u, const Matrix & previous_grad_u, Matrix & sigma, const Matrix & previous_sigma, Matrix & inelas_strain, const Matrix & previous_inelas_strain, const Real & sigma_th, const Real & previous_sigma_th, const Matrix & F_tensor); inline void computeTangentModuliOnQuad( Matrix & tangent, const Matrix & grad_u, const Matrix & previous_grad_u, const Matrix & sigma_tensor, const Matrix & previous_sigma_tensor) const; inline Real computeYieldFunction(const Matrix & sigma); + inline Real computeYieldStress(const Matrix & sigma); + inline void computeDeviatoricStress(const Matrix & sigma, Matrix & sigma_dev); /// rcompute the alpha and k parameters void updateInternalParameters() override; public: // closet point projection method to compute stress state on the // yield surface inline void computeGradientAndPlasticMultplier( const Matrix & sigma_tr, Real & plastic_multiplier_guess, Vector & gradient_f, Vector & delta_inelastic_strain, UInt max_iterations = 100, Real tolerance = 1e-10); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: // Internal friction angle of the material Real phi; // Compressive strength of the material Real fc; // modified friction angle for Drucker-Prager Real alpha; // modified compressive strength for Drucker-Prager Real k; }; } // namespace akantu #include "material_drucker_prager_inline_impl.hh" #endif /*__AKANTU_MATERIAL_DRUCKER_PRAGER_HH__ */ diff --git a/src/model/solid_mechanics/materials/material_plastic/material_drucker_prager_inline_impl.hh b/src/model/solid_mechanics/materials/material_plastic/material_drucker_prager_inline_impl.hh index 1306d7c27..92b682222 100644 --- a/src/model/solid_mechanics/materials/material_plastic/material_drucker_prager_inline_impl.hh +++ b/src/model/solid_mechanics/materials/material_plastic/material_drucker_prager_inline_impl.hh @@ -1,443 +1,460 @@ /** * @file material_linear_isotropic_hardening_inline_impl.hh * * @author Mohit Pundir * * @date creation: Wed Sep 09 2020 * @date last modification: Wed Sep 09 2020 * * @brief Implementation of the inline functions of the material * Drucker-Prager * * * Copyright (©) 2014-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu 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 Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ #include "material_drucker_prager.hh" namespace akantu { /* -------------------------------------------------------------------------- */ /// Deviatoric Stress template inline void MaterialDruckerPrager::computeDeviatoricStress(const Matrix & sigma, Matrix & sigma_dev){ for (UInt i = 0; i < dim; ++i) for (UInt j = 0; j < dim; ++j) sigma_dev(i, j) = sigma(i, j); sigma_dev -= Matrix::eye(dim, sigma.trace() / dim); } +/* -------------------------------------------------------------------------- */ +/// Yield Stress +template +inline Real MaterialDruckerPrager::computeYieldStress(const Matrix & + sigma) { + return this->alpha * sigma.trace() - this->k; +} + + /* -------------------------------------------------------------------------- */ /// Yield function template inline Real MaterialDruckerPrager::computeYieldFunction(const Matrix & sigma) { Matrix sigma_dev(dim, dim, 0); this->computeDeviatoricStress(sigma, sigma_dev); // compute deviatoric invariant J2 Real j2 = (1./ 2.) * sigma_dev.doubleDot(sigma_dev); Real sigma_dev_eff = std::sqrt(3. * j2); - Real modified_yield_stress = this->alpha * sigma.trace() - this->k; + Real modified_yield_stress = computeYieldStress(sigma); return sigma_dev_eff + modified_yield_stress; } /* -------------------------------------------------------------------------- */ template inline void MaterialDruckerPrager::computeGradientAndPlasticMultplier( const Matrix & sigma_trial, Real & plastic_multiplier_guess, Vector & gradient_f, Vector & delta_inelastic_strain, UInt max_iterations, Real tolerance) { UInt size = voigt_h::size; // guess stress state at each iteration, initial guess is the trial state Matrix sigma_guess(sigma_trial); // plastic multiplier guess at each iteration, initial guess is zero plastic_multiplier_guess = 0.; // gradient of yield surface in voigt notation gradient_f.clear(); // plastic strain increment at each iteration delta_inelastic_strain.clear(); // variation in sigma at each iteration Vector delta_sigma(size, 0.); // krocker delta vector in voigt notation Vector kronecker_delta(size, 0.); for(auto i : arange(dim)) kronecker_delta[i] = 1.; // hessian matrix of yield surface Matrix hessian_f(size, size, 0.); // scaling matrix for computing gradient and hessian from voigt notation Matrix scaling_matrix(size, size, 0.); scaling_matrix.eye(1.); for(auto i : arange(dim, size)) for(auto j : arange(dim, size)) scaling_matrix(i, j) *= 2.; // elastic stifnness tensor Matrix De(size, size, 0.); MaterialElastic::computeTangentModuliOnQuad(De); // elastic compliance tensor Matrix Ce(size, size, 0.); Ce.inverse(De); // objective function to be computed Vector f(size, 0.); // yield function value at each iteration Real yield_function; // if sigma is above the threshold value auto above_threshold = [&sigma_guess](Real & k, Real & alpha) { Real I1 = sigma_guess.trace(); return I1 >= k/alpha; }; // to project stress state at origin of yield function if first // invariant is greater than the threshold if(above_threshold(k, alpha) and this->alpha > 0) { auto update_first_obj = [&sigma_guess]() { //const UInt dimension = sigma_guess.cols(); Matrix sigma_dev(dim, dim, 0); for (UInt i = 0; i < dim; ++i) for (UInt j = 0; j < dim; ++j) sigma_dev(i, j) = sigma_guess(i, j); sigma_dev -= Matrix::eye(dim, sigma_guess.trace() / dim); auto error = (1./2) *sigma_dev.doubleDot(sigma_dev); return error; }; auto update_sec_obj = [&sigma_guess](Real & k, Real & alpha) { auto error = alpha*sigma_guess.trace() - k; return error; }; auto projection_error = update_first_obj(); while(tolerance < projection_error) { Matrix delta_sigma(dim, dim); Matrix jacobian(dim, dim); Matrix jacobian_inv(dim, dim); Matrix sigma_dev(dim, dim, 0); this->computeDeviatoricStress(sigma_guess, sigma_dev); jacobian_inv.inverse(sigma_dev); delta_sigma = -projection_error * jacobian_inv; sigma_guess += delta_sigma; projection_error = update_first_obj(); } projection_error = update_sec_obj(k, alpha); while(tolerance < projection_error) { Matrix delta_sigma(dim, dim); Matrix jacobian(dim, dim); Matrix jacobian_inv(dim, dim); jacobian = this->alpha * Matrix::eye(dim, dim); jacobian_inv.inverse(jacobian); delta_sigma += -projection_error * jacobian_inv; sigma_guess += delta_sigma; projection_error = update_sec_obj(k, alpha); } auto delta_sigma_final = sigma_trial - sigma_guess; auto delta_sigma_voigt = voigt_h::matrixToVoigt(delta_sigma_final); delta_inelastic_strain.mul(Ce, delta_sigma_voigt); return; } // lambda function to compute gradient of yield surface in voigt notation auto compute_gradient_f = [&sigma_guess, &scaling_matrix, &kronecker_delta, &gradient_f](Real & alpha){ //const UInt dimension = sigma_guess.cols(); Matrix sigma_dev(dim, dim, 0); for (UInt i = 0; i < dim; ++i) for (UInt j = 0; j < dim; ++j) sigma_dev(i, j) = sigma_guess(i, j); sigma_dev -= Matrix::eye(dim, sigma_guess.trace() / dim); //this->computeDeviatoricStress(sigma_guess, sigma_dev); Vector sigma_dev_voigt = voigt_h::matrixToVoigt(sigma_dev); // compute deviatoric invariant Real j2 = (1./2.) * sigma_dev.doubleDot(sigma_dev); gradient_f.mul(scaling_matrix, sigma_dev_voigt, 3./ (2. * std::sqrt(3. * j2)) ); gradient_f += alpha * kronecker_delta; }; // lambda function to compute hessian matrix of yield surface auto compute_hessian_f = [&sigma_guess, &hessian_f, &scaling_matrix, &kronecker_delta](){ //const UInt dimension = sigma_guess.cols(); Matrix sigma_dev(dim, dim, 0); //this->computeDeviatoricStress(sigma_guess, sigma_dev); for (UInt i = 0; i < dim; ++i) for (UInt j = 0; j < dim; ++j) sigma_dev(i, j) = sigma_guess(i, j); sigma_dev -= Matrix::eye(dim, sigma_guess.trace() / dim); auto sigma_dev_voigt = voigt_h::matrixToVoigt(sigma_dev); // compute deviatoric invariant J2 Real j2 = (1./2.) * sigma_dev.doubleDot(sigma_dev); Vector temp(sigma_dev_voigt.size()); temp.mul(scaling_matrix, sigma_dev_voigt); Matrix id(kronecker_delta.size(), kronecker_delta.size()); id.outerProduct(kronecker_delta, kronecker_delta); id *= -1./3.; id += Matrix::eye(kronecker_delta.size(), 1.); Matrix tmp3(kronecker_delta.size(), kronecker_delta.size()); tmp3.mul(scaling_matrix, id); hessian_f.outerProduct(temp, temp); hessian_f *= -9./(4.* pow(3.*j2, 3./2.)); hessian_f += (3./(2.* pow(3.*j2, 1./2.)))*tmp3; }; /* --------------------------- */ /* init before iteration loop */ /* --------------------------- */ auto update_f = [&f, &sigma_guess, &sigma_trial, &plastic_multiplier_guess, &Ce, &De, &yield_function, &gradient_f, &delta_inelastic_strain, &compute_gradient_f](Real & k, Real & alpha){ // compute gradient compute_gradient_f(alpha); // compute yield function - //yield_function = this->computeYieldFunction(sigma_guess); - //const UInt dimension = sigma_guess.cols(); Matrix sigma_dev(dim, dim, 0); for (UInt i = 0; i < dim; ++i) for (UInt j = 0; j < dim; ++j) sigma_dev(i, j) = sigma_guess(i, j); sigma_dev -= Matrix::eye(dim, sigma_guess.trace() / dim); Real j2 = (1./ 2.) * sigma_dev.doubleDot(sigma_dev); Real sigma_dev_eff = std::sqrt(3. * j2); Real modified_yield_stress = alpha * sigma_guess.trace() - k; yield_function = sigma_dev_eff + modified_yield_stress; // compute increment strain auto sigma_trial_voigt = voigt_h::matrixToVoigt(sigma_trial); auto sigma_guess_voigt = voigt_h::matrixToVoigt(sigma_guess); auto tmp = sigma_trial_voigt - sigma_guess_voigt; delta_inelastic_strain.mul(Ce, tmp); // compute objective function f.mul(De, gradient_f, plastic_multiplier_guess); f = tmp - f; // compute error auto error = std::max(f.norm(), std::abs(yield_function)); return error; }; auto projection_error = update_f(k , alpha); /* --------------------------- */ /* iteration loop */ /* --------------------------- */ UInt iterations{0}; while(tolerance < projection_error and iterations < max_iterations) { // compute hessian at previous step compute_hessian_f(); // compute inverse matrix Xi Matrix xi(size, size); xi = Ce + plastic_multiplier_guess * hessian_f; // compute inverse matrix Xi Matrix xi_inv(size, size); xi_inv.inverse(xi); Vector tmp(size); tmp.mul(xi_inv, gradient_f); auto denominator = gradient_f.dot(tmp); // compute plastic multplier guess Vector tmp1(size); tmp1.mul(xi_inv, delta_inelastic_strain); plastic_multiplier_guess = gradient_f.dot(tmp1); plastic_multiplier_guess += yield_function; plastic_multiplier_guess /= denominator; // compute delta sigma Matrix tmp2(size, size); tmp2.outerProduct(tmp, tmp); tmp2 /= denominator; tmp2 = xi_inv - tmp2; delta_sigma.mul(tmp2, delta_inelastic_strain); delta_sigma -= tmp*yield_function/denominator; // compute sigma_guess Matrix delta_sigma_mat(dim, dim); voigt_h::voigtToMatrix(delta_sigma, delta_sigma_mat); sigma_guess += delta_sigma_mat; projection_error = update_f(k, alpha); iterations++; } } /* -------------------------------------------------------------------------- */ /// Infinitesimal deformations template inline void MaterialDruckerPrager::computeStressOnQuad( const Matrix & grad_u, const Matrix & previous_grad_u, Matrix & sigma, const Matrix & previous_sigma, Matrix & inelastic_strain, const Matrix & previous_inelastic_strain, const Real & sigma_th, const Real & previous_sigma_th) { Real delta_sigma_th = sigma_th - previous_sigma_th; Matrix grad_delta_u(grad_u); grad_delta_u -= previous_grad_u; // Compute trial stress, sigma_tr Matrix sigma_tr(dim, dim); MaterialElastic::computeStressOnQuad(grad_delta_u, sigma_tr, delta_sigma_th); sigma_tr += previous_sigma; - // Compute the yielding function + // Compute the yielding sress + /*Real yield_stress = computeYieldStress(sigma_tr); + + Matrix sigma_tr_dev(dim, dim, 0); + this->computeDeviatoricStress(sigma_tr, sigma_tr_dev); + + Real j2 = (1./ 2.) * sigma_tr_dev.doubleDot(sigma_tr_dev); + Real sigma_tr_dev_eff = std::sqrt(3. * j2); + + bool initial_yielding = ((sigma_tr_dev_eff + yield_stress) > 0);*/ + bool initial_yielding = (this->computeYieldFunction(sigma_tr) > 0); // use closet point projection to compute the plastic multiplier and // gradient and inealstic strain at the surface for the given trial stress state Matrix delta_inelastic_strain(dim, dim, 0.); if(initial_yielding) { UInt size = voigt_h::size; // plastic multiplier Real dp{0.}; // gradient of yield surface in voigt notation Vector gradient_f(size, 0.); // inelastic strain in voigt notation Vector delta_inelastic_strain_voigt(size, 0.); // compute using closet-point projection this->computeGradientAndPlasticMultplier(sigma_tr, dp, gradient_f, delta_inelastic_strain_voigt); for(auto i: arange(dim, size)) delta_inelastic_strain_voigt[i] /= 2.; voigt_h::voigtToMatrix(delta_inelastic_strain_voigt, delta_inelastic_strain); } // Compute the increment in inelastic strain MaterialPlastic::computeStressAndInelasticStrainOnQuad( grad_delta_u, sigma, previous_sigma, inelastic_strain, previous_inelastic_strain, delta_inelastic_strain); } /* -------------------------------------------------------------------------- */ /// Finite deformations template inline void MaterialDruckerPrager::computeStressOnQuad( __attribute__((unused)) const Matrix & grad_u, __attribute__((unused)) const Matrix & previous_grad_u, __attribute__((unused)) Matrix & sigma, __attribute__((unused)) const Matrix & previous_sigma, __attribute__((unused)) Matrix & inelastic_strain, __attribute__((unused)) const Matrix & previous_inelastic_strain, __attribute__((unused)) const Real & sigma_th, __attribute__((unused)) const Real & previous_sigma_th, __attribute__((unused)) const Matrix & F_tensor) { } /* -------------------------------------------------------------------------- */ template inline void MaterialDruckerPrager::computeTangentModuliOnQuad( __attribute__((unused)) Matrix & tangent, __attribute__((unused)) const Matrix & grad_u, __attribute__((unused)) const Matrix & previous_grad_u, __attribute__((unused)) const Matrix & sigma_tensor, __attribute__((unused)) const Matrix & previous_sigma_tensor) const { } }