diff --git a/src/model/solid_mechanics/materials/material_plastic/material_drucker_prager.cc b/src/model/solid_mechanics/materials/material_plastic/material_drucker_prager.cc index 3fe431942..b8033d546 100644 --- a/src/model/solid_mechanics/materials/material_plastic/material_drucker_prager.cc +++ b/src/model/solid_mechanics/materials/material_plastic/material_drucker_prager.cc @@ -1,199 +1,200 @@ /** * @file material_plastic.cc * * @author Mohit Pundir * * @date creation: Wed Sep 09 2020 * @date last modification: Wed Sep 09 2020 * * @brief Implementation of the akantu::MaterialDruckerPrager class * * * 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 { template MaterialDruckerPrager::MaterialDruckerPrager( SolidMechanicsModel & model, const ID & id) : MaterialPlastic(model, id) { AKANTU_DEBUG_IN(); this->initialize(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template MaterialDruckerPrager:: MaterialDruckerPrager(SolidMechanicsModel & model, UInt dim, const Mesh & mesh, FEEngine & fe_engine, const ID & id) : MaterialPlastic(model, dim, mesh, fe_engine, id) { AKANTU_DEBUG_IN(); this->initialize(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialDruckerPrager::initialize() { this->registerParam("phi", phi, Real(0.), _pat_parsable | _pat_modifiable, - "Internal friction angle"); + "Internal friction angle in degrees"); this->registerParam("fc", fc, Real(1.), _pat_parsable | _pat_modifiable, "Compressive strength"); this->updateInternalParameters(); } /* -------------------------------------------------------------------------- */ template void MaterialDruckerPrager::updateInternalParameters() { // compute alpha and k parameters for Drucker-Prager - this->alpha = (6.*sin(this->phi))/(3.-sin(this->phi)); - Real cohesion = this->fc * (1. - sin(this->phi))/(2.*cos(this->phi)); - this->k = (6. * cohesion * cos(this->phi))/(3. - sin(this->phi)); + Real phi_radian = this->phi * M_PI / 180.; + this->alpha = (6.*sin(phi_radian))/(3.-sin(phi_radian)); + Real cohesion = this->fc * (1. - sin(phi_radian))/(2.*cos(phi_radian)); + this->k = (6. * cohesion * cos(phi_radian))/(3. - sin(phi_radian)); } /* -------------------------------------------------------------------------- */ template void MaterialDruckerPrager::computeStress( ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); MaterialThermal::computeStress(el_type, ghost_type); // infinitesimal and finite deformation auto sigma_th_it = this->sigma_th(el_type, ghost_type).begin(); auto previous_sigma_th_it = this->sigma_th.previous(el_type, ghost_type).begin(); auto previous_gradu_it = this->gradu.previous(el_type, ghost_type) .begin(spatial_dimension, spatial_dimension); auto previous_stress_it = this->stress.previous(el_type, ghost_type) .begin(spatial_dimension, spatial_dimension); auto inelastic_strain_it = this->inelastic_strain(el_type, ghost_type) .begin(spatial_dimension, spatial_dimension); auto previous_inelastic_strain_it = this->inelastic_strain.previous(el_type, ghost_type) .begin(spatial_dimension, spatial_dimension); // // Finite Deformations // if (this->finite_deformation) { auto previous_piola_kirchhoff_2_it = this->piola_kirchhoff_2.previous(el_type, ghost_type) .begin(spatial_dimension, spatial_dimension); auto green_strain_it = this->green_strain(el_type, ghost_type) .begin(spatial_dimension, spatial_dimension); MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type); auto & inelastic_strain_tensor = *inelastic_strain_it; auto & previous_inelastic_strain_tensor = *previous_inelastic_strain_it; auto & previous_grad_u = *previous_gradu_it; auto & previous_sigma = *previous_piola_kirchhoff_2_it; auto & green_strain = *green_strain_it; this->template gradUToE(grad_u, green_strain); Matrix previous_green_strain(spatial_dimension, spatial_dimension); this->template gradUToE(previous_grad_u, previous_green_strain); Matrix F_tensor(spatial_dimension, spatial_dimension); this->template gradUToF(grad_u, F_tensor); computeStressOnQuad(green_strain, previous_green_strain, sigma, previous_sigma, inelastic_strain_tensor, previous_inelastic_strain_tensor, *sigma_th_it, *previous_sigma_th_it, F_tensor); ++sigma_th_it; ++inelastic_strain_it; ++previous_sigma_th_it; //++previous_stress_it; ++previous_gradu_it; ++green_strain_it; ++previous_inelastic_strain_it; ++previous_piola_kirchhoff_2_it; MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; } // Infinitesimal deformations else { MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type); auto & inelastic_strain_tensor = *inelastic_strain_it; auto & previous_inelastic_strain_tensor = *previous_inelastic_strain_it; auto & previous_grad_u = *previous_gradu_it; auto & previous_sigma = *previous_stress_it; computeStressOnQuad( grad_u, previous_grad_u, sigma, previous_sigma, inelastic_strain_tensor, previous_inelastic_strain_tensor, *sigma_th_it, *previous_sigma_th_it); ++sigma_th_it; ++inelastic_strain_it; ++previous_sigma_th_it; ++previous_stress_it; ++previous_gradu_it; ++previous_inelastic_strain_it; MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; } AKANTU_DEBUG_OUT(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialDruckerPrager::computeTangentModuli( __attribute__((unused)) const ElementType & el_type, __attribute__((unused)) Array & tangent_matrix, __attribute__((unused)) GhostType ghost_type) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ INSTANTIATE_MATERIAL(plastic_drucker_prager, MaterialDruckerPrager); } // namespace akantu 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 307c39e24..952d9ab80 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,331 +1,334 @@ /** * @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 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; 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; // lambda function to compute gradient of yield surface in voigt notation auto compute_gradient_f = [&sigma_guess, &scaling_matrix, &kronecker_delta, &gradient_f, this](){ const UInt dimension = sigma_guess.cols(); Matrix sigma_dev(dimension, dimension, 0); 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 += this->alpha * kronecker_delta; }; // lambda function to compute hessian matrix of yield surface auto compute_hessian_f = [&sigma_guess, &hessian_f, &scaling_matrix, &kronecker_delta, this](){ const UInt dimension = sigma_guess.cols(); Matrix sigma_dev(dimension, dimension, 0); this->computeDeviatoricStress(sigma_guess, sigma_dev); 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, this](){ // compute gradient compute_gradient_f(); // compute yield function yield_function = this->computeYieldFunction(sigma_guess); // 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(); /* --------------------------- */ /* 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(); 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 bool initial_yielding = - (this->computeYieldFunction(sigma_tr) < 0); + (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 { } }