diff --git a/FEM/CMakeLists.txt b/FEM/CMakeLists.txt index 81f57a6bd..c724e4454 100644 --- a/FEM/CMakeLists.txt +++ b/FEM/CMakeLists.txt @@ -1,5 +1,6 @@ set( HEADERS BoundaryCondition.h + burgers.h PhysicalConstant.h Constrained.h conversion_rate.h @@ -16,9 +17,11 @@ set( HEADERS files0.h GeoInfo.h InitialCondition.h + invariants.h LinearFunctionData.h mathlib.h matrix_class.h + minkley.h Output.h pcs_dm.h problem.h @@ -54,6 +57,7 @@ set( HEADERS set( SOURCES BoundaryCondition.cpp + burgers.cpp CAP_IO.cpp conversion_rate.cpp DistributionInfo.cpp @@ -72,9 +76,11 @@ set( SOURCES files0.cpp GeoInfo.cpp InitialCondition.cpp + invariants.cpp LinearFunctionData.cpp mathlib.cpp matrix_class.cpp + minkley.cpp Output.cpp pcs_dm.cpp problem.cpp diff --git a/FEM/burgers.cpp b/FEM/burgers.cpp new file mode 100644 index 000000000..8911159c3 --- /dev/null +++ b/FEM/burgers.cpp @@ -0,0 +1,156 @@ +/** + * \copyright + * Copyright (c) 2016, OpenGeoSys Community (http://www.opengeosys.org) + * Distributed under a Modified BSD License. + * See accompanying file LICENSE.txt or + * http://www.opengeosys.org/project/license + * + */ + +#include "burgers.h" +#include "PhysicalConstant.h" +namespace Burgers +{ +SolidBurgers::SolidBurgers(const Math_Group::Matrix& data) +{ + GK0 = data(0); // Kelvin shear modulus + mK = data(1); // dependency parameter for " + etaK0 = data(2); // Kelvin viscosity + mvK = data(3); // dependency parameter for " + GM0 = data(4); // Maxwell shear modulus + KM0 = data(5); // Maxwell bulk modulus + etaM0 = data(6); // Maxwell viscosity + mvM = data(7); // dependency parameter for " + m_GM = data(8); // slope of elesticity temperature dependence + m_KM = data(9); // slope of elesticity temperature dependence + T_ref = data(10); // reference temperature dependency parameter for " + B = data(11); // constant factor for Arrhenius term + Q = data(12); // activation energy in Arrhenius term + + GM = GM0; + KM = KM0; + GK = GK0; + etaK = etaK0; + etaM = etaM0; + + if (!T_Process) + { + m_GM = 0.; + m_KM = 0.; + B = 1.; + Q = 0.; // for cutting off Arrhenius term + T_ref = 273.15; + } +} + +/************************************************************************** + FEMLib-Method: Burgers::UpdateBurgersProperties() + Task: Updates BURGERS material parameters in LUBBY2 fashion + Programing: + 07/2014 TN Implementation +**************************************************************************/ +void SolidBurgers::UpdateBurgersProperties(double s_eff, const double Temperature) +{ + const double dT(Temperature - T_ref); + GM = GM0 + m_GM * dT; + KM = KM0 + m_KM * dT; + + s_eff *= GM; + + GK = GK0 * std::exp(mK * s_eff); + etaK = etaK0 * std::exp(mvK * s_eff); + etaM = etaM0 * std::exp(mvM * s_eff) * B + * std::exp(Q * (-dT) / (PhysicalConstant::IdealGasConstant * Temperature * T_ref)); + // if (etaM / etaM0 < 1.e-2) + // std::cout << "WARNING: Maxwell viscosity sank to 100th of original value." << std::endl; +} + +/************************************************************************** + FEMLib-Method: Burgers::CalResidualBurgers() + Task: Calculates the 12x1 residual vector. Implementation fully implicit only. + Programing: + 06/2014 TN Implementation +**************************************************************************/ +void SolidBurgers::CalResidualBurgers(const double dt, const KVec& strain_curr, const KVec& stress_curr, + KVec& strain_Kel_curr, const KVec& strain_Kel_t, KVec& strain_Max_curr, + const KVec& strain_Max_t, Eigen::Matrix& res) +{ + // calculate stress residual + res.block<6, 1>(0, 0) = stress_curr - 2. * (strain_curr - strain_Kel_curr - strain_Max_curr); + + // calculate Kelvin strain residual + res.block<6, 1>(6, 0) = 1. / dt * (strain_Kel_curr - strain_Kel_t) + - 1. / (2. * etaK) * (GM * stress_curr - 2. * GK * strain_Kel_curr); + + // calculate Maxwell strain residual + res.block<6, 1>(12, 0) = 1. / dt * (strain_Max_curr - strain_Max_t) - GM / (2. * etaM) * stress_curr; +} + +/************************************************************************** + FEMLib-Method: Burgers::CalJacobianBurgers() + Task: Calculates the 12x12 Jacobian. Implementation fully implicit only. + Programing: + 06/2014 TN Implementation +**************************************************************************/ +void SolidBurgers::CalJacobianBurgers(const double dt, Eigen::Matrix& Jac, const double s_eff, + const KVec& sig_i, const KVec& eps_K_i) +{ + // Check Dimension of Jacobian + // assert(Jac.cols() == 18 && Jac.rows() && 18); + Jac.setZero(18, 18); + + // build G_11 + Jac.block<6, 6>(0, 0) = SolidMath::ident; + + // build G_12 + Jac.block<6, 6>(0, 6) = 2. * SolidMath::ident; + + // build G_13 + Jac.block<6, 6>(0, 12) = 2. * SolidMath::ident; + + // build G_21 + Jac.block<6, 6>(6, 0) = -GM / (2. * etaK) * SolidMath::ident; + + // build G_22 + Jac.block<6, 6>(6, 6) = (1. / dt + GK / etaK) * SolidMath::ident; + + // nothing to do for G_23 + + // build G_31 + Jac.block<6, 6>(12, 0) = -GM / (2. * etaM) * SolidMath::ident; + + // nothing to do for G_32 + + // build G_33 + Jac.block<6, 6>(12, 12) = 1. / dt * SolidMath::ident; + + if (s_eff > 0.) + { + const KVec dG_K = mK * 3. * GK * GM / (2. * s_eff) * sig_i; + const KVec dmu_vK = mvK * 3. * GM * etaK / (2. * s_eff) * sig_i; + const KVec dmu_vM = mvM * 3. * GM * etaM / (2. * s_eff) * sig_i; + const KVec eps_K_aid = 1. / (etaK * etaK) * (GM * sig_i - 2. * GK * eps_K_i); + + // build G_21 + Jac.block<6, 6>(6, 0) += 0.5 * eps_K_aid * dmu_vK.transpose() + 1. / etaK * eps_K_i * dG_K.transpose(); + + // build G_31 + Jac.block<6, 6>(12, 0) += GM / (2. * etaM * etaM) * sig_i * dmu_vM.transpose(); + } +} + +/************************************************************************** + FEMLib-Method: Burgers::CaldGdEBurgers() + Task: Calculates the 12x6 derivative of the residuals with respect to total strain. Implementation fully implicit +only. + Programing: + 06/2014 TN Implementation +**************************************************************************/ +void SolidBurgers::CaldGdEBurgers(Eigen::Matrix& dGdE) +{ + // Check Dimension of dGdE + // assert(dGdE.cols() == 6 && dGdE.rows() == 18); + dGdE.setZero(18, 6); + dGdE.block<6, 6>(0, 0) = -2. * SolidMath::P_dev; +} +} diff --git a/FEM/burgers.h b/FEM/burgers.h new file mode 100644 index 000000000..e2835e118 --- /dev/null +++ b/FEM/burgers.h @@ -0,0 +1,44 @@ +/** + * \copyright + * Copyright (c) 2016, OpenGeoSys Community (http://www.opengeosys.org) + * Distributed under a Modified BSD License. + * See accompanying file LICENSE.txt or + * http://www.opengeosys.org/project/license + * + */ + +#ifndef BURGERS_H +#define BURGERS_H + +#include +#include "Eigen/Dense" +// FEM-Makros +#include "makros.h" +#include "tools.h" +#include "invariants.h" + +namespace Burgers +{ +typedef Eigen::Matrix KVec; +typedef Eigen::Matrix KMat; + +class SolidBurgers +{ +public: + explicit SolidBurgers(const Math_Group::Matrix& data); + ~SolidBurgers() {} + // basic material parameters + double GK0, GM0, KM0, etaK0, etaM0, mK, mvK, mvM, l0, T_ref, m_GM, m_KM, B, Q; + // solution dependent values + double GM, KM, GK, etaK, etaM; + + void UpdateBurgersProperties(double s_eff, const double Delta_T); + void CalResidualBurgers(const double dt, const KVec& strain_curr, const KVec& stress_curr, KVec& strain_Kel_curr, + const KVec& strain_Kel_t, KVec& strain_Max_curr, const KVec& strain_Max_t, + Eigen::Matrix& res); + void CalJacobianBurgers(const double dt, Eigen::Matrix& Jac, const double s_eff, const KVec& sig_i, + const KVec& eps_K_i); + void CaldGdEBurgers(Eigen::Matrix& dGdE); +}; +} +#endif // BURGERS_H diff --git a/FEM/equation_class.h b/FEM/equation_class.h index 264906c1e..151c9a634 100644 --- a/FEM/equation_class.h +++ b/FEM/equation_class.h @@ -120,7 +120,7 @@ class Linear_EQS double X(const long i) const { return x[i]; } double RHS(const long i) const { return b[i]; } double NormX(); - double NormRHS() { return bNorm; } + double NormRHS() { return Norm(b); } #if defined(USE_MPI) int DOF() { return A->Dof(); } long Size() { return A->Size(); } diff --git a/FEM/fem_ele_vec.cpp b/FEM/fem_ele_vec.cpp index 03ab179cf..b2116393a 100644 --- a/FEM/fem_ele_vec.cpp +++ b/FEM/fem_ele_vec.cpp @@ -960,7 +960,17 @@ void CFiniteElementVec::ComputeMatrix_RHS(const double fkt, const Matrix* p_D) // 2D, in y-direction // 3D, in z-direction i = (ele_dim - 1) * nnodesHQ; - const double coeff = LoadFactor * rho * smat->grav_const * fkt; + double timeFactor = 1.; + if (smat->gravity_ramp) + { + const int curveIndex = smat->grav_curve_id; + const int curveMethod = 0; + + int valid; + timeFactor = GetCurveValue(curveIndex, curveMethod, aktuelle_zeit, &valid); + } + + const double coeff = LoadFactor * rho * smat->grav_const * fkt * timeFactor; for (k = 0; k < nnodesHQ; k++) { (*RHS)[i + k] += coeff * shapefctHQ[k]; @@ -2420,6 +2430,127 @@ void CFiniteElementVec::GlobalAssembly_RHS() for (i = 0; i < ns; i++) // JT: This was commented. It shouldn't be. dstrain[i] += strain_ne[i]; } + + if (smat->Creep_mode == 1001) // BURGERS. + { + // get total strains at time t and current iteration i + std::vector strain_t(6), strain_curr(6); + + // get stresses as well as internal variables at time t and current iteration i (equal before local + // Newton) + std::vector stress_curr(6), eps_K_curr(6), eps_M_curr(6); + for (int compnt(0); compnt < ns; compnt++) + { + strain_t[compnt] = (*eleV_DM->Strain_t_ip)(compnt, gp); + strain_curr[compnt] = strain_t[compnt] + dstrain[compnt]; + stress_curr[compnt] = (*eleV_DM->Stress)(compnt, gp); + eps_K_curr[compnt] = (*eleV_DM->Strain_Kel)(compnt, gp); + eps_M_curr[compnt] = (*eleV_DM->Strain_Max)(compnt, gp); + } + // Fill remaining 3D components with 0 if only 2D + if (ns == 4) + for (int compnt(4); compnt < 6; compnt++) + { + strain_t[compnt] = 0.0; + strain_curr[compnt] = 0.0; + stress_curr[compnt] = 0.0; + eps_K_curr[compnt] = 0.0; + eps_M_curr[compnt] = 0.0; + } + + // 6x6 tangent + Matrix ConsD(6, 6); + // Pass as 6D vectors, i.e. set stress and strain [4] and [5] to zero for 2D and AXI as well as + // strain[3] to zero for 2D (plane strain) + double local_res; + smat->LocalNewtonBurgers(dt, strain_curr, stress_curr, eps_K_curr, eps_M_curr, ConsD, t1, + local_res); + + // Then update (and reduce for 2D) stress increment vector and reduce (for 2D) ConsistDep, update + // internal variables + for (int compnt(0); compnt < ns; compnt++) + { + dstress[compnt] = stress_curr[compnt]; + if (update > 0) + { + (*eleV_DM->Strain_Kel)(compnt, gp) = eps_K_curr[compnt]; + (*eleV_DM->Strain_Max)(compnt, gp) = eps_M_curr[compnt]; + (*eleV_DM->Strain_t_ip)(compnt, gp) = strain_curr[compnt]; + (*eleV_DM->ev_loc_nr_res)(gp) = local_res; + } + for (int compnt2(0); compnt2 < ns; compnt2++) + (*De)(compnt, compnt2) = ConsD(compnt, compnt2); + (*eleV_DM->Strain)(compnt, gp) = strain_curr[compnt]; + } + } + + if (smat->Creep_mode == 1002) // MINKLEY + { + // get total strains at time t and current iteration i + std::vector strain_t(6), strain_curr(6); + + std::vector stress_curr(6), eps_K_curr(6), eps_M_curr(6), eps_pl_curr(6); + for (int compnt(0); compnt < ns; compnt++) + { + stress_curr[compnt] = (*eleV_DM->Stress)(compnt, gp); + eps_K_curr[compnt] = (*eleV_DM->Strain_Kel)(compnt, gp); + eps_M_curr[compnt] = (*eleV_DM->Strain_Max)(compnt, gp); + eps_pl_curr[compnt] = (*eleV_DM->Strain_pl)(compnt, gp); + strain_t[compnt] = (*eleV_DM->Strain_t_ip)(compnt, gp); + strain_curr[compnt] = strain_t[compnt] + dstrain[compnt]; + } + // Fill remaining 3D components with 0 if only 2D + if (ns == 4) + for (int compnt(4); compnt < 6; compnt++) + { + strain_t[compnt] = 0.0; + strain_curr[compnt] = 0.0; + stress_curr[compnt] = 0.0; + eps_K_curr[compnt] = 0.0; + eps_M_curr[compnt] = 0.0; + eps_pl_curr[compnt] = 0.0; + } + + double e_pl_v = (*eleV_DM->e_pl)(gp); + double e_pl_eff = (*eleV_DM->pStrain)(gp); + double lam = (*eleV_DM->lambda_pl)(gp);//NOTE: May set starting value to zero in case of trouble with load reversals + + // 6x6 tangent + Matrix ConsD(6, 6); + + // Pass as 6D vectors, i.e. set stress and strain [4] and [5] to zero for 2D and AXI as well as + // strain[3] to zero for 2D (plane strain) + double local_res(0.); + smat->LocalNewtonMinkley(dt, strain_curr, stress_curr, eps_K_curr, eps_M_curr, eps_pl_curr, e_pl_v, + e_pl_eff, lam, ConsD, t1, local_res); + + // Then update (and reduce for 2D) stress increment vector and reduce (for 2D) ConsistDep, update + // internal variables + for (int compnt(0); compnt < ns; compnt++) + { + dstress[compnt] = stress_curr[compnt]; + if (update > 0) + { + (*eleV_DM->Strain_Kel)(compnt, gp) = eps_K_curr[compnt]; + (*eleV_DM->Strain_Max)(compnt, gp) = eps_M_curr[compnt]; + (*eleV_DM->Strain_pl)(compnt, gp) = eps_pl_curr[compnt]; + (*eleV_DM->Strain_t_ip)(compnt, gp) = strain_curr[compnt]; + } + for (int compnt2(0); compnt2 < ns; compnt2++) + (*De)(compnt, compnt2) = ConsD(compnt, compnt2); + (*eleV_DM->Strain)(compnt, gp) = strain_curr[compnt]; + } + + if (update > 0) + { + (*eleV_DM->e_pl)(gp) = e_pl_v; + (*eleV_DM->lambda_pl)(gp) = lam; + (*eleV_DM->pStrain)(gp) = e_pl_eff; + (*eleV_DM->ev_loc_nr_res)(gp) = local_res; + } + } + + // Fluid coupling; S_Water = 1.0; if (Flow_Type > 0 && Flow_Type != 10) @@ -3659,7 +3790,7 @@ void CFiniteElementVec::GlobalAssembly_RHS() ---------------------- -----------------------------------------------------------------*/ ElementValue_DM::ElementValue_DM(CElem * ele, const int NGP, bool HM_Staggered) - : NodesOnPath(NULL), orientation(NULL) + : NodesOnPath(NULL), orientation(NULL), Strain(NULL), Strain_t_ip(NULL) { int Plastic = 1; const int LengthMat = 7; // Number of material parameter of SYS model. @@ -3680,6 +3811,12 @@ void CFiniteElementVec::GlobalAssembly_RHS() ele_dim = ele->GetDimension(); sdp = msp_vector[ele->GetPatchIndex()]; Plastic = sdp->Plastictity(); + Strain_Kel = NULL; + Strain_Max = NULL; + Strain_pl = NULL; + e_pl = NULL; + lambda_pl = NULL; + ev_loc_nr_res = NULL; if (ele_dim == 2) LengthBS = 4; @@ -3733,6 +3870,46 @@ void CFiniteElementVec::GlobalAssembly_RHS() xi = new Matrix(LengthBS); *xi = 0.0; } + if (sdp->CreepModel() == 1001) // Burgers + { + Strain_Kel = new Matrix(6, NGPoints); // Only 3D size for now. Will work with all. Separation into special + // cases (LengthBS) may follow later + *Strain_Kel = 0.0; + Strain_Max = new Matrix(6, NGPoints); // Only 3D size for now. Will work with all. Separation into special + // cases (LengthBS) may follow later + *Strain_Max = 0.0; + Strain_t_ip = new Matrix(6, NGPoints); + *Strain_t_ip = 0.0; + ev_loc_nr_res = new Matrix(NGPoints); + *ev_loc_nr_res = 0.; + Strain = new Matrix(6, NGPoints); + *Strain = 0.0; + } + + if (sdp->CreepModel() == 1002) // Minkley + { + Strain_Kel = new Matrix(6, NGPoints); // Only 3D size for now. Will work with all. Separation into special + // cases (LengthBS) may follow later + *Strain_Kel = 0.0; + Strain_Max = new Matrix(6, NGPoints); // Only 3D size for now. Will work with all. Separation into special + // cases (LengthBS) may follow later + *Strain_Max = 0.0; + Strain_pl = new Matrix(6, NGPoints); // Only 3D size for now. Will work with all. Separation into special + // cases (LengthBS) may follow later + *Strain_pl = 0.0; + Strain_t_ip = new Matrix(6, NGPoints); + *Strain_t_ip = 0.0; + e_pl = new Matrix(NGPoints); + *e_pl = 0.; + pStrain = new Matrix(NGPoints); + *pStrain = 0.; + lambda_pl = new Matrix(NGPoints); + *lambda_pl = 0.; + ev_loc_nr_res = new Matrix(NGPoints); + *ev_loc_nr_res = 0.; + Strain = new Matrix(6, NGPoints); + *Strain = 0.0; + } disp_j = 0.0; tract_j = 0.0; Localized = false; @@ -3760,6 +3937,22 @@ void CFiniteElementVec::GlobalAssembly_RHS() y_surface->Write_BIN(os); if (xi) xi->Write_BIN(os); + if (Strain) // NB + Strain->Write_BIN(os); + if (Strain_Kel) // TN Burgers + Strain_Kel->Write_BIN(os); + if (Strain_Max) // TN Burgers + Strain_Max->Write_BIN(os); + if (Strain_pl) // TN + Strain_pl->Write_BIN(os); + if (Strain_t_ip) + Strain_t_ip->Write_BIN(os); + if (e_pl) + e_pl->Write_BIN(os); + if (lambda_pl) + lambda_pl->Write_BIN(os); + if (ev_loc_nr_res) + ev_loc_nr_res->Write_BIN(os); if (MatP) MatP->Write_BIN(os); if (prep0) @@ -3785,6 +3978,22 @@ void CFiniteElementVec::GlobalAssembly_RHS() y_surface->Read_BIN(is); if (xi) xi->Read_BIN(is); + if (Strain) // NB + Strain->Read_BIN(is); + if (Strain_Kel) // TN Burgers + Strain_Kel->Read_BIN(is); + if (Strain_Max) // TN Burgers + Strain_Max->Read_BIN(is); + if (Strain_pl) // TN + Strain_pl->Read_BIN(is); + if (Strain_t_ip) + Strain_t_ip->Read_BIN(is); + if (e_pl) + e_pl->Read_BIN(is); + if (lambda_pl) + lambda_pl->Read_BIN(is); + if (ev_loc_nr_res) + ev_loc_nr_res->Read_BIN(is); if (MatP) MatP->Read_BIN(is); if (prep0) @@ -3863,6 +4072,22 @@ void CFiniteElementVec::GlobalAssembly_RHS() delete scalar_aniso_comp; // WX:09.2011 if (scalar_aniso_tens) delete scalar_aniso_tens; + if (Strain) + delete Strain; + if (Strain_Kel) // TN Strain in Kelvin element + delete Strain_Kel; + if (Strain_Max) // TN Strain in Maxwell element + delete Strain_Max; + if (Strain_pl) // TN - Minkley plastic deviatoric strain + delete Strain_pl; + if (Strain_t_ip) + delete Strain_t_ip; + if (e_pl) + delete e_pl; + if (lambda_pl) + delete lambda_pl; + if (ev_loc_nr_res) + delete ev_loc_nr_res; NodesOnPath = NULL; orientation = NULL; @@ -3878,6 +4103,14 @@ void CFiniteElementVec::GlobalAssembly_RHS() MatP = NULL; scalar_aniso_comp = NULL; scalar_aniso_tens = NULL; + Strain = NULL; + Strain_Kel = NULL; + Strain_Max = NULL; + Strain_pl = NULL; + Strain_t_ip = NULL; + e_pl = NULL; + lambda_pl = NULL; + ev_loc_nr_res = NULL; } /*************************************************************************** diff --git a/FEM/fem_ele_vec.h b/FEM/fem_ele_vec.h index 112617dc6..72f562998 100644 --- a/FEM/fem_ele_vec.h +++ b/FEM/fem_ele_vec.h @@ -79,7 +79,14 @@ class ElementValue_DM // Variables of single yield surface model Matrix* xi; // Rotational hardening variables Matrix* MatP; // Material parameters - + Matrix* Strain_Kel; // TN - Burgers model dev. Kelvin strain + Matrix* Strain_Max; // TN - Burgers model dev. Maxwell strain + Matrix* Strain_pl; // TN - Minkley - deviatoric plastic strain + Matrix* Strain_t_ip; // TN - introduced to save previous strain in all integration points + Matrix* e_pl; + Matrix* ev_loc_nr_res; + Matrix* lambda_pl; // TN - Minkley - volumetric plastic strain, plastic arc length, plastic multiplier + Matrix* Strain; // NB - Strain tensor for reload-feature // Discontinuity double disp_j; double tract_j; diff --git a/FEM/invariants.cpp b/FEM/invariants.cpp new file mode 100644 index 000000000..c4ab55521 --- /dev/null +++ b/FEM/invariants.cpp @@ -0,0 +1,227 @@ +/** + * \copyright + * Copyright (c) 2016, OpenGeoSys Community (http://www.opengeosys.org) + * Distributed under a Modified BSD License. + * See accompanying file LICENSE.txt or + * http://www.opengeosys.org/project/license + * + */ + +#include "invariants.h" + +namespace SolidMath +{ +KMat Initialise_ident() +{ + KMat I; + I.setIdentity(); + return I; +} + +KMat Initialise_P_dev() +{ + KMat P; + P.setIdentity(); + const double third(1. / 3.); + for (size_t i = 0; i < 3; i++) + for (size_t j = 0; j < 3; j++) + { + P(i, j) -= third; + } + return P; +} + +KMat Initialise_P_sph() +{ + KMat P; + P.setZero(6, 6); + const double third(1. / 3.); + for (size_t i = 0; i < 3; i++) + for (size_t j = 0; j < 3; j++) + { + P(i, j) = third; + } + return P; +} + +KVec Initialise_ivec() +{ + KVec I; + for (size_t j = 0; j < 3; j++) + { + I(j) = 1.; + I(j + 3) = 0.; + } + return I; +} + +/************************************************************************** + Voigt_to_Kelvin_Strain + Task: Maps a strain vector in Voigt notation into on in Kelvin notation + This is an auxilliary routine that will not be needed when the entire FE + code is set up in Kelvin notation + Programing: + 06/2015 TN Implementation +**************************************************************************/ +KVec Voigt_to_Kelvin_Strain(const std::vector& voigt_strain) +{ + KVec kelvin_strain; + for (size_t i = 0; i < 3; i++) + { + // Normal components + kelvin_strain(i) = voigt_strain[i]; + // Shear components + kelvin_strain(i + 3) = voigt_strain[i + 3] / sqrt(2.); + } + return kelvin_strain; +} + +/************************************************************************** + Voigt_to_Kelvin_Stress + Task: Maps a stress vector in Voigt notation into on in Kelvin notation + This is an auxilliary routine that will not be needed when the entire FE + code is set up in Kelvin notation + Programing: + 06/2015 TN Implementation +**************************************************************************/ +KVec Voigt_to_Kelvin_Stress(const std::vector& voigt_stress) +{ + KVec kelvin_stress; + for (size_t i = 0; i < 3; i++) + { + // Normal components + kelvin_stress(i) = voigt_stress[i]; + // Shear components + kelvin_stress(i + 3) = voigt_stress[i + 3] * sqrt(2.); + } + return kelvin_stress; +} + +/************************************************************************** + Kelvin_to_Voigt_Strain() + Task: Maps a strain vector in Kelvin notation into on in Voigt notation + This is an auxilliary routine that will not be needed when the entire FE + code is set up in Kelvin notation + Programing: + 06/2015 TN Implementation +**************************************************************************/ +void Kelvin_to_Voigt_Strain(const KVec& kelvin_strain, std::vector& voigt_strain) +{ + for (size_t i = 0; i < 3; i++) + { + // Normal components + voigt_strain[i] = kelvin_strain(i); + // Shear components + voigt_strain[i + 3] = kelvin_strain(i + 3) * sqrt(2.); + } + return; +} + +/************************************************************************** + Kelvin_to_Voigt_Stress() + Task: Maps a stress vector in Kelvin notation into on in Voigt notation + This is an auxilliary routine that will not be needed when the entire FE + code is set up in Kelvin notation + Programing: + 06/2015 TN Implementation +**************************************************************************/ +void Kelvin_to_Voigt_Stress(const KVec& kelvin_stress, std::vector& voigt_stress) +{ + for (size_t i = 0; i < 3; i++) + { + // Normal components + voigt_stress[i] = kelvin_stress(i); + // Shear components + voigt_stress[i + 3] = kelvin_stress(i + 3) / sqrt(2.); + } + return; +} + +// Maps a 6D Kelvin vector back into 3D Tensor coordinates +Eigen::Matrix KelvinVectorToTensor(const KVec& vec) +{ + Eigen::Matrix tens; + tens(0, 0) = vec(0); + tens(1, 1) = vec(1); + tens(2, 2) = vec(2); + tens(0, 1) = tens(1, 0) = vec(3) / std::sqrt(2.); + tens(1, 2) = tens(2, 1) = vec(4) / std::sqrt(2.); + tens(0, 2) = tens(2, 0) = vec(5) / std::sqrt(2.); + return tens; +} + +// Maps a 2nd order 3D Tensor into Kelvin representation +KVec TensorToKelvinVector(const Eigen::Matrix& tens) +{ + KVec vec; + vec(0) = tens(0, 0); + vec(1) = tens(1, 1); + vec(2) = tens(2, 2); + vec(3) = tens(0, 1) * std::sqrt(2.); + vec(4) = tens(1, 2) * std::sqrt(2.); + vec(5) = tens(0, 2) * std::sqrt(2.); + return vec; +} + +// Task: calculates second deviatoric invariant. Note that this routine requires a +// Kelvin mapped DEVIATORIC vector. +// A deviatoric mapping was not done here in order to avoid unnecessary calculations +double CalJ2(const KVec& dev_vec) +{ + double s(0.); + s = dev_vec.transpose() * dev_vec; // Kelvin mapping, deviator + return 0.5 * s; +} + +// Task: calculates effective stress. Note that this routine requires a +// Kelvin mapped DEVIATORIC stress vector +double CalEffectiveStress(const KVec& dev_stress) +{ + double s(0.); + s = 3. * CalJ2(dev_stress); // Kelvin mapped deviatoric stress has to be used + return sqrt(s); +} + +// Task: calculates third deviatoric invariant. Note that this routine requires a +// Kelvin mapped DEVIATORIC vector +// A deviatoric mapping was not done here in order to avoid unnecessary calculations +double CalJ3(const KVec& dev_vec) +{ + Eigen::Matrix tens; + tens = KelvinVectorToTensor(dev_vec); // TN: Not strictly necessary. Can be written explicitly for vector + // coordinates + return tens.determinant(); +} + +// Takes 2nd Order tensor in Kelvin representation and returns its +// inverse in Kelvin representation +KVec InvertVector(const KVec& vec) +{ + Eigen::Matrix tens = KelvinVectorToTensor(vec); + return TensorToKelvinVector(tens.inverse()); +} + +// calculates Lode angle. Note that this routine requires a +// Kelvin mapped DEVIATORIC vector +// A deviatoric mapping was not done here in order to avoid unnecessary calculations +double CalLodeAngle(const KVec& dev_vec) +{ + const double J2(CalJ2(dev_vec)); + double theta, thetaR; + + if (std::abs(J2) > DBL_EPSILON) // catch zero-stress state + theta = -3. * std::sqrt(3.) * CalJ3(dev_vec) / (2. * std::pow(J2, 1.5)); + else + theta = 0.; + + thetaR = std::min(std::max(theta, -1. + DBL_EPSILON), 1. - DBL_EPSILON); // strict limits for trigonometric + // functions + return 1. / 3. * std::asin(thetaR); +} + +// calculates first invariant (trace) of a vector-mapped 2nd order tensor. +double CalI1(const KVec& vec) +{ + return vec(0) + vec(1) + vec(2); +} +} diff --git a/FEM/invariants.h b/FEM/invariants.h new file mode 100644 index 000000000..1d30d94d2 --- /dev/null +++ b/FEM/invariants.h @@ -0,0 +1,72 @@ +/** + * \copyright + * Copyright (c) 2016, OpenGeoSys Community (http://www.opengeosys.org) + * Distributed under a Modified BSD License. + * See accompanying file LICENSE.txt or + * http://www.opengeosys.org/project/license + * + */ + +#ifndef INVARIANTS_H +#define INVARIANTS_H + +#include +#include "Eigen/Dense" +// FEM-Makros +#include "makros.h" +#include "tools.h" + +namespace SolidMath +{ +typedef Eigen::Matrix KVec; +typedef Eigen::Matrix KMat; + +KMat Initialise_ident(); +KMat Initialise_P_dev(); +KMat Initialise_P_sph(); +KVec Initialise_ivec(); + +static const KMat ident(Initialise_ident()); // identity matrix +static const KMat P_dev(Initialise_P_dev()); // deviatoric projection matrix +static const KMat P_sph(Initialise_P_sph()); // spherical projection matrix +static const KVec ivec(Initialise_ivec()); // Kelvin mapping of 2nd order identity + +// Kelvin/Voigt mapping routines for 6D vectors +KVec Voigt_to_Kelvin_Stress(const std::vector& voigt_stress); +KVec Voigt_to_Kelvin_Strain(const std::vector& voigt_strain); +void Kelvin_to_Voigt_Stress(const KVec& kelvin_stress, std::vector& voigt_stress); +void Kelvin_to_Voigt_Strain(const KVec& kelvin_strain, std::vector& voigt_strain); + +// Maps a 6D Kelvin vector back into 3D Tensor coordinates +Eigen::Matrix KelvinVectorToTensor(const KVec& vec); + +// Maps a 2nd order 3D Tensor into Kelvin representation +KVec TensorToKelvinVector(const Eigen::Matrix& tens); + +// Task: calculates second deviatoric invariant. Note that this routine requires a +// Kelvin mapped DEVIATORIC vector. +// A deviatoric mapping was not done here in order to avoid unnecessary calculations +double CalJ2(const KVec& dev_vec); + +// Task: calculates effective stress. Note that this routine requires a +// Kelvin mapped DEVIATORIC stress vector +double CalEffectiveStress(const KVec& dev_stress); + +// Task: calculates third deviatoric invariant. Note that this routine requires a +// Kelvin mapped DEVIATORIC vector +// A deviatoric mapping was not done here in order to avoid unnecessary calculations +double CalJ3(const KVec& dev_vec); + +// Takes 2nd Order tensor in Kelvin representation and returns its +// inverse in Kelvin representation +KVec InvertVector(const KVec& vec); + +// calculates Lode angle. Note that this routine requires a +// Kelvin mapped DEVIATORIC vector +// A deviatoric mapping was not done here in order to avoid unnecessary calculations +double CalLodeAngle(const KVec& dev_vec); + +// calculates first invariant (trace) of a vector-mapped 2nd order tensor. +double CalI1(const KVec& vec); +} +#endif // INVARIANTS_H diff --git a/FEM/minkley.cpp b/FEM/minkley.cpp new file mode 100644 index 000000000..ee97dd29e --- /dev/null +++ b/FEM/minkley.cpp @@ -0,0 +1,658 @@ +/** + * \copyright + * Copyright (c) 2016, OpenGeoSys Community (http://www.opengeosys.org) + * Distributed under a Modified BSD License. + * See accompanying file LICENSE.txt or + * http://www.opengeosys.org/project/license + * + */ + +#include "minkley.h" + +namespace Minkley +{ +// Template to implement Signum Function +template +int sgn(T val) +{ + return (T(0) < val) - (val < T(0)); +} + +SolidMinkley::SolidMinkley(const Math_Group::Matrix& data) +{ + GK0 = data(0); // Kelvin shear modulus + etaK0 = data(1); // Kelvin viscosity + GM0 = data(2); // Maxwell shear modulus + KM0 = data(3); // Maxwell bulk modulus + etaM0 = data(4); // Maxwell viscosity + if (std::abs(data(5)) < 1.e-3) + { + mvM = std::log(1. + std::sqrt(2.)); // m -- effective stress factor in viscosity relationship + nvM = 0.; // n -- effective stress exponent in viscosity relationship + std::cout << "WARNING. Viscosity parameter m was chosen lower than 1e-3. " + "Setting model parameters to achieve constant viscosity."; + } + else + { + mvM = data(5); // m -- effective stress factor in viscosity relationship + nvM = data(6); // n -- effective stress exponent in viscosity relationship + } + coh0 = data(7); // initial cohesion + hard = data(8); // hardening/softening modulus + phi = data(9) * PI / 180.; // friction angle + psi = data(10) * PI / 180.; // dilatancy angle + thetaT = data(11) * PI / 180.; // transition angle + eta_reg = data(12); // viscosity for viscoplastic regularisation + m_GM = data(13); // slope of elesticity temperature dependence + m_KM = data(14); // slope of elesticity temperature dependence + T_ref = data(15); // reference temperature dependency parameter for " + Bt = data(16); // constant factor for Arrhenius term + Q = data(17); // activation energy in Arrhenius term + hard2 = data(18); //second order hardening term + hard4 = data(19); //second order hardening term + + etaM = etaM0; + coh = coh0; + GM = GM0; + KM = KM0; + + if (!T_Process) + { + m_GM = 0.; + m_KM = 0.; + Bt = 1.; + Q = 0.; // for cutting off Arrhenius term + T_ref = 273.15; + } +} + +/************************************************************************** + FEMLib-Method: Minkley::UpdateMinkleyProperties() + Task: Updates BURGERS material parameters in LUBBY2 fashion + Programing: + 06/2015 TN Implementation +**************************************************************************/ +void SolidMinkley::UpdateMinkleyProperties(double s_eff, const double eps_p_eff, double Temperature) +{ + const double dT(Temperature - T_ref); + GM = GM0 + m_GM * dT; + KM = KM0 + m_KM * dT; + + s_eff *= GM; + + if (s_eff > DBL_EPSILON) + etaM = etaM0 / std::sinh(mvM * std::pow(s_eff, nvM)); // viscosity function update + else + etaM = etaM0; + etaM *= Bt * std::exp(Q * (-dT) / (PhysicalConstant::IdealGasConstant * Temperature * T_ref)); + + coh = coh0 + * (1. + + eps_p_eff + * (hard + + eps_p_eff + * (hard2 + eps_p_eff * eps_p_eff * hard4))); // fourth order isotropic hardening/softening + // if (etaM / etaM0 < 1.e-2) + // std::cout << "WARNING: Maxwell viscosity sank to 100th of original value." << std::endl; +} + +/************************************************************************** + FEMLib-Method: Minkley::A() + Task: Expression A in Sloan's yield function + Programing: + 06/2015 TN Implementation +**************************************************************************/ +double SolidMinkley::A(const double theta, const double alpha) +{ + double A; + A = std::cos(thetaT) / 3. + * (3. + std::tan(thetaT) * std::tan(3. * thetaT) + + 1. / std::sqrt(3.) * sgn(theta) * (std::tan(3. * thetaT) - 3. * std::tan(thetaT)) * std::sin(alpha)); + return A; +} + +/************************************************************************** + FEMLib-Method: Minkley::B() + Task: Expression B in Sloan's yield function + Programing: + 06/2015 TN Implementation +**************************************************************************/ +double SolidMinkley::B(const double theta, const double alpha) +{ + double B; + B = 1. / (3. * std::cos(3. * thetaT)) + * (sgn(theta) * std::sin(thetaT) + std::sin(alpha) * std::cos(thetaT) / std::sqrt(3.)); + return B; +} + +/************************************************************************** + FEMLib-Method: Minkley::YieldMohrCoulomb() + Task: Yield function Mohr Coulomb with corner smoothing (Sloan et al.) + Programing: + 06/2015 TN Implementation +**************************************************************************/ +double SolidMinkley::YieldMohrCoulomb(const KVec& sig) +{ + double F(0.); + const KVec sigd(SolidMath::P_dev * sig); + const double theta(SolidMath::CalLodeAngle(sigd)); + + if (std::abs(theta) < thetaT) + { + F = SolidMath::CalI1(sig) / 3. * std::sin(phi) + + std::sqrt(SolidMath::CalJ2(sigd)) + * (std::cos(theta) - 1. / std::sqrt(3.) * std::sin(phi) * std::sin(theta)) + - coh * std::cos(phi); + } + else + F = SolidMath::CalI1(sig) / 3. * std::sin(phi) + + std::sqrt(SolidMath::CalJ2(sigd)) * (A(theta, phi) - B(theta, phi) * std::sin(3. * theta)) + - coh * std::cos(phi); + + return F; +} + +/************************************************************************** + FEMLib-Method: Minkley::DetaM_Dsigma() + Task: Derivative of Maxwell viscosity with respect to normalised deviatoric stress + Programing: + 06/2015 TN Implementation +**************************************************************************/ +KVec SolidMinkley::DetaM_Dsigma(double sig_eff, const KVec& sigd_i) +{ + KVec res; + if (sig_eff < DBL_EPSILON || std::abs(nvM) < DBL_EPSILON) + return sigd_i * 0.; + else + { + res = 3. / 2. * sigd_i * GM; // sig_eff in denominator lumped into pow function (-2 instead of -1) + res *= -etaM * mvM * nvM * std::pow(sig_eff, nvM - 2.) //etaM = etaM0/sinh(..) + / std::tanh(mvM * std::pow(sig_eff, nvM)); + return res; + } +} + +/************************************************************************** + FEMLib-Method: Minkley::CalViscoelasticResidual() + Task: Calculates the 12x1 residual vector. Implementation fully implicit only. + Programing: + 06/2015 TN Implementation +**************************************************************************/ +void SolidMinkley::CalViscoelasticResidual(const double dt, const KVec& dstrain_curr, const double e_curr, + const double e_p_curr, const KVec& stress_curr, const KVec& dstrain_Kel_curr, + const KVec& dstrain_Kel_t, const KVec& dstrain_Max_curr, + const KVec& dstrain_Max_t, const KVec& dstrain_p_curr, + Eigen::Matrix& res) +{ + const KVec dstress_curr(GM * SolidMath::P_dev * stress_curr); + + // calculate stress residual + res.block<6, 1>(0, 0) = stress_curr - (2. * (dstrain_curr - dstrain_Kel_curr - dstrain_Max_curr - dstrain_p_curr) + + KM / GM * (e_curr - e_p_curr) * SolidMath::ivec); + // calculate Kelvin strain residual + res.block<6, 1>(6, 0) = (dstrain_Kel_curr - dstrain_Kel_t) / dt + - (dstress_curr - 2. * GK0 * dstrain_Kel_curr) / (2. * etaK0); + // calculate Kelvin strain residual + res.block<6, 1>(12, 0) = (dstrain_Max_curr - dstrain_Max_t) / dt - dstress_curr / (2. * etaM); +} + +/************************************************************************** + FEMLib-Method: Minkley::CalViscoelasticJacobian() + Task: Calculates the 12x2 Jacobian matrix. Implementation fully implicit only. + Programing: + 06/2015 TN Implementation +**************************************************************************/ +void SolidMinkley::CalViscoelasticJacobian(const double dt, const KVec& stress_curr, const double sig_eff, + Eigen::Matrix& Jac) +{ + // 6x6 submatrices of the Jacobian + const KVec sigd_curr(GM * SolidMath::P_dev * stress_curr); + const KVec dmu_vM = DetaM_Dsigma(sig_eff * GM, sigd_curr); + + // Check Dimension of Jacobian + // assert(Jac.cols() == 18 && Jac.rows() == 18); + Jac.setZero(18, 18); + + // build G_11 + Jac.block<6, 6>(0, 0) = SolidMath::ident; + + // build G_12 + Jac.block<6, 6>(0, 6) = 2. * SolidMath::ident; + + // build G_13 + Jac.block<6, 6>(0, 12) = 2. * SolidMath::ident; + + // build G_21 + Jac.block<6, 6>(6, 0) = -GM / (2. * etaK0) * SolidMath::P_dev; + + // build G_22 + Jac.block<6, 6>(6, 6) = (1. / dt + GK0 / etaK0) * SolidMath::ident; + + // nothing to do for G_23 + + // build G_31 + Jac.block<6, 6>(12, 0) = -1. / (2. * etaM) * (GM * SolidMath::P_dev - sigd_curr * dmu_vM.transpose() / etaM); + + // nothing to do for G_32 + + // build G_33 + Jac.block<6, 6>(12, 12) = 1. / dt * SolidMath::ident; +} + +/************************************************************************** + FEMLib-Method: Burgers::CaldGdE() + Task: Calculates the 12x6 derivative of the residuals with respect to total strain. Implementation fully implicit +only. + Programing: + 06/2015 TN Implementation +**************************************************************************/ +void SolidMinkley::CaldGdE(Eigen::Matrix& dGdE) +{ + // Check Dimension of dGdE + // assert(dGdE.cols() == 6 && dGdE.rows() == 18); + dGdE.setZero(18, 6); + dGdE.block<6, 6>(0, 0) = -2. * SolidMath::P_dev - 3. * KM / GM * SolidMath::P_sph; +} + +/************************************************************************** + FEMLib-Method: Minkley::DG_DI1() + Task: \partial G / \partial I_1 + Programing: + 06/2015 TN Implementation +**************************************************************************/ +double SolidMinkley::DG_DI1(const double alpha) +{ + return std::sin(alpha) / 3.; +} + +/************************************************************************** + FEMLib-Method: Minkley::DG_DJ2() + Task: \partial G / \partial J2 + Programing: + 06/2015 TN Implementation +**************************************************************************/ +double SolidMinkley::DG_DJ2(const double theta, const double J2, const double alpha) +{ + if (std::abs(theta) < thetaT) + return (std::cos(theta) - std::sin(alpha) * std::sin(theta) / std::sqrt(3.)) / (2. * std::sqrt(J2)); + else + return (A(theta, alpha) - B(theta, alpha) * std::sin(3. * theta)) / (2. * std::sqrt(J2)); +} + +/************************************************************************** + FEMLib-Method: Minkley::DG_Dtheta() + Task: \partial G / \partial theta + Programing: + 06/2015 TN Implementation +**************************************************************************/ +double SolidMinkley::DG_Dtheta(const double theta, const double J2, const double alpha) +{ + if (std::abs(theta) < thetaT) + return -std::sqrt(J2) * (std::sin(theta) + std::sin(alpha) * std::cos(theta) / std::sqrt(3.)); + else + return -std::sqrt(J2) * 3. * B(theta, alpha) * std::cos(3. * theta); +} + +/************************************************************************** + FEMLib-Method: Minkley::Dtheta_DJ2() + Task: \partial theta / \partial J_2 + Programing: + 06/2015 TN Implementation +**************************************************************************/ +double SolidMinkley::Dtheta_DJ2(const double theta, const double J2) +{ + return -std::tan(3. * theta) / (2. * J2); +} + +/************************************************************************** + FEMLib-Method: Minkley::Dtheta_DJ3() + Task: \partial theta / \partial J_3 + Programing: + 06/2015 TN Implementation +**************************************************************************/ +double SolidMinkley::Dtheta_DJ3(const double theta, const double J3) +{ + return std::tan(3. * theta) / (3. * J3); +} + +/************************************************************************** + FEMLib-Method: Minkley::CalViscoplasticResidual() + Task: Calculates the 21x1 residual vector. Implementation fully implicit only. + Programing: + 06/2015 TN Implementation +**************************************************************************/ +void SolidMinkley::CalViscoplasticResidual(const double dt, const KVec& dstrain_curr, const double e_curr, + const KVec& stress_curr, const KVec& dstrain_Kel_curr, + const KVec& dstrain_Kel_t, const KVec& dstrain_Max_curr, + const KVec& dstrain_Max_t, const KVec& dstrain_pl_curr, + const KVec& dstrain_pl_t, const double e_pl_vol_curr, + const double e_pl_vol_t, const double e_pl_eff_curr, const double e_pl_eff_t, + const double lam_curr, Eigen::Matrix& res) +{ + const KVec sigd_curr(GM * SolidMath::P_dev * stress_curr); + const double J_2(SolidMath::CalJ2(sigd_curr)), J_3(SolidMath::CalJ3(sigd_curr)), + theta(SolidMath::CalLodeAngle(sigd_curr)); + const KVec dev_sigd_curr_inv(SolidMath::P_dev * SolidMath::InvertVector(sigd_curr)); + const double vol_flow(3. * DG_DI1(psi)); + + KVec dev_flow; + if (std::abs(J_3) > 0.) + dev_flow = (DG_DJ2(theta, J_2, psi) + DG_Dtheta(theta, J_2, psi) * Dtheta_DJ2(theta, J_2)) * sigd_curr + + (DG_Dtheta(theta, J_2, psi) * Dtheta_DJ3(theta, J_3) * J_3) * dev_sigd_curr_inv; + else + dev_flow.setZero(6); + + // calculate stress residual + res.block<6, 1>(0, 0) = stress_curr - (2. * (dstrain_curr - dstrain_Kel_curr - dstrain_Max_curr - dstrain_pl_curr) + + KM / GM * (e_curr - e_pl_vol_curr) * SolidMath::ivec); + + // calculate deviatoric Kelvin strain residual + res.block<6, 1>(6, 0) = (dstrain_Kel_curr - dstrain_Kel_t) / dt + - (sigd_curr - 2. * GK0 * dstrain_Kel_curr) / (2. * etaK0); + + // calculate deviatoric Maxwell strain residual + res.block<6, 1>(12, 0) = (dstrain_Max_curr - dstrain_Max_t) / dt - sigd_curr / (2. * etaM); + + // calculate deviatoric plastic strain residual + res.block<6, 1>(18, 0) = (dstrain_pl_curr - dstrain_pl_t) / dt - lam_curr * dev_flow; + + // calculate volumetric plastic strain residual + res.block<1, 1>(24, 0)(0) = (e_pl_vol_curr - e_pl_vol_t) / dt - lam_curr * vol_flow; + + // calculate effective plastic strain residual + res.block<1, 1>(25, 0)(0) = (e_pl_eff_curr - e_pl_eff_t) / dt + - std::sqrt(2. / 3. * lam_curr * lam_curr * (double)(dev_flow.transpose() * dev_flow)); + + // yield function with viscoplastic regularisation + res.block<1, 1>(26, 0)(0) = YieldMohrCoulomb(stress_curr * GM) / GM - lam_curr * eta_reg; +} + +/************************************************************************** + FEMLib-Method: Minkley::s_odot_s() + Task: vec \odot vec (all in Kelvin mapping) + Programing: + 06/2015 TN Implementation +**************************************************************************/ +/*Kelvin mapping of odot changed according to symbolic conversion. +Note: the factors of 2 and sqrt(2) come from the fact that the +incoming quantities are transformed to actual tensor coordinates +and the resulting 4th order tensor is then remapped with the +Kelvin scheme*/ +KMat SolidMinkley::s_odot_s(const KVec& vec) +{ + KMat odot; + + odot(0, 0) = -vec(0) * vec(0); + odot(0, 1) = odot(1, 0) = -vec(3) * vec(3) / 2.; + odot(0, 2) = odot(2, 0) = -vec(5) * vec(5) / 2.; + odot(0, 3) = odot(3, 0) = -vec(0) * vec(3); + odot(0, 4) = odot(4, 0) = -vec(3) * vec(5) / std::sqrt(2.); + odot(0, 5) = odot(5, 0) = -vec(0) * vec(5); + + odot(1, 1) = -vec(1) * vec(1); + odot(1, 2) = odot(2, 1) = -vec(4) * vec(4) / 2.; + odot(1, 3) = odot(3, 1) = -vec(3) * vec(1); + odot(1, 4) = odot(4, 1) = -vec(1) * vec(4); + odot(1, 5) = odot(5, 1) = -vec(3) * vec(4) / std::sqrt(2.); + + odot(2, 2) = -vec(2) * vec(2); + odot(2, 3) = odot(3, 2) = -vec(5) * vec(4) / std::sqrt(2.); + odot(2, 4) = odot(4, 2) = -vec(4) * vec(2); + odot(2, 5) = odot(5, 2) = -vec(5) * vec(2); + + odot(3, 3) = -vec(0) * vec(1) - vec(3) * vec(3) / 2.; + odot(3, 4) = odot(4, 3) = -vec(3) * vec(4) / 2. - vec(5) * vec(1) / std::sqrt(2.); + odot(3, 5) = odot(5, 3) = -vec(0) * vec(4) / std::sqrt(2.) - vec(3) * vec(5) / 2.; + + odot(4, 4) = -vec(1) * vec(2) - vec(4) * vec(4) / 2.; + odot(4, 5) = odot(5, 4) = -vec(3) * vec(2) / std::sqrt(2.) - vec(5) * vec(4) / 2.; + + odot(5, 5) = -vec(0) * vec(2) - vec(5) * vec(5) / 2.; + return odot; +} + +/************************************************************************** + FEMLib-Method: Minkley::DDG_DDJ2() + Task: \partial^2 G / \partial J_2^2 + Programing: + 06/2015 TN Implementation +**************************************************************************/ +double SolidMinkley::DDG_DDJ2(const double theta, const double J2, const double alpha) +{ + if (std::abs(theta) < thetaT) + return (std::cos(theta) - std::sin(alpha) * std::sin(theta) / std::sqrt(3.)) / (-4. * std::pow(J2, 1.5)); + else + return (A(theta, alpha) - B(theta, alpha) * std::sin(3. * theta)) / (-4. * std::pow(J2, 1.5)); +} + +/************************************************************************** + FEMLib-Method: Minkley::DDG_DJ2_Dtheta() + Task: \partial^2 G / (\partial J_2 \partial theta) + Programing: + 06/2015 TN Implementation +**************************************************************************/ +double SolidMinkley::DDG_DJ2_Dtheta(const double theta, const double J2, const double alpha) +{ + if (std::abs(theta) < thetaT) + return (std::sin(theta) + std::sin(alpha) * std::cos(theta) / std::sqrt(3.)) / (-2. * std::sqrt(J2)); + else + return 3. * B(theta, alpha) * std::cos(3. * theta) / (-2. * std::sqrt(J2)); +} + +/************************************************************************** + FEMLib-Method: Minkley::DDG_DDtheta() + Task: \partial^2 G / \partial theta^2 + Programing: + 06/2015 TN Implementation +**************************************************************************/ +double SolidMinkley::DDG_DDtheta(const double theta, const double J2, const double alpha) +{ + if (std::abs(theta) < thetaT) + return (std::cos(theta) - std::sin(alpha) * std::sin(theta) / std::sqrt(3.)) * (-std::sqrt(J2)); + else + return 9. * B(theta, alpha) * std::sin(3. * theta) * std::sqrt(J2); +} + +/************************************************************************** + FEMLib-Method: Minkley::DDtheta_DDJ2() + Task: \partial^2 theta / \partial J_2^2 + Programing: + 06/2015 TN Implementation +**************************************************************************/ +double SolidMinkley::DDtheta_DDJ2(const double theta, const double J2) +{ + return std::tan(3. * theta) / (2. * J2 * J2); +} + +/************************************************************************** + FEMLib-Method: Minkley::DDtheta_DJ2_Dtheta() + Task: \partial^2 theta / (\partial J_2 \partial theta) + Programing: + 06/2015 TN Implementation +**************************************************************************/ +double SolidMinkley::DDtheta_DJ2_Dtheta(const double theta, const double J2) +{ + return -3. / (2. * J2 * std::pow(std::cos(3. * theta), 2.)); +} + +/************************************************************************** + FEMLib-Method: Minkley::DDtheta_DJ3_Dtheta() + Task: \partial^2 theta / (\partial J_3 \partial theta) + Programing: + 06/2015 TN Implementation +**************************************************************************/ +double SolidMinkley::DDtheta_DJ3_Dtheta(const double theta, const double J3) +{ + return 1. / (J3 * std::pow(std::cos(3. * theta), 2.)); +} + +/************************************************************************** + FEMLib-Method: Minkley::DDtheta_DDJ3() + Task: \partial^2 theta / \partial J_3^2 + Programing: + 06/2015 TN Implementation +**************************************************************************/ +double SolidMinkley::DDtheta_DDJ3(const double theta, const double J3) +{ + return -std::tan(3. * theta) / (3. * J3 * J3); +} + +/************************************************************************** + FEMLib-Method: Minkley::CalViscoplasticJacobian() + Task: Calculates the 27x27 Jacobian matrix. Implementation fully implicit only. + Programing: + 06/2015 TN Implementation +**************************************************************************/ +void SolidMinkley::CalViscoplasticJacobian(const double dt, const KVec& stress_curr, const double sig_eff, + const double lam_curr, const double e_eff_i, Eigen::Matrix& Jac) +{ + // submatrices of the Jacobian + const KVec sigd_curr(GM * SolidMath::P_dev * stress_curr); + const double J_2(SolidMath::CalJ2(sigd_curr)), J_3(SolidMath::CalJ3(sigd_curr)), + theta(SolidMath::CalLodeAngle(sigd_curr)); + const KVec sigd_curr_inv(SolidMath::InvertVector(sigd_curr)); + const KVec dev_sigd_curr_inv(SolidMath::P_dev * sigd_curr_inv); + const double vol_flow(3. * DG_DI1(psi)); + const KVec dmu_vM = DetaM_Dsigma(sig_eff * GM, sigd_curr); + KVec dev_flow; + KMat Ddev_flowDsigma; + const double DthetaDJ2(Dtheta_DJ2(theta, J_2)); + const double DthetaDJ3(Dtheta_DJ3(theta, J_3)); + + // Check Dimension of Jacobian + // assert(Jac.cols() == 27 && Jac.rows() == 27); + Jac.setZero(27, 27); + + if (std::abs(J_3) < 0.) + { + dev_flow.Zero(6, 1); + Ddev_flowDsigma.Zero(6, 6); + } + else + { + const double DGDtheta(DG_Dtheta(theta, J_2, psi)); + dev_flow = (DG_DJ2(theta, J_2, psi) + DGDtheta * DthetaDJ2) * sigd_curr + + (DGDtheta * DthetaDJ3 * J_3) * dev_sigd_curr_inv; + Ddev_flowDsigma + = ((DG_DJ2(theta, J_2, psi) + DGDtheta * DthetaDJ2) * SolidMath::P_dev + + (DGDtheta * DthetaDJ3 * J_3) * SolidMath::P_dev * s_odot_s(sigd_curr_inv) * SolidMath::P_dev + + (DDG_DDJ2(theta, J_2, psi) + DDG_DJ2_Dtheta(theta, J_2, psi) * DthetaDJ2 + + DGDtheta * DDtheta_DDJ2(theta, J_2)) * sigd_curr * sigd_curr.transpose() + + (DDG_DJ2_Dtheta(theta, J_2, psi) + DDG_DDtheta(theta, J_2, psi) * DthetaDJ2 + + DGDtheta * DDtheta_DJ2_Dtheta(theta, J_2)) + * (DthetaDJ2 * sigd_curr * sigd_curr.transpose() + + DthetaDJ3 * J_3 * sigd_curr * dev_sigd_curr_inv.transpose()) + + DDG_DJ2_Dtheta(theta, J_2, psi) * DthetaDJ3 * J_3 * dev_sigd_curr_inv * sigd_curr.transpose() + + DGDtheta * J_3 * (DthetaDJ3 + DDtheta_DDJ3(theta, J_3) * J_3) * dev_sigd_curr_inv + * dev_sigd_curr_inv.transpose() + + J_3 * (DDG_DDtheta(theta, J_2, psi) * DthetaDJ3 + DGDtheta * DDtheta_DJ3_Dtheta(theta, J_3)) + * (DthetaDJ2 * dev_sigd_curr_inv * sigd_curr.transpose() + + DthetaDJ3 * J_3 * dev_sigd_curr_inv * dev_sigd_curr_inv.transpose())) * GM; + } + + const double eff_flow = std::sqrt(2. / 3. * lam_curr * lam_curr * (double)(dev_flow.transpose() * dev_flow)); + + // build G_11 + // G_66 = SolidMath::ident/dt + + Jac.block<6, 6>(0, 0) = SolidMath::ident; + + // build G_12, G_13 and G_14 + Jac.block<6, 6>(0, 6) = 2. * SolidMath::ident; + Jac.block<6, 6>(0, 12) = 2. * SolidMath::ident; + Jac.block<6, 6>(0, 18) = 2. * SolidMath::ident; + + // build G_15 + Jac.block<6, 1>(0, 24) = KM / GM * SolidMath::ivec; + + // G_16 and G_17 remain zeros and are not set separately + + // build G_21 + Jac.block<6, 6>(6, 0) = -GM / (2. * etaK0) * SolidMath::P_dev; + + // build G_22 + Jac.block<6, 6>(6, 6) = (1. / dt + GK0 / etaK0) * SolidMath::ident; + + // G_23 through G_27 are zero + + // build G_31 + Jac.block<6, 6>(12, 0) = -1. / (2. * etaM) * (GM * SolidMath::P_dev - sigd_curr * dmu_vM.transpose() / etaM); + + // G_32 is 0 + + // G_33 + Jac.block<6, 6>(12, 12) = SolidMath::ident / dt; + + // G_34 through G_37 are zero + + // G_41 + Jac.block<6, 6>(18, 0) = -lam_curr * Ddev_flowDsigma; + + // build G_42 and G_43 are zero + + // build G_44 + Jac.block<6, 6>(18, 18) = SolidMath::ident / dt; + + // G_45 and G_46 are zero + + // build G_47 + Jac.block<6, 1>(18, 26) = -dev_flow; + + // G_51 to G_54 are zero + + // build G_55 + Jac.block<1, 1>(24, 24)(0) = 1. / dt; + + // G_56 is zero + + // G_57 + Jac.block<1, 1>(24, 26)(0) = -vol_flow; + + // Conditional build of G_61 und G_67 + if (std::abs(eff_flow) > 0.) + { + Jac.block<1, 6>(25, 0) = -2. * lam_curr * lam_curr / (3. * eff_flow) * dev_flow.transpose() + * Ddev_flowDsigma.transpose(); + Jac.block<1, 1>(25, 26)(0) = -2. / (3. * eff_flow) * std::abs(lam_curr) + * (double)(dev_flow.transpose() * dev_flow); + } + // G_62 to G_64 and G_65 are zero + + // build G_66 + Jac.block<1, 1>(25, 25)(0) = 1. / dt; + + // Yield surface derivatives + if (std::abs(J_3) > 0.) + dev_flow = (DG_DJ2(theta, J_2, phi) + DG_Dtheta(theta, J_2, phi) * DthetaDJ2) * sigd_curr + + (DG_Dtheta(theta, J_2, phi) * DthetaDJ3 * J_3) * dev_sigd_curr_inv; + + else + dev_flow.Zero(6, 1); + + // build G_71 + Jac.block<1, 6>(26, 0) = (DG_DI1(phi) * SolidMath::ivec + dev_flow).transpose(); + + // G_72 - G_75 zero + + // build G_76 + Jac.block<1, 1>(26, 25)(0) = -coh0 * (hard + 2. * e_eff_i * (hard2 + 2.*e_eff_i*e_eff_i*hard4)) * std::cos(phi) / GM; + + // build G_77 + Jac.block<1, 1>(26, 26)(0) = -eta_reg; +} + +/************************************************************************** + FEMLib-Method: Burgers::CalEPdGdE() + Task: Calculates the 27x6 derivative of the residuals with respect to total strain. Implementation fully implicit +only. + Programing: + 06/2015 TN Implementation +**************************************************************************/ +void SolidMinkley::CalEPdGdE(Eigen::Matrix& dGdE) +{ + // Check Dimension of dGdE + // assert(dGdE.cols() == 6 && dGdE.rows() == 27); + dGdE.setZero(27, 6); + dGdE.block<6, 6>(0, 0) = -2. * SolidMath::P_dev - 3. * KM / GM * SolidMath::P_sph; +} + +} // namespace Minkley ends diff --git a/FEM/minkley.h b/FEM/minkley.h new file mode 100644 index 000000000..894ac234c --- /dev/null +++ b/FEM/minkley.h @@ -0,0 +1,84 @@ +/** + * \copyright + * Copyright (c) 2016, OpenGeoSys Community (http://www.opengeosys.org) + * Distributed under a Modified BSD License. + * See accompanying file LICENSE.txt or + * http://www.opengeosys.org/project/license + * + */ + +#ifndef MINKLEY_H +#define MINKLEY_H + +#include +#include "Eigen/Dense" +// FEM-Makros +#include "makros.h" +#include "tools.h" +#include "rf_msp_new.h" +#include "invariants.h" +#include "PhysicalConstant.h" + +namespace Minkley +{ +typedef Eigen::Matrix KVec; +typedef Eigen::Matrix KMat; + +class SolidMinkley +{ +public: + explicit SolidMinkley(const Math_Group::Matrix& data); + ~SolidMinkley() {} + // basic material parameters + double GK0, GM0, KM0, etaK0, etaM0, mvM; + // solution dependent values + double etaM, GM, KM; + // Parameters for Minkley model + double nvM, coh0, hard, phi, psi, thetaT, eta_reg, hard2, hard4; + // solution dependent + double coh; + double m_GM; // slope of elesticity temperature dependence + double m_KM; // slope of elesticity temperature dependence + double T_ref; // reference temperature dependency parameter for " + double Bt; // constant factor for Arrhenius term + double Q; // activation energy in Arrhenius term + + void UpdateMinkleyProperties(double s_eff, const double eps_p_eff, double Temperature); + double YieldMohrCoulomb(const KVec& sig); + void CalViscoelasticResidual(const double dt, const KVec& dstrain_curr, const double e_curr, const double e_p_curr, + const KVec& stress_curr, const KVec& dstrain_Kel_curr, const KVec& dstrain_Kel_t, + const KVec& dstrain_Max_curr, const KVec& dstrain_Max_t, const KVec& dstrain_p_curr, + Eigen::Matrix& res); + void CalViscoelasticJacobian(const double dt, const KVec& stress_curr, const double sig_eff, + Eigen::Matrix& Jac); + void CaldGdE(Eigen::Matrix& dGdE); + void CalViscoplasticResidual(const double dt, const KVec& dstrain_curr, const double e_curr, + const KVec& stress_curr, const KVec& dstrain_Kel_curr, const KVec& dstrain_Kel_t, + const KVec& dstrain_Max_curr, const KVec& dstrain_Max_t, const KVec& dstrain_pl_curr, + const KVec& dstrain_pl_t, const double e_pl_vol_curr, const double e_pl_vol_t, + const double e_pl_eff_curr, const double e_pl_eff_t, const double lam_curr, + Eigen::Matrix& res); + void CalViscoplasticJacobian(const double dt, const KVec& stress_curr, const double sig_eff, const double lam_curr, + const double e_eff_i, Eigen::Matrix& Jac); + void CalEPdGdE(Eigen::Matrix& dGdE); + +private: + double A(const double theta, const double alpha); + double B(const double theta, const double alpha); + KVec DetaM_Dsigma(double sig_eff, const KVec& sigd_i); + double DG_DI1(const double alpha); + double DG_DJ2(const double theta, const double J2, const double alpha); + double DG_Dtheta(const double theta, const double J2, const double alpha); + double Dtheta_DJ2(const double theta, const double J2); + double Dtheta_DJ3(const double theta, const double J3); + KMat s_odot_s(const KVec& vec); + double DDG_DDJ2(const double theta, const double J2, const double alpha); + double DDG_DJ2_Dtheta(const double theta, const double J2, const double alpha); + double DDG_DDtheta(const double theta, const double J2, const double alpha); + double DDtheta_DDJ2(const double theta, const double J2); + double DDtheta_DJ2_Dtheta(const double theta, const double J2); + double DDtheta_DJ3_Dtheta(const double theta, const double J3); + double DDtheta_DDJ3(const double theta, const double J3); +}; +} +#endif // MINKLEY_H diff --git a/FEM/pcs_dm.cpp b/FEM/pcs_dm.cpp index 020621d0d..a7296a287 100644 --- a/FEM/pcs_dm.cpp +++ b/FEM/pcs_dm.cpp @@ -772,19 +772,13 @@ double CRFProcessDeformation::Execute(int loop_process_number) Error = Norm / InitialNorm; ErrorU = NormU / InitialNormU0; - if (Norm < Tolerance_global_Newton && Error > Norm) - Error = Norm; - // if(Norm1.0e-1) damping=0.5; - if (Error / Error1 > 1.0e-1 || ErrorU / ErrorU1 > 1.0e-1) - damping = 0.5; - if (ErrorU < Error) - Error = ErrorU; + if (Error / Error1 > m_num->newton_damping_tolerance + || ErrorU / ErrorU1 > m_num->newton_damping_tolerance) + damping = m_num->newton_damping_factor; + #if defined(NEW_EQS) && defined(JFNK_H2M) /// If JFNK, get w from the buffer if (m_num->nls_method == 2) @@ -811,27 +805,30 @@ double CRFProcessDeformation::Execute(int loop_process_number) if (myrank == 0) { #endif - //Screan printing: - std::cout<<" -->End of Newton-Raphson iteration: "<End of Newton-Raphson iteration: " << ite_steps << "/" << MaxIteration + << "\n"; + cout.width(8); + cout.precision(2); + cout.setf(ios::scientific); + std::cout << " DeltaU/DeltaU0" + << " " + << "DeltaU" + << " " + << "DeltaF/DeltaF0 " + << " " + << "DeltaF" + << " " + << "Damping" + << "\n"; + std::cout << " " << ErrorU << " " << NormU << " " << Error << " " << Norm << " " + << " " << damping << "\n"; + std::cout << " ------------------------------------------------" + << "\n"; #if defined(USE_MPI) || defined(USE_PETSC) } #endif - if (Error > 100.0 && ite_steps > 1) - { - printf("\n Attention: Newton-Raphson step is diverged. Programme halt!\n"); - exit(1); - } - if (InitialNorm < 10 * Tolerance_global_Newton) - break; - if (Norm < 0.001 * InitialNorm) - break; - if (Error <= Tolerance_global_Newton) + if (Error <= Tolerance_global_Newton)//testing only relative residual so far { if (ite_steps == 1) // WX:05.2012 { @@ -842,6 +839,7 @@ double CRFProcessDeformation::Execute(int loop_process_number) break; } } + // w = w+dw for Newton-Raphson UpdateIterativeStep(damping, 0); // w = w+dw } // Newton-Raphson iteration diff --git a/FEM/rf_msp_new.cpp b/FEM/rf_msp_new.cpp index e0c4bbeac..edc15190b 100644 --- a/FEM/rf_msp_new.cpp +++ b/FEM/rf_msp_new.cpp @@ -23,6 +23,8 @@ //#include #include +#include "Eigen/Dense" + // FEM-Makros #include "makros.h" #include "rf_pcs.h" @@ -40,6 +42,9 @@ #include "tools.h" // GetLineFromFile #include "PhysicalConstant.h" +#include "minkley.h" +#include "burgers.h" + std::vector msp_vector; std::vector msp_key_word_vector; // OK @@ -575,6 +580,73 @@ std::ios::pos_type CSolidProperties::Read(std::ifstream* msp_file) in_sd >> (*data_Creep)(i, 0); in_sd.clear(); } + if (line_string.find("BURGERS") != string::npos) + { + Creep_mode = 1001; + // data_Creep: + // 0: G_K0 + // 1: m_K + // 2: mu_K + // 3: m_vK + // 4: G_M0 + // 5: K_M0 + // 6: mu_M + // 7: m_vM + // 8: m_GM // slope of elesticity temperature dependence of G_M + // 9: m_KM // slope of elesticity temperature dependence of K_M + // 10: T_ref // reference temperature; + // 11: B // constant factor for Arrhenius term + // 12: Q // activation energy in Arrhenius term + data_Creep = new Matrix(13); + in_sd.str(GetLineFromFile1(msp_file)); + for (i = 0; i < 13; i++) + { + (*data_Creep)(i) = 0.; + in_sd >> (*data_Creep)(i); + } + in_sd.clear(); + + // Local Newton scheme for Burgers model. TN 06.06.2014 + // initialised fully 3D + material_burgers = new Burgers::SolidBurgers(*data_Creep); + } + if (line_string.find("MINKLEY") != string::npos) + { + Creep_mode = 1002; + // data_Creep: + // 0: Kelvin shear modulus + // 1: Kelvin shear viscosity + // 2: Maxwell shear modulus + // 3: Maxwell bulk modulus + // 4: Maxwell shear viscosity + // 5: m -- effective stress factor in viscosity relationship + // 6: n -- effective stress exponent in viscosity relationship + // 7: initial cohesion + // 8: hardening/softening modulus + // 9: friction angle + // 10: dilatancy angle + // 11: transition angle for corner smoothing + // 12: viscosity for viscoplastic regularisation + // 13: m_GM // slope of elesticity temperature dependence of G_M + // 14: m_KM // slope of elesticity temperature dependence of K_M + // 15: T_ref // reference temperature; + // 16: B // constant factor for Arrhenius term + // 17: Q // activation energy in Arrhenius term + // 18: h2 // second order hardening term + // 19: h4 // fourth order hardening term + + data_Creep = new Matrix(20); + in_sd.str(GetLineFromFile1(msp_file)); + for (i = 0; i < 20; i++){ + (*data_Creep)(i) = 0.; + in_sd >> (*data_Creep)(i); + } + in_sd.clear(); + + // Local Newton scheme for Burgers model. TN 06.06.2014 + // initialised fully 3D + material_minkley = new Minkley::SolidMinkley(*data_Creep); + } } // WX:10.2012, threshold dev. stress for Lubby2 if (line_string.find("$THRESHOLD_DEV_STR") != string::npos) @@ -642,6 +714,14 @@ std::ios::pos_type CSolidProperties::Read(std::ifstream* msp_file) in_sd >> grav_const; in_sd.clear(); } + if (line_string.find("$GRAVITY_RAMP") != string::npos) + { + in_sd.str(GetLineFromFile1(msp_file)); + in_sd >> grav_curve_id; + in_sd >> grav_const; + in_sd.clear(); + gravity_ramp = 1; + } //.................................................................... // subkeyword found if (line_string.find("$PLASTICITY") != string::npos) @@ -945,6 +1025,7 @@ CSolidProperties::CSolidProperties() T_0 = 0.0; Creep_mode = -1; grav_const = 9.81; // WW + gravity_ramp = 0; excavation = -1; // 12.2009. WW excavated = false; // To be ..... 12.2009. WW E_Function_Model = -1; // WX:06.2012 E dependence @@ -1025,6 +1106,8 @@ CSolidProperties::CSolidProperties() reaction_entropy = 0.0; non_reactive_solid_volume_fraction = 0.0; non_reactive_solid_density = 0.0; + material_minkley = NULL; + material_burgers = NULL; T_ref_enthalpy_correction = 573.15; specific_heat_source = 0.0; @@ -1162,6 +1245,9 @@ CSolidProperties::~CSolidProperties() rhs_l = NULL; x_l = NULL; Li = NULL; + material_minkley = NULL; + material_burgers = NULL; + smath = NULL; } //---------------------------------------------------------------------------- @@ -1579,6 +1665,304 @@ double CSolidProperties::Kronecker(const int ii, const int jj) return delta; } +/************************************************************************** + FEMLib-Method: CSolidProperties::ExtractConsistentTangent() + Task: general routine to get 6x6 consistent tangent from local Newton iteration of material functionals + Programing: + 06/2014 TN Implementation +**************************************************************************/ +void CSolidProperties::ExtractConsistentTangent(const Eigen::MatrixXd& Jac, const Eigen::MatrixXd& dGdE, + const bool pivoting, Eigen::Matrix& dsigdE) +{ + const unsigned int local_dim(Jac.cols()); + // Check Dimensions + assert(local_dim == dGdE.rows() && dGdE.cols() == 6); + + Eigen::MatrixXd dzdE(local_dim, 6); + // solve linear system + if (pivoting) + dzdE = Jac.fullPivHouseholderQr().solve(-1.0 * dGdE); // Could consider moving to different Eigen solver. + //dzdE = Jac.fullPivLu().solve(-1.0 * dGdE); // Could consider moving to different Eigen solver. + else + dzdE = Jac.householderQr().solve(-1.0 * dGdE); // Could consider moving to different Eigen solver. + // in-built Gauss elimination solver was at least 4 OoM more inaccurate. + + // Extract matrix part relevant for global tangent + dsigdE = dzdE.block<6, 6>(0, 0); +} + +/************************************************************************** + FEMLib-Method: CSolidProperties::LocalNewtonBurgers() + Task: general local Newton routine to integrate inelastic material models + Programing: + 06/2014 TN Implementation + 03/2015 NB Modified +**************************************************************************/ +void CSolidProperties::LocalNewtonBurgers(const double dt, const std::vector& strain_curr, + std::vector& stress_curr, std::vector& strain_K_curr, + std::vector& strain_M_curr, Math_Group::Matrix& Consistent_Tangent, + double Temperature, double& local_res) +{ + // stress, strain, internal variable + KVec sig_j, eps_K_j, eps_M_j; + // deviatoric stress, strain + KVec sigd_j; + // local residual vector and Jacobian + Eigen::Matrix res_loc, inc_loc; + Eigen::Matrix K_loc; + double sig_eff(1.); + + // initialisation of Kelvin vectors + // Note: Can be done in one loop instead of 5 if done right here. + const KVec eps_i(SolidMath::Voigt_to_Kelvin_Strain(strain_curr)); + const KVec eps_K_t(SolidMath::Voigt_to_Kelvin_Strain(strain_K_curr)); + const KVec eps_M_t(SolidMath::Voigt_to_Kelvin_Strain(strain_M_curr)); + eps_M_j = eps_M_t; + eps_K_j = eps_K_t; + + if (!T_Process) + Temperature = material_burgers->T_ref; + + // calculation of deviatoric and spherical parts + const double e_i = eps_i(0) + eps_i(1) + eps_i(2); + const KVec epsd_i(SolidMath::P_dev * eps_i); + + // dimensionless stresses + sigd_j = 2.0 * (epsd_i - eps_M_t - eps_K_t); // initial guess as elastic predictor + + // Calculate effective stress and update material properties + sig_eff = SolidMath::CalEffectiveStress(sigd_j); + + material_burgers->UpdateBurgersProperties(sig_eff, Temperature); + + // initial evaluation of residual and Jacobian + material_burgers->CalResidualBurgers(dt, epsd_i, sigd_j, eps_K_j, eps_K_t, eps_M_j, eps_M_t, res_loc); + // initial evaluation of Jacobian + material_burgers->CalJacobianBurgers( + dt, K_loc, sig_eff, sigd_j, + eps_K_j); // Note - With constant properties a single evaluation would be sufficient. + + // Loop variables + int counter = 0; + const int counter_max(20); + + // for (int counter(0); counter local_tolerance; ++counter) + while (res_loc.norm() > Tolerance_Local_Newton && counter < counter_max) + { + counter++; + // Solve linear system + // inc_loc = K_loc.fullPivHouseholderQr().solve(res_loc); //other linear solvers (faster but less accurate) can + // be considered. + // inc_loc = K_loc.fullPivLu().solve(res_loc); //other linear solvers (faster but less accurate) can be + // considered + // inc_loc = K_loc.colPivHouseholderQr().solve(res_loc); + inc_loc = K_loc.householderQr().solve(-res_loc); + // increment solution vectors + sigd_j += inc_loc.block<6, 1>(0, 0); + eps_K_j += inc_loc.block<6, 1>(6, 0); + eps_M_j += inc_loc.block<6, 1>(12, 0); + // Calculate effective stress and update material properties + sig_eff = SolidMath::CalEffectiveStress(sigd_j); + material_burgers->UpdateBurgersProperties(sig_eff, Temperature); + // evaluation of new residual + material_burgers->CalResidualBurgers(dt, epsd_i, sigd_j, eps_K_j, eps_K_t, eps_M_j, eps_M_t, res_loc); + // Get Jacobian + material_burgers->CalJacobianBurgers(dt, K_loc, sig_eff, sigd_j, eps_K_j); // for solution dependent Jacobians + } + // if (counter == counter_max) + // std::cout << "WARNING: Maximum iteration number needed in LocalNewtonBurgers. Convergence not guaranteed." + // << std::endl; + local_res = res_loc.norm(); + + // dGdE matrix and dsigdE matrix + Eigen::Matrix dGdE; + Eigen::Matrix dsigdE; + + // Calculate dGdE for time step + material_burgers->CaldGdEBurgers(dGdE); + // get dsigdE matrix + ExtractConsistentTangent(K_loc, dGdE, false, dsigdE); + + // add hydrostatic part to stress and tangent + sig_j = material_burgers->GM * sigd_j + material_burgers->KM * e_i * SolidMath::ivec; + dsigdE = material_burgers->GM * dsigdE * SolidMath::P_dev + 3. * material_burgers->KM * SolidMath::P_sph; + + // Sort into Consistent Tangent matrix for global Newton iteration and into standard OGS arrays + SolidMath::Kelvin_to_Voigt_Stress(sig_j, stress_curr); + SolidMath::Kelvin_to_Voigt_Strain(eps_K_j, strain_K_curr); + SolidMath::Kelvin_to_Voigt_Strain(eps_M_j, strain_M_curr); + for (size_t i = 0; i < 3; i++) + { + for (size_t j = 0; j < 3; j++) + { + Consistent_Tangent(i, j) = dsigdE(i, j); + Consistent_Tangent(i, j + 3) = dsigdE(i, j + 3) + / sqrt(2.); // from local to global shear components (Kelvin to Voigt) + Consistent_Tangent(i + 3, j) = dsigdE(i + 3, j) + / sqrt(2.); // from local to global shear components (Kelvin to Voigt) + Consistent_Tangent(i + 3, j + 3) = dsigdE(i + 3, j + 3) + / 2.; // from local to global shear components (Kelvin to Voigt) + } + } +} + +/************************************************************************** + FEMLib-Method: CSolidProperties::LocalNewtonMinkley() + Task: general local Newton routine to integrate inelastic material models + Programing: + 06/2015 TN Implementation +**************************************************************************/ +void CSolidProperties::LocalNewtonMinkley(const double dt, const std::vector& strain_curr, + std::vector& stress_curr, std::vector& eps_K_curr, + std::vector& eps_M_curr, std::vector& eps_pl_curr, + double& e_pl_v, double& e_pl_eff, double& lam, + Math_Group::Matrix& Consistent_Tangent, double Temperature, + double& local_res) +{ + // stress, strain, internal variable + KVec sig_j, eps_K_j, eps_M_j, eps_pl_j; + // deviatoric stress, strain + KVec sigd_j; + const double e_pl_v_t = e_pl_v, e_pl_eff_t = e_pl_eff; + // local residual vector and Jacobian + Eigen::Matrix res_loc, inc_loc; + Eigen::Matrix K_loc; + double sig_eff; + Eigen::Matrix dsigdE; + + // initialisation of Kelvin vectors + // Note: Can be done in one loop instead of 5 if done right here. + const KVec eps_i(SolidMath::Voigt_to_Kelvin_Strain(strain_curr)); + const KVec eps_K_t(SolidMath::Voigt_to_Kelvin_Strain(eps_K_curr)); + const KVec eps_M_t(SolidMath::Voigt_to_Kelvin_Strain(eps_M_curr)); + const KVec eps_pl_t(SolidMath::Voigt_to_Kelvin_Strain(eps_pl_curr)); + eps_M_j = eps_M_t; + eps_K_j = eps_K_t; + eps_pl_j = eps_pl_t; + + if (!T_Process) + Temperature = material_minkley->T_ref; + + // calculation of deviatoric and spherical parts + const double e_i = eps_i(0) + eps_i(1) + eps_i(2); + const KVec epsd_i(SolidMath::P_dev * eps_i); + + // dimensionless stresses + sigd_j = 2.0 * (epsd_i - eps_M_j - eps_K_j - eps_pl_j); // initial guess as elastic predictor + // Calculate effective stress and update material properties + sig_eff = SolidMath::CalEffectiveStress(sigd_j); + material_minkley->UpdateMinkleyProperties(sig_eff, e_pl_eff, Temperature); + sig_j = sigd_j + material_minkley->KM / material_minkley->GM * (e_i - e_pl_v) * SolidMath::ivec; + // std::cout << "Stress guesstimate_zz " << sig_j(2)*material_minkley->GM0 << std::endl; + + // initial evaluation of residual and Jacobian + material_minkley->CalViscoelasticResidual(dt, epsd_i, e_i, e_pl_v, sig_j, eps_K_j, eps_K_t, eps_M_j, eps_M_t, + eps_pl_j, res_loc); + // initial evaluation of Jacobian + material_minkley->CalViscoelasticJacobian(dt, sig_j, sig_eff, K_loc); + + // Loop variables + int counter = 0; + const int counter_max(20); + + while (res_loc.norm() > Tolerance_Local_Newton && counter < counter_max) + { + counter++; + // Solve linear system; Choice of solver can be influenced by material properties/property ratios + //inc_loc = K_loc.fullPivHouseholderQr().solve(-res_loc); + //others: fullPivLu() and colPivHouseholderQr() + inc_loc = K_loc.fullPivLu().solve(-res_loc); + // increment solution vectors + sig_j += inc_loc.block<6, 1>(0, 0); + eps_K_j += inc_loc.block<6, 1>(6, 0); + eps_M_j += inc_loc.block<6, 1>(12, 0); + // Calculate effective stress and update material properties + sig_eff = SolidMath::CalEffectiveStress(SolidMath::P_dev * sig_j); + material_minkley->UpdateMinkleyProperties(sig_eff, e_pl_eff, Temperature); + // evaluation of new residual + material_minkley->CalViscoelasticResidual(dt, epsd_i, e_i, e_pl_v, sig_j, eps_K_j, eps_K_t, eps_M_j, eps_M_t, + eps_pl_j, res_loc); + // Get Jacobian + material_minkley->CalViscoelasticJacobian(dt, sig_j, sig_eff, K_loc); + } + if (!(material_minkley->YieldMohrCoulomb(sig_j * material_minkley->GM) < Tolerance_Local_Newton)) + { + Eigen::Matrix res_loc_p, inc_loc_p; + Eigen::Matrix K_loc_p; + material_minkley->CalViscoplasticResidual(dt, epsd_i, e_i, sig_j, eps_K_j, eps_K_t, eps_M_j, eps_M_t, eps_pl_j, + eps_pl_t, e_pl_v, e_pl_v_t, e_pl_eff, e_pl_eff_t, lam, res_loc_p); + material_minkley->CalViscoplasticJacobian(dt, sig_j, sig_eff, lam, e_pl_eff, K_loc_p); + while (res_loc_p.norm() > Tolerance_Local_Newton && counter < 2*counter_max) + { + counter++; + // Solve linear system; Choice of solver can be influenced by material properties/property ratios + inc_loc_p = K_loc_p.fullPivLu().solve(-res_loc_p); + //inc_loc_p = K_loc_p.householderQr().solve(-res_loc_p); + // increment solution vectors + sig_j += inc_loc_p.block<6, 1>(0, 0); + eps_K_j += inc_loc_p.block<6, 1>(6, 0); + eps_M_j += inc_loc_p.block<6, 1>(12, 0); + eps_pl_j += inc_loc_p.block<6, 1>(18, 0); + e_pl_v += inc_loc_p.block<1, 1>(24, 0)(0); + e_pl_eff += inc_loc_p.block<1, 1>(25, 0)(0); + lam += inc_loc_p.block<1, 1>(26, 0)(0); + // Calculate effective stress and update material properties + sig_eff = SolidMath::CalEffectiveStress(SolidMath::P_dev * sig_j); + material_minkley->UpdateMinkleyProperties(sig_eff, e_pl_eff, Temperature); + // evaluation of new residual + material_minkley->CalViscoplasticResidual(dt, epsd_i, e_i, sig_j, eps_K_j, eps_K_t, eps_M_j, eps_M_t, + eps_pl_j, eps_pl_t, e_pl_v, e_pl_v_t, e_pl_eff, e_pl_eff_t, lam, + res_loc_p); + // Get Jacobian + material_minkley->CalViscoplasticJacobian(dt, sig_j, sig_eff, lam, e_pl_eff, K_loc_p); + } + // dGdE matrix and dsigdE matrix + Eigen::Matrix dGdE; + SolidMath::Kelvin_to_Voigt_Strain(eps_pl_j, eps_pl_curr); + // Calculate dGdE for time step + material_minkley->CalEPdGdE(dGdE); + // get dsigdE matrix + ExtractConsistentTangent(K_loc_p, dGdE, true, dsigdE); // Full pivoting needed for global convergence + local_res = res_loc_p.norm(); + } + else + { + // dGdE matrix and dsigdE matrix + Eigen::Matrix dGdE; + // Calculate dGdE for time step + material_minkley->CaldGdE(dGdE); + // get dsigdE matrix + ExtractConsistentTangent(K_loc, dGdE, true, dsigdE); // Full pivoting needed for global convergence + // if (counter == counter_max) + // std::cout << "WARNING: Maximum iteration number needed in LocalNewtonMinkley. Convergence not + // guaranteed." + // << std::endl; + // << std::endl;local_res = res_loc.norm(); + local_res = -res_loc.norm(); //negative to imply that viscoelastic routine was used. + } + // add hydrostatic part to stress and tangent + sig_j *= material_minkley->GM; + dsigdE *= material_minkley->GM; + // Sort into Consistent Tangent matrix for global Newton iteration and into standard OGS arrays + SolidMath::Kelvin_to_Voigt_Stress(sig_j, stress_curr); + SolidMath::Kelvin_to_Voigt_Strain(eps_K_j, eps_K_curr); + SolidMath::Kelvin_to_Voigt_Strain(eps_M_j, eps_M_curr); + // plastic strain dealt with further up + for (size_t i = 0; i < 3; i++) + { + for (size_t j = 0; j < 3; j++) + { + Consistent_Tangent(i, j) = dsigdE(i, j); + Consistent_Tangent(i, j + 3) = dsigdE(i, j + 3) + / sqrt(2.); // from local to global shear components (Kelvin to Voigt) + Consistent_Tangent(i + 3, j) = dsigdE(i + 3, j) + / sqrt(2.); // from local to global shear components (Kelvin to Voigt) + Consistent_Tangent(i + 3, j + 3) = dsigdE(i + 3, j + 3) + / 2.; // from local to global shear components (Kelvin to Voigt) + } + } +} + /************************************************************************** FEMLib-Method: CSolidProperties::Calculate_Lame_Constant() Task: Get density diff --git a/FEM/rf_msp_new.h b/FEM/rf_msp_new.h index 466592a93..ade669edb 100644 --- a/FEM/rf_msp_new.h +++ b/FEM/rf_msp_new.h @@ -22,6 +22,8 @@ last modified: //#include //#include +#include "invariants.h" + #define MSP_FILE_EXTENSION ".msp" extern bool MSPRead(const std::string& given_file_base_name); @@ -44,9 +46,23 @@ namespace process { class CRFProcessDeformation; } +namespace Minkley +{ +class SolidMinkley; +} + +namespace Burgers +{ +class SolidBurgers; +} +namespace SolidMath +{ +class Invariants; +} namespace SolidProp { +typedef Eigen::Matrix KVec; class CSolidProperties { public: @@ -136,6 +152,17 @@ class CSolidProperties // Set value for solid reactive system - TN void setSolidReactiveSystem(FiniteElement::SolidReactiveSystem reactive_system); + // general routine to get consistent tangent from local Newton iteration of material functionals + void ExtractConsistentTangent(const Eigen::MatrixXd& Jac, const Eigen::MatrixXd& dGdE, const bool pivoting, + Eigen::Matrix& dsigdE); + // general local Newton routines to integrate inelastic material models + void LocalNewtonBurgers(const double dt, const std::vector& strain_curr, std::vector& stress_curr, + std::vector& strain_K_curr, std::vector& strain_M_curr, + Math_Group::Matrix& Consistent_Tangent, double Temperature, double& local_res); + void LocalNewtonMinkley(const double dt, const std::vector& strain_curr, std::vector& stress_curr, + std::vector& eps_K_curr, std::vector& eps_M_curr, + std::vector& eps_pl_curr, double& e_pl_v, double& e_pl_eff, double& lam, + Math_Group::Matrix& Consistent_Tangent, double Temperature, double& local_res); private: // CMCD FiniteElement::CFiniteElementStd* Fem_Ele_Std; @@ -158,6 +185,8 @@ class CSolidProperties double bishop_model_value; // 05.2011 WX double threshold_dev_str; // 12.2012 WX double grav_const; // WW + int grav_curve_id; // NB + bool gravity_ramp; // NB Math_Group::Matrix* data_Density; // Math_Group::Matrix* data_Capacity; @@ -391,6 +420,9 @@ class CSolidProperties // Get solid reactive system - TN FiniteElement::SolidReactiveSystem getSolidReactiveSystem() const; + Minkley::SolidMinkley* material_minkley; + Burgers::SolidBurgers* material_burgers; + SolidMath::Invariants* smath; FiniteElement::SolidReactiveSystem _reactive_system; // Friends that can access to this data explicitly diff --git a/FEM/rf_num_new.cpp b/FEM/rf_num_new.cpp index 2cf35d8cb..cebf37e6e 100644 --- a/FEM/rf_num_new.cpp +++ b/FEM/rf_num_new.cpp @@ -126,6 +126,8 @@ CNumerics::CNumerics(string name) fct_method = -1; // NW fct_prelimiter_type = 0; // NW fct_const_alpha = -1.0; // NW + newton_damping_factor = 1.0; + newton_damping_tolerance = 1.e3; //---------------------------------------------------------------------- // Deformation GravityProfile = 0; @@ -384,6 +386,8 @@ ios::pos_type CNumerics::Read(ifstream* num_file) { line.str(GetLineFromFile1(num_file)); line >> nls_plasticity_local_tolerance; + line.clear(); + continue; } //.................................................................... // subkeyword found ($NON_LINEAR_ITERATION -or- $NON_LINEAR_ITERATIONS) @@ -698,6 +702,18 @@ ios::pos_type CNumerics::Read(ifstream* num_file) << "\n"; continue; } + // Automatic damping of Newton scheme + if (line_string.find("$NEWTON_DAMPING") != string::npos) + { + line.str(GetLineFromFile1(num_file)); + line >> newton_damping_tolerance; //if NR error decreases by less than this factor, the next step will be dampened + line >> newton_damping_factor; //dampened by this factor + line.clear(); + std::cout << "NR step will be damped by " << newton_damping_factor + << " if relative residual or relative unknown increment decrease by less than " + << newton_damping_tolerance << " from one iteration to the next." << std::endl; + continue; + } //.................................................................... /* diff --git a/FEM/rf_num_new.h b/FEM/rf_num_new.h index b5310018b..fca14a8cf 100644 --- a/FEM/rf_num_new.h +++ b/FEM/rf_num_new.h @@ -103,6 +103,8 @@ class CNumerics double lag_min_weight; int lag_use_matrix; int lag_vel_method; + double newton_damping_factor; + double newton_damping_tolerance; // // Configure void NumConfigure(bool overall_coupling_exists); // JT2012 diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index e53299341..05923223e 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -8,6 +8,7 @@ add_subdirectory( gtest ) set ( SOURCES testrunner.cpp testBase.cpp + testSolidProps.cpp ) # Add tests here if they need testdata diff --git a/tests/testSolidProps.cpp b/tests/testSolidProps.cpp new file mode 100644 index 000000000..ddd798c4b --- /dev/null +++ b/tests/testSolidProps.cpp @@ -0,0 +1,750 @@ +/** + * \copyright + * Copyright (c) 2016, OpenGeoSys Community (http://www.opengeosys.org) + * Distributed under a Modified BSD License. + * See accompanying file LICENSE.txt or + * http://www.opengeosys.org/project/license + * + */ + +// ** INCLUDES ** +//#include "stdio.h" +#include + +#include "Eigen/Dense" + +#include "gtest.h" +#include "fem_ele_vec.h" +#include "rf_msp_new.h" +#include "makros.h" +#include +#include "burgers.h" +#include "minkley.h" +#include "invariants.h" + +TEST(SolidProps, EffectiveStress) +{ + // stress vectors + Eigen::Matrix sig, sigd; + sigd.setZero(6); + sig.setZero(6); + + // stress values + std::vector stress(6); + stress[0] = 0.01; + stress[1] = 0.02; + stress[2] = 0.03; + stress[3] = 0.04; stress[4] = 0.05; stress[5] = 0.06; + sig = SolidMath::Voigt_to_Kelvin_Stress(stress); // Kelvin mapping + sigd = SolidMath::P_dev * sig; + + const double s_eff(SolidMath::CalEffectiveStress(sig)); + const double sd_eff(SolidMath::CalEffectiveStress(sigd)); + + ASSERT_NEAR(s_eff, 0.158745078664, 1.e-10); + ASSERT_NEAR(sd_eff, 0.152970585408, 1.e-10); +} + +TEST(SolidProps, LodeAngle) +{ + // stress vectors + Eigen::Matrix sig, sigd; + sigd.setZero(6); + sig.setZero(6); + + // uniaxial tension + sig(0) = 2.; + sigd = SolidMath::P_dev * sig; + + double theta(SolidMath::CalLodeAngle(sigd)); + ASSERT_NEAR(theta, -30. * PI / 180., 1.e-8); + + // uniaxial compression + sig(0) = -2.; + sigd = SolidMath::P_dev * sig; + + theta = SolidMath::CalLodeAngle(sigd); + ASSERT_NEAR(theta, 30. * PI / 180., 1.e-8); + + // equibiaxial tension + sig(0) = 2.; + sig(2) = 2.; + sigd = SolidMath::P_dev * sig; + + theta = SolidMath::CalLodeAngle(sigd); + ASSERT_NEAR(theta, 30. * PI / 180., 1.e-8); + + // biaxial tension-compression + sig(0) = -2.; + sig(2) = 2.; + sigd = SolidMath::P_dev * sig; + + theta = SolidMath::CalLodeAngle(sigd); + ASSERT_NEAR(theta, 0., 1.e-8); + + // hydrostatic loading + sig(0) = 2.; + sig(1) = 2.; + sig(2) = 2.; + sigd = SolidMath::P_dev * sig; + + theta = SolidMath::CalLodeAngle(sigd); + ASSERT_NEAR(theta, 0., 1.e-8); +} + +TEST(SolidProps, Lubby2JacobianNumeric) +{ + Math_Group::Matrix* data; + data = new Math_Group::Matrix(13); + + // set Constants + (*data)(0) = 2.0; // Kelvin shear modulus + (*data)(1) = 10.0; // dependency parameter for " + (*data)(2) = 2.0; // Kelvin viscosity + (*data)(3) = 10.0; // dependency parameter for " + (*data)(4) = 1.7; // Maxwell shear modulus + (*data)(5) = 2.0; // Maxwell bulk modulus + (*data)(6) = 10.0; // Maxwell viscosity + (*data)(7) = 10.0; // dependency parameter for " + (*data)(8) = 0.; // slope of elesticity temperature dependence + (*data)(9) = 0.; // slope of elesticity temperature dependence + (*data)(10) = 273.; // reference temperature dependency parameter for " + (*data)(11) = 1.; // constant factor for Arrhenius term + (*data)(12) = 0.; // activation energy in Arrhenius term + Burgers::SolidBurgers material(*data); + + // state and trial variables + Eigen::Matrix eps_i, eps_K_t, eps_K_j, eps_M_t, eps_M_j; + Eigen::Matrix sigd_t, sigd_j, epsd_i, epsd_t; + eps_i.setZero(6); + eps_K_t.setZero(6); + eps_K_j.setZero(6); + eps_M_t.setZero(6); + eps_M_j.setZero(6); + sigd_t.setZero(6); + sigd_j.setZero(6); + epsd_t.setZero(6); + epsd_i.setZero(6); + + // set a strain increment and perform Kelvin mapping + std::vector strain_bc(6); + strain_bc[0] = 0.01; + strain_bc[1] = 0.02; + strain_bc[2] = 0.03; + strain_bc[3] = 0.04 * 2.; + strain_bc[4] = 0.05 * 2.; + strain_bc[5] = 0.06 * 2.; + eps_i = SolidMath::Voigt_to_Kelvin_Strain(strain_bc); + + epsd_i = SolidMath::P_dev * eps_i; + + // guess stress increment (dimensionless) + sigd_j = 2.0 * epsd_i; + + // Update Material parameters + material.UpdateBurgersProperties(SolidMath::CalEffectiveStress(sigd_j), 273.); + + // set nontrivial internal variables + eps_K_j = 0.1 * epsd_i; + eps_M_j = 0.15 * epsd_i; + + Eigen::Matrix Jacobian; + Jacobian.setZero(18, 18); + + // Calculate Jacobian for time step + const double dt(0.01); + material.CalJacobianBurgers(dt, Jacobian, SolidMath::CalEffectiveStress(sigd_j), sigd_j, eps_K_j); + + // Residual vector + Eigen::Matrix residual; + + // Numerically calculate Jacobian + Eigen::Matrix Jac_num; + Jac_num.setZero(18, 18); + + Eigen::Matrix upper, lower; + const double pertub(1.e-6); + double up, low; + + for (size_t i = 0; i < 18; i++) + for (size_t j = 0; j < 18; j++) + { + if (j < 6) + { + upper = sigd_j; + lower = sigd_j; + upper(j) += pertub; + lower(j) -= pertub; + material.UpdateBurgersProperties(SolidMath::CalEffectiveStress(upper), 273.); + material.CalResidualBurgers(dt, epsd_i, upper, eps_K_j, eps_K_t, eps_M_j, eps_M_t, residual); + up = residual(i); + material.UpdateBurgersProperties(SolidMath::CalEffectiveStress(lower), 273.); + material.CalResidualBurgers(dt, epsd_i, lower, eps_K_j, eps_K_t, eps_M_j, eps_M_t, residual); + low = residual(i); + } + else if (j < 12) + { + upper = eps_K_j; + lower = eps_K_j; + upper(j - 6) += pertub; + lower(j - 6) -= pertub; + // Update Material parameters + material.UpdateBurgersProperties(SolidMath::CalEffectiveStress(sigd_j), 273.); + material.CalResidualBurgers(dt, epsd_i, sigd_j, upper, eps_K_t, eps_M_j, eps_M_t, residual); + up = residual(i); + material.CalResidualBurgers(dt, epsd_i, sigd_j, lower, eps_K_t, eps_M_j, eps_M_t, residual); + low = residual(i); + } + else + { + upper = eps_M_j; + lower = eps_M_j; + upper(j - 12) += pertub; + lower(j - 12) -= pertub; + // Update Material parameters + material.UpdateBurgersProperties(SolidMath::CalEffectiveStress(sigd_j), 273.); + material.CalResidualBurgers(dt, epsd_i, sigd_j, eps_K_j, eps_K_t, upper, eps_M_t, residual); + up = residual(i); + material.CalResidualBurgers(dt, epsd_i, sigd_j, eps_K_j, eps_K_t, lower, eps_M_t, residual); + low = residual(i); + } + Jac_num(i, j) = (up - low) / (2. * pertub); + ASSERT_NEAR(Jacobian(i, j), Jac_num(i, j), 1.e-8); + } +} + +TEST(SolidProps, YieldMohrCoulomb) +{ + Math_Group::Matrix* data; + data = new Math_Group::Matrix(20); + + // set Constants + for (int i(0); i < 20; i++) + (*data)(i) = (double)i; + + (*data)(7) = 2.; // cohesion + (*data)(9) = 30.; // friction angle + const double phi((*data)(9) * PI / 180.); + (*data)(11) = 29.9; // transition angle + + Minkley::SolidMinkley material(*data); + + // state and trial variables + Eigen::Matrix sig_i; + sig_i.setZero(6, 1); + + // MC formulation + sig_i(0) = -2.; // sigma_1 + sig_i(2) = sig_i(0) * (1. + std::sin(phi)) / (1. - std::sin(phi)) + - 2. * (*data)(7) * std::cos(phi) / (1. - std::sin(phi)); // sigma_3 + sig_i(1) = 0.5 * (sig_i(0) + sig_i(2)); // sigma_2 + + ASSERT_NEAR(0., material.YieldMohrCoulomb(sig_i), 1.e-12); + + // Corner smoothed region + sig_i(0) = -2.; // sigma_1 + sig_i(2) = sig_i(0) * (1. + std::sin(phi)) / (1. - std::sin(phi)) + - 2. * (*data)(7) * std::cos(phi) / (1. - std::sin(phi)); // sigma_3 + sig_i(1) = sig_i(0); // sigma_2 + + // Large tolerance wanted -- yield surfaces deviate in this region. + ASSERT_NEAR(0., material.YieldMohrCoulomb(sig_i), 1.e-2); +} + +TEST(SolidProps, MinkleyJacobianNumeric) +{ + Math_Group::Matrix* data; + data = new Math_Group::Matrix(20); + + // set Constants + (*data)(0) = 63.e3; // Kelvin shear modulus + (*data)(1) = 14.e6; // Kelvin viscosity + (*data)(2) = 12.e3; // Maxwell shear modulus + (*data)(3) = 18.e3; // Maxwell bulk modulus + (*data)(4) = 10.e10; // Maxwell viscosity + (*data)(5) = 4.9; // dependency parameter for " (m) + (*data)(6) = 0.33; // dependency parameter for " (n) + (*data)(7) = 2.; // cohesion + (*data)(8) = 1.; // hardening + (*data)(9) = 30.; // friction angle + (*data)(10) = 10.; // dilatancy angle + (*data)(11) = 28.; // transition angle + (*data)(12) = 0.1; // regularisation + (*data)(13) = 0.; // slope of elesticity temperature dependence + (*data)(14) = 0.; // slope of elesticity temperature dependence + (*data)(15) = 273.; // reference temperature dependency parameter for " + (*data)(16) = 1.; // constant factor for Arrhenius term + (*data)(17) = 0.; // activation energy in Arrhenius term + (*data)(18) = 0.; //second order hardening + (*data)(19) = 0.; //fourth order hardening + Minkley::SolidMinkley material(*data); + + // state and trial variables + Eigen::Matrix eps_i, eps_K_t, eps_K_j, eps_M_t, eps_M_j, eps_pl_j, eps_pl_t; + Eigen::Matrix sig_t, sigd_j, sig_j, epsd_i, epsd_t; + eps_i.setZero(6); + eps_K_t.setZero(6); + eps_K_j.setZero(6); + eps_M_t.setZero(6); + eps_M_j.setZero(6); + eps_pl_t.setZero(6); + eps_pl_j.setZero(6); + sig_t.setZero(6); + sigd_j.setZero(6); + sig_j.setZero(6); + epsd_t.setZero(6); + epsd_i.setZero(6); + + // set a strain increment and perform Kelvin mapping + std::vector strain_bc(6); + strain_bc[0] = 0.01; + strain_bc[1] = 0.02; + strain_bc[2] = 0.03; + strain_bc[3] = 0.04 * 2.; + strain_bc[4] = 0.05 * 2.; + strain_bc[5] = 0.06 * 2.; + eps_i = SolidMath::Voigt_to_Kelvin_Strain(strain_bc); + eps_i *= 0.01; + + epsd_i = SolidMath::P_dev * eps_i; + + // set nontrivial internal variables + eps_K_j = 0.1 * epsd_i; + eps_M_j = 0.15 * epsd_i; + eps_pl_j = 0.1 * epsd_i; + + const double e_i(SolidMath::CalI1(eps_i)); + const double e_pv_j(0.1 * e_i); + + // guess stress increment (dimensionless) + sigd_j = 2.0 * epsd_i; + sig_j = sigd_j + material.KM0 / material.GM0 * e_i * SolidMath::ivec; + + // Update Material parameters + material.UpdateMinkleyProperties(SolidMath::CalEffectiveStress(sigd_j), 0., 273.); + + Eigen::Matrix Jacobian; + Jacobian.setZero(18, 18); + + // Calculate Jacobian for time step + const double dt(1.5); + material.CalViscoelasticJacobian(dt, sig_j, SolidMath::CalEffectiveStress(sigd_j), Jacobian); + + // Residual vector + Eigen::Matrix residual; + + // Numerically calculate Jacobian + Eigen::Matrix Jac_num; + Jac_num.setZero(18, 18); + + Eigen::Matrix upper, lower; + const double pertub(1.e-8); + double up, low; + + for (size_t i = 0; i < 18; i++) + for (size_t j = 0; j < 18; j++) + { + if (j < 6) + { + upper = sig_j; + lower = sig_j; + upper(j) += pertub; + lower(j) -= pertub; + // Update Material parameters + material.UpdateMinkleyProperties(SolidMath::CalEffectiveStress(SolidMath::P_dev * upper), 0., 273.); + material.CalViscoelasticResidual(dt, epsd_i, e_i, e_pv_j, upper, eps_K_j, eps_K_t, eps_M_j, eps_M_t, + eps_pl_j, residual); + up = residual(i); + // Update Material parameters + material.UpdateMinkleyProperties(SolidMath::CalEffectiveStress(SolidMath::P_dev * lower), 0., 273.); + material.CalViscoelasticResidual(dt, epsd_i, e_i, e_pv_j, lower, eps_K_j, eps_K_t, eps_M_j, eps_M_t, + eps_pl_j, residual); + low = residual(i); + } + else if (j < 12) + { + upper = eps_K_j; + lower = eps_K_j; + upper(j - 6) += pertub; + lower(j - 6) -= pertub; + // Update Material parameters + material.UpdateMinkleyProperties(SolidMath::CalEffectiveStress(SolidMath::P_dev * sig_j), 0., 273.); + material.CalViscoelasticResidual(dt, epsd_i, e_i, e_pv_j, sig_j, upper, eps_K_t, eps_M_j, eps_M_t, + eps_pl_j, residual); + up = residual(i); + material.CalViscoelasticResidual(dt, epsd_i, e_i, e_pv_j, sig_j, lower, eps_K_t, eps_M_j, eps_M_t, + eps_pl_j, residual); + low = residual(i); + } + else + { + upper = eps_M_j; + lower = eps_M_j; + upper(j - 12) += pertub; + lower(j - 12) -= pertub; + // Update Material parameters + material.UpdateMinkleyProperties(SolidMath::CalEffectiveStress(SolidMath::P_dev * sig_j), 0., 273.); + material.CalViscoelasticResidual(dt, epsd_i, e_i, e_pv_j, sig_j, eps_K_j, eps_K_t, upper, eps_M_t, + eps_pl_j, residual); + up = residual(i); + material.CalViscoelasticResidual(dt, epsd_i, e_i, e_pv_j, sig_j, eps_K_j, eps_K_t, lower, eps_M_t, + eps_pl_j, residual); + low = residual(i); + } + Jac_num(i, j) = (up - low) / (2. * pertub); + ASSERT_NEAR(Jacobian(i, j), Jac_num(i, j), 1.e-8); + } +} + +TEST(SolidProps, TensorInversion) +{ + // state and trial variables + Eigen::Matrix vec, inv, inv_t; + vec(0) = 1.7; + vec(1) = -1.3; + vec(2) = 0.8; + vec(3) = 1.1 * std::sqrt(2.); + vec(4) = -2.0 * std::sqrt(2.); + vec(5) = 0.3 * std::sqrt(2.); + + // Python result + inv(0) = 0.4693174411; + inv(1) = -0.1182605457; + inv(2) = 0.3184654065; + inv(3) = 0.1378154391 * std::sqrt(2.); + inv(4) = -0.3473321538 * std::sqrt(2.); + inv(5) = 0.1685445572 * std::sqrt(2.); + + // OGS result + inv_t = SolidMath::InvertVector(vec); + + for (size_t i = 0; i < 6; i++) + ASSERT_NEAR(inv_t(i), inv(i), 1.e-10); +} + +TEST(SolidProps, MinkleyFullResidual) +{ + Math_Group::Matrix* data; + data = new Math_Group::Matrix(20); + + // set Constants + (*data)(0) = 63.e3; // Kelvin shear modulus + (*data)(1) = 8.37e4; // Kelvin viscosity + (*data)(2) = 12.e3; // Maxwell shear modulus + (*data)(3) = 18.e3; // Maxwell bulk modulus + (*data)(4) = 4.03e7; // Maxwell viscosity + (*data)(5) = 4.9; // dependency parameter for " (m) + (*data)(6) = 0.33; // dependency parameter for " (n) + (*data)(7) = 2.; // cohesion + (*data)(8) = 100.; // hardening + (*data)(9) = 30.; // friction angle + (*data)(10) = 10.; // dilatancy angle + (*data)(11) = 29.; // transition angle + (*data)(12) = 0.01; // regularisation + (*data)(13) = 0.; // slope of elesticity temperature dependence + (*data)(14) = 0.; // slope of elesticity temperature dependence + (*data)(15) = 273.15; // reference temperature dependency parameter for " + (*data)(16) = 1.; // constant factor for Arrhenius term + (*data)(17) = 0.; // activation energy in Arrhenius term + (*data)(18) = 0.; //second order hardening + (*data)(19) = 0.; //fourth order hardening + Minkley::SolidMinkley material(*data); + + // state and trial variables + Eigen::Matrix eps_i, eps_K_t, eps_K_j, eps_M_t, eps_M_j, eps_pl_t, eps_pl_j; + Eigen::Matrix sig_t, sigd_j, sig_j, epsd_i, epsd_t; + eps_i.setZero(6); + eps_K_t.setZero(6); + eps_K_j.setZero(6); + eps_M_t.setZero(6); + eps_M_j.setZero(6); + eps_pl_t.setZero(6); + eps_pl_j.setZero(6); + sig_t.setZero(6); + sigd_j.setZero(6); + sig_j.setZero(6); + epsd_t.setZero(6); + epsd_i.setZero(6); + + // Residual vector + Eigen::Matrix residual; + + // set a strain increment and perform Kelvin mapping + std::vector strain_bc(6); + strain_bc[0] = 0.01; + strain_bc[1] = 0.02; + strain_bc[2] = 0.03; + strain_bc[3] = 0.04 * 2.; + strain_bc[4] = 0.05 * 2.; + strain_bc[5] = 0.06 * 2.; + eps_i = SolidMath::Voigt_to_Kelvin_Strain(strain_bc); + eps_i *= -0.01; + + epsd_i = SolidMath::P_dev * eps_i; + + const double e_i(SolidMath::CalI1(eps_i)); + + // guess stress increment (dimensionless) + sigd_j = 2.0 * epsd_i; + sig_j = sigd_j + material.KM / material.GM * e_i * SolidMath::ivec; + + // Calculate residual for time step + const double dt(10.); + // Modify plastic arrays to ensure full check + // do get nonzero flow driven residuals + const double lam_j(0.01); + eps_K_j = 0.1 * epsd_i; + eps_M_j = 0.15 * epsd_i; + eps_pl_j = 0.2 * epsd_i; + const double e_v_i = 0.2 * e_i; + const double e_eff_i = 1.2 * e_i; + + // Update Material parameters + material.UpdateMinkleyProperties(SolidMath::CalEffectiveStress(sigd_j), e_eff_i, 273.15); + + material.CalViscoplasticResidual(dt, epsd_i, e_i, sig_j, eps_K_j, eps_K_t, eps_M_j, eps_M_t, eps_pl_j, eps_pl_t, + e_v_i, 0., e_eff_i, 0., lam_j, residual); + + Eigen::Matrix residual_Python; + + // These values are taken from a constituive model test script coded in Python + residual_Python(0) = -9e-05; + residual_Python(1) = -0.00018; + residual_Python(2) = -0.00027; + residual_Python(3) = -0.000509116882454; + residual_Python(4) = -0.000636396103068; + residual_Python(5) = -0.000763675323681; + residual_Python(6) = -5.81003584229e-06; + residual_Python(7) = 8.19126399401e-21; + residual_Python(8) = 5.81003584229e-06; + residual_Python(9) = 3.28665259442e-05; + residual_Python(10) = 4.10831574302e-05; + residual_Python(11) = 4.92997889163e-05; + residual_Python(12) = -0.144920047348; + residual_Python(13) = 1.17842582492e-16; + residual_Python(14) = 0.144920047348; + residual_Python(15) = 0.819791585675; + residual_Python(16) = 1.02473948209; + residual_Python(17) = 1.22968737851; + residual_Python(18) = -0.00140743973352; + residual_Python(19) = 0.00173295588137; + residual_Python(20) = -0.000325516147857; + residual_Python(21) = 0.00127776404536; + residual_Python(22) = 0.00237699541911; + residual_Python(23) = 0.00614979374969; + residual_Python(24) = -0.00174848177667; + residual_Python(25) = -0.00587444554726; + residual_Python(26) = 0.000736664715603; + + for (size_t i = 0; i < 27; i++) + { + ASSERT_NEAR(residual(i), residual_Python(i), 1.e-10); + } +} + +TEST(SolidProps, MinkleyFullJacobian) +{ + Math_Group::Matrix* data; + data = new Math_Group::Matrix(20); + + // set Constants + (*data)(0) = 63.e3; // Kelvin shear modulus + (*data)(1) = 8.37e4; // Kelvin viscosity + (*data)(2) = 12.e3; // Maxwell shear modulus + (*data)(3) = 18.e3; // Maxwell bulk modulus + (*data)(4) = 4.03e7; // Maxwell viscosity + (*data)(5) = 4.9; // dependency parameter for " (m) + (*data)(6) = 0.33; // dependency parameter for " (n) + (*data)(7) = 2.; // cohesion + (*data)(8) = 100.; // hardening + (*data)(9) = 30.; // friction angle + (*data)(10) = 10.; // dilatancy angle + (*data)(11) = 29.; // transition angle + (*data)(12) = 0.01; // regularisation + (*data)(13) = 0.; // slope of elesticity temperature dependence + (*data)(14) = 0.; // slope of elesticity temperature dependence + (*data)(15) = 273.15; // reference temperature dependency parameter for " + (*data)(16) = 1.; // constant factor for Arrhenius term + (*data)(17) = 0.; // activation energy in Arrhenius term + (*data)(18) = -1.e4; //second order hardening + (*data)(19) = -1.e6; //fourth order hardening + Minkley::SolidMinkley material(*data); + + // state and trial variables + Eigen::Matrix eps_i, eps_K_t, eps_K_j, eps_M_t, eps_M_j, eps_pl_t, eps_pl_j; + Eigen::Matrix sig_t, sigd_j, sig_j, epsd_i, epsd_t; + eps_i.setZero(6); + eps_K_t.setZero(6); + eps_K_j.setZero(6); + eps_M_t.setZero(6); + eps_M_j.setZero(6); + eps_pl_t.setZero(6); + eps_pl_j.setZero(6); + sig_t.setZero(6); + sigd_j.setZero(6); + sig_j.setZero(6); + epsd_t.setZero(6); + epsd_i.setZero(6); + + // Residual vector + Eigen::Matrix residual; + + // set a strain increment and perform Kelvin mapping + std::vector strain_bc(6); + strain_bc[0] = 0.01; + strain_bc[1] = 0.02; + strain_bc[2] = 0.03; + strain_bc[3] = 0.04 * 2.; + strain_bc[4] = 0.05 * 2.; + strain_bc[5] = 0.06 * 2.; + eps_i = SolidMath::Voigt_to_Kelvin_Strain(strain_bc); + eps_i *= -0.01; + + epsd_i = SolidMath::P_dev * eps_i; + + const double e_i(SolidMath::CalI1(eps_i)); + + // guess stress increment (dimensionless) + sigd_j = 2.0 * epsd_i; + sig_j = sigd_j + material.KM / material.GM * e_i * SolidMath::ivec; + + // Calculate residual for time step + const double dt(10.); + // Modify plastic arrays to ensure full check + // do get nonzero flow driven residuals + const double lam_j(0.01); + eps_K_j = 0.1 * epsd_i; + eps_M_j = 0.15 * epsd_i; + eps_pl_j = 0.2 * epsd_i; + const double e_v_i = 0.2 * e_i; + const double e_eff_i = 1.2 * e_i; + + // Update Material parameters + material.UpdateMinkleyProperties(SolidMath::CalEffectiveStress(sigd_j), e_eff_i, 273.); + + // Calculate Jacobian for time step + Eigen::Matrix Jacobian; + material.CalViscoplasticJacobian(dt, sig_j, SolidMath::CalEffectiveStress(sigd_j), lam_j, e_eff_i, Jacobian); + + // Numerically calculate Jacobian + Eigen::Matrix Jac_num; + Jac_num.setZero(27, 27); + + Eigen::Matrix upper, lower; + const double pertub(1.e-8); + double up, low; + + for (size_t i = 0; i < 27; i++) + for (size_t j = 0; j < 27; j++) + { + if (j < 6) + { + upper = sig_j; + lower = sig_j; + upper(j) += pertub; + lower(j) -= pertub; + // Update Material parameters + material.UpdateMinkleyProperties(SolidMath::CalEffectiveStress(SolidMath::P_dev * upper), e_eff_i, + 273.); + material.CalViscoplasticResidual(dt, epsd_i, e_i, upper, eps_K_j, eps_K_t, eps_M_j, eps_M_t, eps_pl_j, + eps_pl_t, e_v_i, 0., e_eff_i, 0., lam_j, residual); + up = residual(i); + // Update Material parameters + material.UpdateMinkleyProperties(SolidMath::CalEffectiveStress(SolidMath::P_dev * lower), e_eff_i, + 273.); + material.CalViscoplasticResidual(dt, epsd_i, e_i, lower, eps_K_j, eps_K_t, eps_M_j, eps_M_t, eps_pl_j, + eps_pl_t, e_v_i, 0., e_eff_i, 0., lam_j, residual); + low = residual(i); + } + else if (j < 12) + { + upper = eps_K_j; + lower = eps_K_j; + upper(j - 6) += pertub; + lower(j - 6) -= pertub; + // Update Material parameters + material.UpdateMinkleyProperties(SolidMath::CalEffectiveStress(SolidMath::P_dev * sig_j), e_eff_i, + 273.); + material.CalViscoplasticResidual(dt, epsd_i, e_i, sig_j, upper, eps_K_t, eps_M_j, eps_M_t, eps_pl_j, + eps_pl_t, e_v_i, 0., e_eff_i, 0., lam_j, residual); + up = residual(i); + material.CalViscoplasticResidual(dt, epsd_i, e_i, sig_j, lower, eps_K_t, eps_M_j, eps_M_t, eps_pl_j, + eps_pl_t, e_v_i, 0., e_eff_i, 0., lam_j, residual); + low = residual(i); + } + else if (j < 18) + { + upper = eps_M_j; + lower = eps_M_j; + upper(j - 12) += pertub; + lower(j - 12) -= pertub; + // Update Material parameters + material.UpdateMinkleyProperties(SolidMath::CalEffectiveStress(SolidMath::P_dev * sig_j), e_eff_i, + 273.); + material.CalViscoplasticResidual(dt, epsd_i, e_i, sig_j, eps_K_j, eps_K_t, upper, eps_M_t, eps_pl_j, + eps_pl_t, e_v_i, 0., e_eff_i, 0., lam_j, residual); + up = residual(i); + material.CalViscoplasticResidual(dt, epsd_i, e_i, sig_j, eps_K_j, eps_K_t, lower, eps_M_t, eps_pl_j, + eps_pl_t, e_v_i, 0., e_eff_i, 0., lam_j, residual); + low = residual(i); + } + else if (j < 24) + { + upper = eps_pl_j; + lower = eps_pl_j; + upper(j - 18) += pertub; + lower(j - 18) -= pertub; + // Update Material parameters + material.UpdateMinkleyProperties(SolidMath::CalEffectiveStress(SolidMath::P_dev * sig_j), e_eff_i, + 273.); + material.CalViscoplasticResidual(dt, epsd_i, e_i, sig_j, eps_K_j, eps_K_t, eps_M_j, eps_M_t, upper, + eps_pl_t, e_v_i, 0., e_eff_i, 0., lam_j, residual); + up = residual(i); + material.CalViscoplasticResidual(dt, epsd_i, e_i, sig_j, eps_K_j, eps_K_t, eps_M_j, eps_M_t, lower, + eps_pl_t, e_v_i, 0., e_eff_i, 0., lam_j, residual); + low = residual(i); + } + else if (j < 25) + { + // Update Material parameters + material.UpdateMinkleyProperties(SolidMath::CalEffectiveStress(SolidMath::P_dev * sig_j), e_eff_i, + 273.); + material.CalViscoplasticResidual(dt, epsd_i, e_i, sig_j, eps_K_j, eps_K_t, eps_M_j, eps_M_t, eps_pl_j, + eps_pl_t, e_v_i + pertub, 0., e_eff_i, 0., lam_j, residual); + up = residual(i); + material.CalViscoplasticResidual(dt, epsd_i, e_i, sig_j, eps_K_j, eps_K_t, eps_M_j, eps_M_t, eps_pl_j, + eps_pl_t, e_v_i - pertub, 0., e_eff_i, 0., lam_j, residual); + low = residual(i); + } + else if (j < 26) + { + // Update Material parameters + material.UpdateMinkleyProperties(SolidMath::CalEffectiveStress(SolidMath::P_dev * sig_j), + e_eff_i + pertub, 273.); + material.CalViscoplasticResidual(dt, epsd_i, e_i, sig_j, eps_K_j, eps_K_t, eps_M_j, eps_M_t, eps_pl_j, + eps_pl_t, e_v_i, 0., e_eff_i + pertub, 0., lam_j, residual); + up = residual(i); + // Update Material parameters + material.UpdateMinkleyProperties(SolidMath::CalEffectiveStress(SolidMath::P_dev * sig_j), + e_eff_i - pertub, 273.); + material.CalViscoplasticResidual(dt, epsd_i, e_i, sig_j, eps_K_j, eps_K_t, eps_M_j, eps_M_t, eps_pl_j, + eps_pl_t, e_v_i, 0., e_eff_i - pertub, 0., lam_j, residual); + low = residual(i); + } + else if (j < 27) + { + // Update Material parameters + material.UpdateMinkleyProperties(SolidMath::CalEffectiveStress(SolidMath::P_dev * sig_j), e_eff_i, + 273.); + material.CalViscoplasticResidual(dt, epsd_i, e_i, sig_j, eps_K_j, eps_K_t, eps_M_j, eps_M_t, eps_pl_j, + eps_pl_t, e_v_i, 0., e_eff_i, 0., lam_j + pertub, residual); + up = residual(i); + material.CalViscoplasticResidual(dt, epsd_i, e_i, sig_j, eps_K_j, eps_K_t, eps_M_j, eps_M_t, eps_pl_j, + eps_pl_t, e_v_i, 0., e_eff_i, 0., lam_j - pertub, residual); + low = residual(i); + } + Jac_num(i, j) = (up - low) / (2. * pertub); + // std::cout << i << " " << j << " " << Jac_num(i,j) << " " << Jacobian(i,j) << std::endl; + ASSERT_NEAR(Jacobian(i, j), Jac_num(i, j), 1.e-6); + } +}