Select Git revision
-
Ludovic Noels authoredLudovic Noels authored
j2plast.c 66.52 KiB
#ifndef J2PLAST_C
#define J2PLAST_C 1
//#ifndef max
//#define max(a,b) ( ((a) > (b)) ? (a) : (b) )
//#endif
#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include "j2plast.h"
#include "matrix_operations.h"
#include "elasticlcc.h"
#ifdef NONLOCALGMSH
#include "ID.h"
using namespace MFH;
#endif
#define DP_MIN 1e-60
//ELASTO-PLASTIC CONSTITUTIVE BOX
int constbox_ep(double E, double nu, double sy0, int htype, double hmod1, double hmod2, double hp0,
double hexp, double *dstrn, double *strs_n, double *strs, double* pstrn_n, double* pstrn,
double p_n, double *p, double *twomu_iso, double* dtwomu_iso, double** Calgo, double*** dCalgo, double* dpdE, bool residual)
{
int i,j,k,it,error;
double strseqtr, trstrstr,strseqtrNor0; //Variables for the elastic prediction
double R, dR, ddR, strseq, G; //Variables for the correction phase
double kp;
double Dp;
double h, tmp1, tmp2;
double cp=0.;
double tol = 1.e-12; //tolerance on the residual
double res, res0; //residual
double mu_iso, k_iso;
double G3, h3;
double sq2;
sq2 = sqrt(2.);
R=0.;
dR=0.;
ddR=0.;
G = E/(2.*(1.+nu));
k_iso = E/(3.*(1.-(2.*nu)));
G3 = G*G*G;
static bool initialized=false;
// initialize only once!
static double **Cel=NULL;
static double **NoN=NULL;
static double **Idev=NULL;
static double **mat1=NULL;
static double *strstr=NULL;
static double *Str=NULL; //used to compute the modified normal
static double *Ntr=NULL; //used to compute the modified normal
static double *Stau=NULL;
static double *dstrs=NULL;
static double *strstrNor0=NULL; //used to compute the modified normal
bool stay=true;
if(!initialized)
{
mallocmatrix(&Cel, 6, 6);
mallocmatrix(&NoN,6,6);
mallocmatrix(&Idev,6,6); mallocmatrix(&mat1,6,6);
mallocvector(&strstr,6);
mallocvector(&Str,6);
mallocvector(&Ntr,6);
mallocvector(&Stau,6);
mallocvector(&dstrs,6);
mallocvector(&strstrNor0,6);
initialized = true;
}
cleartens4(Cel);
cleartens4(NoN);
cleartens4(Idev);
cleartens4(mat1);
cleartens2(strstr);
cleartens2(Str);
cleartens2(Ntr);
cleartens2(Stau);
cleartens2(dstrs);
cleartens2(strstrNor0);
//deviatoric operator
for(i=0;i<3;i++){
for(j=0;j<3;j++){
Idev[i][j]=-1./3.;
}
Idev[i][i]+=1.;
Idev[i+3][i+3]=1.;
}
//reset tangent operators to zero
for(i=0;i<6;i++){
for(j=0;j<6;j++){
Calgo[i][j]=0.;
}
}
//elastic predictor
compute_Hooke(E,nu,Cel);
//compute trial stress
contraction42(Cel, dstrn, dstrs);
addtens2(1., strs_n, 1., dstrs, strstr);
strseqtr = vonmises(strstr);
strseqtrNor0 = vonmises(dstrs);
//normal to the yield surface
if(!residual)
{
dev(strstr, Str);
if(strseqtr > 0){
addtens2(1.5/strseqtr, Str, 0, Str, Ntr);
produit_tensoriel(Ntr,Ntr,NoN);
}
}
else
{
dev(dstrs, Str);
if(strseqtrNor0 > 0){
addtens2(1.5/strseqtrNor0, Str, 0, Str, Ntr);
produit_tensoriel(Ntr,Ntr,NoN);
}
}
//elastic increment
Dp=DP_MIN;
j2hard (p_n+Dp, sy0, htype, hmod1, hmod2, hp0, hexp, &R, &dR,&ddR);
strseq= sy0 + R;
h = (3.*G) + dR;
kp = strseq + 3.*G*Dp - strseqtr;
if(kp>=0.){
//printf("elastic: !");
copyvect(strstr, strs, 6); copyvect(pstrn_n, pstrn,6);
*p = p_n;
Dp=0.;
copymatr(Cel,Calgo,6,6);
for(i=0;i<6;i++){
for(j=0;j<6;j++){
for(k=0;k<6;k++){
dCalgo[i][j][k]=0.;
}
}
}
//***** By Wu Ling 06/04/2011************************************************************************************
//compute dpdE accumulate plastice strain p derivative w.r.t. strain increment
cleartens2(dpdE);
//******* By Wu Ling 06/04/2011**********************************************************************************
}
//Plastic correction: return mapping algorithm
else{
Dp = fmax(DP_MIN,-kp/h); //correction must be such that Dp >= DP_MIN
it=0;
res=1.;
res0=res;
while(res>tol && stay){
it++;
if(it>20){
printf("no convergence in return mapping algorithm, p: %lf\n",*p);
stay=false;
Dp=1.;
//return(1);
}
*p=p_n+Dp;
j2hard (*p, sy0, htype, hmod1, hmod2, hp0, hexp, &R, &dR,&ddR);
kp = 3.*G*Dp + R + sy0 - strseqtr;
res = fabs(kp)/G;
if(it==1||res<res0){
res0=res;
cp = - kp/(3.*G+dR);
tmp1 = 0.99*(DP_MIN-Dp);
if(cp<tmp1) {
//printf("new Dp below Dpmin : correction!\n");
cp=tmp1;
}
Dp+=cp;
}
else{
//printf("residual is not decreasing: correction!\n");
cp=cp/2.;
Dp=Dp-cp;
res0=res;
}
//printf("j2plast iteration: %d, res: %.10e, p: %.10e\n",it,res,*p);
}
//printf("Number of iterations in radial return algorithm: %d\n",it);
strseq = sy0 + R;
addtens2(1., strstr,-2.*G*(Dp), Ntr, strs);
h=(3.*G) + dR;
addtens2(1.,pstrn_n, Dp, Ntr, pstrn);
h3 = h*h*h;
tmp1 = 4.*G*G/h;
tmp2 = 4.*G*G*(*p-p_n)/strseqtr;
addtens4(1.,Cel,(-tmp1+tmp2),NoN, Calgo);
addtens4(1.,Calgo,-1.5*tmp2,Idev,Calgo);
addtens4(1.5,Idev,-1.,NoN,mat1);
for(i=0;i<6;i++){
for(j=0;j<6;j++){
for(k=0;k<6;k++){
dCalgo[i][j][k] = (8.*G3*ddR/h3)*NoN[i][j]*Ntr[k] - (16.*G3/(h*strseqtr))*Ntr[i]*mat1[j][k] - mat1[i][j]*8.*G3*Ntr[k]*((1./(h*strseqtr))-(Dp/(strseqtr*strseqtr))) + (16.*G3*Dp/(strseqtr*strseqtr))*Ntr[i]*mat1[j][k];
}
}
}
//***** By Wu Ling 06/04/2011************************************************************************************
//compute dpdE accumulate plastice strain p derivative w.r.t. strain increment
addtens2((2.*G/h), Ntr, 0, Ntr, dpdE);
//******* By Wu Ling 06/04/2011**********************************************************************************
} //end case plastic correction
//compute isoJ2 tangent operator and its derivative w.r.t. strain increment
if(Dp<DP_MIN){ //that is, Dp=0
mu_iso=G;
for(i=0;i<6;i++){
dtwomu_iso[i]=0.;
}
}
else {
mu_iso = G*(1.-(3.*G/h));
for(i=0;i<6;i++){
dtwomu_iso[i] = 12.*pow(G/h,3)*ddR*Ntr[i];
}
}
*twomu_iso = 2.*mu_iso;
// no free as initialized once
// freematrix(Cel, 6);
// freematrix(Idev,6);
// freematrix(NoN,6);
// freematrix(mat1,6);
if(!stay) return(1); //no convergence
return(0);
}
int constbox_ep_res(double E, double nu, double sy0, int htype, double hmod1, double hmod2, double hp0,
double hexp, double *dstrn, double *strs_n, double *strs, double* pstrn_n, double* pstrn,
double p_n, double *p, double *twomu_iso, double* dtwomu_iso, double** Calgo, double*** dCalgo, double* dpdE, bool residual)
{
int i,j,k,it,error;
double strseqtr, strseqtrNor0, trstrstr; //Variables for the elastic prediction
double R, dR, ddR, strseq, G; //Variables for the correction phase
double kp;
double Dp;
double h, tmp1, tmp2;
double cp=0.;
double tol = 1.e-12; //tolerance on the residual
double res, res0; //residual
double mu_iso, k_iso; double G3, h3;
double sq2;
sq2 = sqrt(2.);
R=0.;
dR=0.;
ddR=0.;
G = E/(2.*(1.+nu));
k_iso = E/(3.*(1.-(2.*nu)));
G3 = G*G*G;
static bool initialized=false;
// initialize only once!
static double **Cel=NULL;
static double **NoN=NULL;
static double **Idev=NULL;
static double **mat1=NULL;
static double *strstr=NULL;
static double *strstrNor0=NULL; //used to compute the modified normal
static double *Str=NULL; //used to compute the modified normal
static double *Ntr=NULL; //used to compute the modified normal
static double *Stau=NULL;
static double *dstrs=NULL;
bool stay=true;
if(!initialized)
{
mallocmatrix(&Cel, 6, 6);
mallocmatrix(&NoN,6,6);
mallocmatrix(&Idev,6,6);
mallocmatrix(&mat1,6,6);
mallocvector(&strstr,6);
mallocvector(&strstrNor0,6);
mallocvector(&Str,6);
mallocvector(&Ntr,6);
mallocvector(&Stau,6);
mallocvector(&dstrs,6);
initialized = true;
}
cleartens4(Cel);
cleartens4(NoN);
cleartens4(Idev);
cleartens4(mat1);
cleartens2(strstr);
cleartens2(strstrNor0);
cleartens2(Str);
cleartens2(Ntr);
cleartens2(Stau);
cleartens2(dstrs);
//deviatoric operator
for(i=0;i<3;i++){
for(j=0;j<3;j++){
Idev[i][j]=-1./3.;
}
Idev[i][i]+=1.;
Idev[i+3][i+3]=1.;
}
//reset tangent operators to zero
for(i=0;i<6;i++){
for(j=0;j<6;j++){
Calgo[i][j]=0.;
}
}
//elastic predictor
compute_Hooke(E,nu,Cel);
//compute trial stress contraction42(Cel, dstrn, dstrs);
addtens2(1.0, strs_n, 1.0, dstrs, strstr);
strseqtrNor0 = vonmises(dstrs);
strseqtr = vonmises(strstr);
//normal to the yield surface
dev(dstrs, Str);
if(strseqtrNor0 > 0){
addtens2(1.5/strseqtrNor0, Str, 0, Str, Ntr);
produit_tensoriel(Ntr,Ntr,NoN);
}
//elastic increment
Dp=DP_MIN;
j2hard (p_n+Dp, sy0, htype, hmod1, hmod2, hp0, hexp, &R, &dR,&ddR);
strseq= sy0 + R;
h = (3.*G) + dR;
kp = strseq + 3.*G*Dp - strseqtr;
if(kp>=0.){
//printf("elastic: !");
copyvect(strstr, strs, 6);
copyvect(pstrn_n, pstrn,6);
*p = p_n;
Dp=0.;
copymatr(Cel,Calgo,6,6);
for(i=0;i<6;i++){
for(j=0;j<6;j++){
for(k=0;k<6;k++){
dCalgo[i][j][k]=0.;
}
}
}
//***** By Wu Ling 06/04/2011************************************************************************************
//compute dpdE accumulate plastice strain p derivative w.r.t. strain increment
cleartens2(dpdE);
//******* By Wu Ling 06/04/2011**********************************************************************************
}
//Plastic correction: return mapping algorithm
else{
Dp = fmax(DP_MIN,-kp/h); //correction must be such that Dp >= DP_MIN
it=0;
res=1.;
res0=res;
while(res>tol && stay){
it++;
if(it>20){
printf("no convergence in return mapping algorithm, p: %lf\n",*p);
//return(1);
stay=false;
Dp=1.;
}
*p=p_n+Dp;
j2hard (*p, sy0, htype, hmod1, hmod2, hp0, hexp, &R, &dR,&ddR);
kp = 3.*G*Dp + R + sy0 - strseqtr;
res = fabs(kp)/G;
if(it==1||res<res0){
res0=res;
cp = - kp/(3.*G+dR);
tmp1 = 0.99*(DP_MIN-Dp);
if(cp<tmp1) {
//printf("new Dp below Dpmin : correction!\n");
cp=tmp1;
}
Dp+=cp;
}
else{
//printf("residual is not decreasing: correction!\n");
cp=cp/2.;
Dp=Dp-cp;
res0=res;
}
//printf("j2plast iteration: %d, res: %.10e, p: %.10e\n",it,res,*p);
}
//printf("Number of iterations in radial return algorithm: %d\n",it);
strseq = sy0 + R;
addtens2(1., strstr,-2.*G*(Dp), Ntr, strs);
h=(3.*G) + dR;
addtens2(1.,pstrn_n, Dp, Ntr, pstrn);
h3 = h*h*h;
tmp1 = 4.*G*G/h;
tmp2 = 4.*G*G*(*p-p_n)/strseqtr;
addtens4(1.,Cel,(-tmp1+tmp2),NoN, Calgo);
addtens4(1.,Calgo,-1.5*tmp2,Idev,Calgo);
addtens4(1.5,Idev,-1.,NoN,mat1);
for(i=0;i<6;i++){
for(j=0;j<6;j++){
for(k=0;k<6;k++){
dCalgo[i][j][k] = (8.*G3*ddR/h3)*NoN[i][j]*Ntr[k] - (16.*G3/(h*strseqtr))*Ntr[i]*mat1[j][k] - mat1[i][j]*8.*G3*Ntr[k]*((1./(h*strseqtr))-(Dp/(strseqtr*strseqtr))) + (16.*G3*Dp/(strseqtr*strseqtr))*Ntr[i]*mat1[j][k];
}
}
}
//***** By Wu Ling 06/04/2011************************************************************************************
//compute dpdE accumulate plastice strain p derivative w.r.t. strain increment
addtens2((2.*G/h), Ntr, 0, Ntr, dpdE);
//******* By Wu Ling 06/04/2011**********************************************************************************
} //end case plastic correction
//compute isoJ2 tangent operator and its derivative w.r.t. strain increment
if(Dp<DP_MIN){ //that is, Dp=0
mu_iso=G;
for(i=0;i<6;i++){
dtwomu_iso[i]=0.;
}
}
else {
mu_iso = G*(1.-(3.*G/h));
for(i=0;i<6;i++){
dtwomu_iso[i] = 12.*pow(G/h,3)*ddR*Ntr[i];
}
}
*twomu_iso = 2.*mu_iso;
// no free as initialized once
// freematrix(Cel, 6);
// freematrix(Idev,6);
// freematrix(NoN,6);
// freematrix(mat1,6);
if(!stay) return(1); //no convergence
return(0);
}
//PRESSURE SENSITIVE ELASTO-PLASTIC CONSTITUTIVE BOX
int constbox_ep(double E, double nu, double sy0, int htype, double hmod1, double hmod2, double hp0, double alpha_DP, double m_DP, double nup, double hexp, double *dstrn, double *strs_n, double *strs, double* pstrn_n, double* pstrn, double p_n, double *p, double *twomu_iso, double* dtwomu_iso, double *twok_iso, double* dtwok_iso, double** Calgo, double*** dCalgo, double* dpdE, bool residual)
{
int i,j,k,it,error;
double strseqtr, trstrstr,strseqtrNor0; //Variables for the elastic prediction
double R, dR, ddR, strseq, G, kappa_el, dkp_dDp, dDp_dGamma, dkp_dGamma, Deriv_kp; //Variables for the correction phase
double kp;
double beta;
double Dp, DGAMMA_MIN, Gamma, phitr, phires;
double h, tmp1, tmp2, ka, NR1_prime_eq, NR2_prime, NR1_eq, NR2;
double cGamma=0.;
double tol = 1.e-8; //tolerance on the residual
double res, res0; //residual
double mu_iso, k_iso;
double G3, h3;
double sq2;
sq2 = sqrt(2.);
R=0.;
dR=0.;
ddR=0.;
G = E/(2.*(1.+nu));
kappa_el = E/(3.*(1.-(2.*nu)));
ka = 1./sqrt(1.+2.*pow(nup,2));
G3 = G*G*G;
beta = (9.-18.*nup)/(2.*nup+2.);
static bool initialized=false;
// initialize only once
static double **Cel=NULL;
static double **tmp_tens4=NULL;
static double **tmp_tens42=NULL;
static double **NoN=NULL;
static double **Idev=NULL;
static double **Ivol=NULL;
static double **mat1=NULL;
static double *Identity=NULL;
static double *tmp_tens2=NULL;
static double *tmp_tens22=NULL;
static double *dkp_dE=NULL;
static double *dGamma_dEr=NULL;
static double *NR1_prime=NULL;
static double *NR1=NULL;
static double *NR1_dev=NULL;
static double *NR1_prime_dev=NULL;
static double *strstr=NULL;
static double *StrVol=NULL;
static double *Str=NULL;
static double *StrVoltmp=NULL;
static double *Strtmp=NULL;
static double *Ntr=NULL;
static double *Stau=NULL;
static double *dstrs=NULL;
static double *strstrNor0=NULL;
bool stay=true; if(!initialized)
{
mallocmatrix(&Cel, 6, 6);
mallocmatrix(&tmp_tens4, 6, 6);
mallocmatrix(&tmp_tens42, 6, 6);
mallocmatrix(&NoN,6,6);
mallocmatrix(&Idev,6,6);
mallocmatrix(&Ivol,6,6);
mallocmatrix(&mat1,6,6);
mallocvector(&Identity,6);
mallocvector(&tmp_tens2,6);
mallocvector(&tmp_tens22,6);
mallocvector(&dkp_dE,6);
mallocvector(&dGamma_dEr,6);
mallocvector(&NR1_prime,6);
mallocvector(&NR1,6);
mallocvector(&NR1_dev,6);
mallocvector(&NR1_prime_dev,6);
mallocvector(&strstr,6);
mallocvector(&StrVol,6);
mallocvector(&Str,6);
mallocvector(&StrVoltmp,6);
mallocvector(&Strtmp,6);
mallocvector(&Ntr,6);
mallocvector(&Stau,6);
mallocvector(&dstrs,6);
mallocvector(&strstrNor0,6);
initialized = true;
}
cleartens4(Cel);
cleartens4(tmp_tens4);
cleartens4(tmp_tens42);
cleartens4(NoN);
cleartens4(Idev);
cleartens4(Ivol);
cleartens4(mat1);
cleartens2(Identity);
cleartens2(tmp_tens2);
cleartens2(tmp_tens22);
cleartens2(dkp_dE);
cleartens2(dGamma_dEr);
cleartens2(NR1_prime);
cleartens2(NR1);
cleartens2(NR1_dev);
cleartens2(NR1_prime_dev);
cleartens2(strstr);
cleartens2(StrVol);
cleartens2(Str);
cleartens2(Ntr);
cleartens2(Stau);
cleartens2(dstrs);
cleartens2(strstrNor0);
//deviatoric operator and Identity
for(i=0;i<3;i++){
for(j=0;j<3;j++){
Idev[i][j]=-1./3.;
Ivol[i][j]=1./3.;
}
Idev[i][i]+=1.;
Idev[i+3][i+3]=1.;
Identity[i]=1.;
}
//reset tangent operators to zero
for(i=0;i<6;i++){
for(j=0;j<6;j++){
Calgo[i][j]=0.;
}
}
//elastic predictor compute_Hooke(E,nu,Cel);
//compute trial stress
contraction42(Cel, dstrn, dstrs);
addtens2(1., strs_n, 1., dstrs, strstr);
strseqtr = vonmises(strstr);
strseqtrNor0 = vonmises(dstrs);
//normal to the yield surface
if(!residual)
{
dev(strstr, Str);
addtens2(1./3.*trace(strstr),Identity,0.,Identity,StrVol);
if(strseqtr > 0){
addtens2(3., Str, 2.*beta/3., StrVol, Ntr);
produit_tensoriel(Ntr,Ntr,NoN);
}
}
else
{
dev(dstrs, Str);
addtens2(1./3.*trace(dstrs),Identity,0.,Identity,StrVol);
if(strseqtrNor0 > 0){
addtens2(3., Str, 2.*beta/3., StrVol, Ntr);
produit_tensoriel(Ntr,Ntr,NoN);
}
}
//elastic increment
Dp=DP_MIN;
DGAMMA_MIN = DP_MIN/100.;
j2hard (p_n+Dp, sy0, htype, hmod1, hmod2, hp0, hexp, &R, &dR,&ddR);
strseq= sy0 + R;
phitr = 1./3.*trace(strstr);
kp = pow(strseqtr,alpha_DP)/pow(strseq,alpha_DP)-3.*(pow(m_DP,alpha_DP)-1.)/((m_DP+1.)*strseq)*phitr - (pow(m_DP,alpha_DP)+m_DP)/(m_DP+1.);
if(kp<=0.){
//printf("elastic: !");
copyvect(strstr, strs, 6);
copyvect(pstrn_n, pstrn,6);
*p = p_n;
Dp=0.;
copymatr(Cel,Calgo,6,6);
for(i=0;i<6;i++){
for(j=0;j<6;j++){
for(k=0;k<6;k++){
dCalgo[i][j][k]=0.;
}
}
}
//***** By Juan Manuel Calleja 12/2020************************************************************************************
//compute dpdE accumulate plastic strain p derivative w.r.t. strain increment
cleartens2(dpdE);
//******* By Juan Manuel Calleja 12/2020**********************************************************************************
}
//Plastic correction: return mapping algorithm
else{
Gamma = DGAMMA_MIN;
it=0;
res=1.;
res0=res; *p=p_n;
while(res>tol && stay){
if(it>20)
printf("Return mapping algorithm, p: %.10e, Gamma: %.10e, res: %.10e\n",*p, Gamma, res);
it++;
if(it>200){
printf("no convergence in return mapping algorithm, p: %lf\n",*p);
stay=false;
//Gamma +=0.5 ;
//return(1);
}
if(!residual)
{
phires = 0.;
addtens2(1./(1.+6.*G*(Gamma)),strstr,0.,strs_n,NR1_prime);
//dev(NR1_prime, NR1_prime_dev);
dev(strstr, NR1_prime_dev);
addtens2(1./(1.+6.*G*(Gamma)),NR1_prime_dev,0.,strs_n,NR1_prime_dev);
NR1_prime_eq = vonmises(NR1_prime);
NR2_prime = phitr/(1.+2.*kappa_el*(Gamma)*beta);
}
else
{
phires = 1./3.*trace(strs_n);
addtens2(1./(1.+6.*G*(Gamma)),dstrs,0.,strs_n,NR1_prime);
//dev(NR1_prime, NR1_prime_dev);
dev(dstrs, NR1_prime_dev);
addtens2(1./(1.+6.*G*(Gamma)),NR1_prime_dev,0.,strs_n,NR1_prime_dev);
NR1_prime_eq = vonmises(NR1_prime);
NR2_prime = (phitr-phires)/(1.+2.*kappa_el*(Gamma)*beta);
}
Dp = ka*(Gamma)*sqrt(6.*pow(NR1_prime_eq,2.)+4.*pow(beta,2.)/3.*pow(NR2_prime,2.));
*p = p_n+Dp;
j2hard (*p, sy0, htype, hmod1, hmod2, hp0, hexp, &R, &dR,&ddR);
strseq = sy0 + R;
if(!residual)
{
addtens2(1./(1.+6.*G*(Gamma)),strstr,0.,strstr,NR1);
NR1_eq = vonmises(NR1);
NR2 = phitr/(1.+2.*kappa_el*(Gamma)*beta);
}
else
{
addtens2(1./(1.+6.*G*(Gamma)),dstrs,1.,strs_n,NR1);
NR1_eq = vonmises(NR1);
NR2 = (phitr-phires)/(1.+2.*kappa_el*(Gamma)*beta)+phires;
}
tmp1 = pow(NR1_eq,alpha_DP)/pow(strseq,alpha_DP)-3.*NR2*(pow(m_DP,alpha_DP)-1.)/((m_DP+1.)*strseq);
kp = tmp1 - (pow(m_DP,alpha_DP)+m_DP)/(m_DP+1.);
res = fabs(kp);
if(it==1 || res<res0){
res0=res;
dkp_dDp = -pow(NR1_eq,alpha_DP)*alpha_DP/pow(strseq,alpha_DP+1.)*dR+3.*(pow(m_DP,alpha_DP)-1.)/((m_DP+1.)*pow(strseq,2.))*dR*NR2;
if(!residual)
{
tmp1 = 72.*G*pow(strseqtr,2.)/pow(1.+6.*G*(Gamma),3.)+16.*pow(beta,3.)*kappa_el/3.*pow(phitr,2.)/pow(1.+2.*kappa_el*(Gamma)*beta,3.);
} else
{
tmp1 = 72.*G*pow(strseqtrNor0,2.)/pow(1.+6.*G*(Gamma),3.)+16.*pow(beta,3.)*kappa_el/3.*pow(phitr-phires,2.)/pow(1.+2.*kappa_el*(Gamma)*beta,3.);
}
dDp_dGamma = ka*sqrt(6.*pow(NR1_prime_eq,2.)+4.*pow(beta,2.)/3.*pow(NR2_prime,2.))-ka*(Gamma)/(2.*sqrt(6.*pow(NR1_prime_eq,2.)+4.*pow(beta,2.)/3.*pow(NR2_prime,2.)))*tmp1;
dev(NR1,NR1_dev);
dkp_dGamma = -alpha_DP*pow(NR1_eq/strseq,alpha_DP-1.)*(9.*G*contraction22(Str,NR1_dev))/(strseq*pow(1.+6.*G*(Gamma),2.)*NR1_eq)+3.*(pow(m_DP,alpha_DP)-1.)/((m_DP+1.)*strseq)*(1./3.*trace(StrVol)/pow(1.+2.*kappa_el*(Gamma)*beta,2.)*2.*kappa_el*beta);
Deriv_kp = dkp_dDp*dDp_dGamma+dkp_dGamma;
cGamma = - kp/Deriv_kp; // -kp / derivative of kp w.r.t. dp. Newton-Raphson correction.
tmp1 = 0.99*(DGAMMA_MIN-Gamma);
if(cGamma<tmp1) {
printf("new Dp below DGAMMA_MIN : correction!\n");
cGamma=tmp1;
}
Gamma+=cGamma;
}
else{
res0=res;
printf("residual is not decreasing: correction!\n");
cGamma=0.5*cGamma;
Gamma-=cGamma;
}
//printf("j2plast iteration: %d, res: %.10e, p: %.10e\n",it,res,*p);
}
//printf("Number of iterations in radial return algorithm: %d\n",it);
strseq = sy0 + R;
addtens2(1./(1.+6.*G*Gamma),Str, 0., Str, Strtmp);
addtens2(1./(1.+(2.*kappa_el*Gamma*beta)),StrVol, 0., StrVol, StrVoltmp);
addtens2(1.,Strtmp, 1., StrVoltmp, strs);
if(residual)
{
addtens2(1.,strs, 1., strs_n, strs);
}
//addtens2(1.,pstrn_n, 3.*Gamma, NR1_prime, pstrn);
addtens2(1.,pstrn_n, 3.*Gamma, NR1_prime_dev, pstrn);
addtens2(1.,pstrn, 2.*beta/3.*Gamma*NR2_prime, Identity, pstrn);
//Computation of Calgo
contraction44(Idev,Cel,tmp_tens4); //ORDER 4
addtens4(1./(1.+6.*G*(Gamma)),tmp_tens4,0.,tmp_tens4,tmp_tens4); //ORDER 4
contraction42(tmp_tens4,NR1_dev,tmp_tens2); //ORDER 2
addtens2(3./2.*alpha_DP * (pow(NR1_eq,alpha_DP-1.)/NR1_eq),tmp_tens2,0.,tmp_tens2,tmp_tens2); //ORDER 2 58.1
addtens2(kappa_el/(1.+2.*kappa_el*(Gamma)*beta),Identity,0.,Identity,tmp_tens22); //ORDER 2 58.2
addtens2(1/pow(strseq,alpha_DP),tmp_tens2,-3.*(pow(m_DP,alpha_DP)-1.)/((m_DP+1.)*strseq),tmp_tens22, dkp_dE); //ORDER 2 57
addtens2(-1./Deriv_kp,dkp_dE,0.,dkp_dE,dGamma_dEr); //ORDER 2 56
//*************************************************************
addtens2(1.+6.*G*(Gamma),dGamma_dEr,-6.*G*(Gamma),dGamma_dEr,tmp_tens2); //ORDER 2
addtens2(1./pow(1.+6.*G*(Gamma),2.),tmp_tens2,0.,tmp_tens2,tmp_tens2); //ORDER 2
produit_tensoriel(Str,tmp_tens2,tmp_tens4); //ORDER 4
addtens2(1.+2.*kappa_el*beta*(Gamma),dGamma_dEr,-2.*kappa_el*beta*(Gamma),dGamma_dEr,tmp_tens2); //ORDER 2
addtens2(1./pow(1.+2.*kappa_el*beta*(Gamma),2.),tmp_tens2,0.,tmp_tens2,tmp_tens2); //ORDER 2
produit_tensoriel(StrVol,tmp_tens2,tmp_tens42); //ORDER 4
addtens4(2.*G*(Gamma)/(1.+6.*G*(Gamma)),Idev,1.,tmp_tens4,tmp_tens4); //ORDER 4
addtens4(3.*kappa_el*(Gamma)/(1.+2.*beta*kappa_el*(Gamma)),Ivol,1.,tmp_tens42,tmp_tens42); //ORDER 4
//************************************************************
addtens4(1.,Cel,-6.*G,tmp_tens4,Calgo);
addtens4(1.,Calgo,-2.*kappa_el*beta,tmp_tens42,Calgo);
//*****************************************************************************************
//compute dpdE accumulate plastice strain p derivative w.r.t. strain increment
addtens2(dDp_dGamma,dGamma_dEr, 0., dGamma_dEr, dpdE); //ORDER 2
//*****************************************************************************************
} //end case plastic correction
//compute isoJ2 tangent operator and its derivative w.r.t. strain increment
if(Dp<DP_MIN){ //that is, Dp=0
mu_iso=G;
k_iso=kappa_el;
for(i=0;i<6;i++){
dtwomu_iso[i]=0.;
dtwok_iso[i]=0.;
}
}
else {
mu_iso = G*(1-(6.*G*Gamma)/(1.+6.*G*Gamma));
k_iso= kappa_el*(1.-(2.*beta*kappa_el*Gamma)/(1.+2.*beta*kappa_el*Gamma));
addtens2(-6.*pow(G,2.)/pow(1.+6.*G*Gamma,2.),dGamma_dEr,0.,dGamma_dEr,dtwok_iso);
addtens2(-2.*beta*pow(kappa_el,2.)/pow(1.+2.*kappa_el*Gamma*beta,2.),dGamma_dEr,0.,dGamma_dEr,dtwok_iso);
}
*twomu_iso = 2.*mu_iso;
*twok_iso = 2.*k_iso;
// no free as initialized once
// freematrix(Cel, 6);
// freematrix(Idev,6);
// freematrix(NoN,6);
// freematrix(mat1,6);
if(!stay) return(1); //no convergence
return(0);
}
//****************************************************
/**********J2 plasticity for large deformation***/
//******************************************************
int constbox_ep(double E, double nu, double sy0, int htype, double hmod1, double hmod2, double hp0, double hexp, double *strstr, double *Rstrs_n, double *strs, double* pstrn_n, double* pstrn, double p_n, double *p, double** Calgo, double*** dCalgo, double* dpdE, bool residual)
{
int i,j,k,it,error;
double strseqtr,strseqtr0, trstrstr; //Variables for the elastic prediction
double R, dR, ddR, strseq, G; //Variables for the correction phase
double kp;
double Dp;
double h, tmp1, tmp2;
double cp=0.;
double tol = 1.e-12; //tolerance on the residual
double res, res0; //residual
double mu_iso, k_iso;
double G3, h3;
double sq2;
sq2 = sqrt(2.);
R=0.;
dR=0.; ddR=0.;
G = E/(2.*(1.+nu));
k_iso = E/(3.*(1.-(2.*nu)));
G3 = G*G*G;
static bool initialized=false;
// initialize only once!
static double **Cel;
static double **NoN;
static double **Idev;
static double **mat1;
static double *strstrNor0=NULL; //used to compute the modified normal
static double *Str=NULL; //used to compute the modified normal
static double *Ntr=NULL; //used to compute the modified normal
static double *Stau=NULL;
static double *dstrs=NULL;
bool stay=true;
if(!initialized)
{
mallocmatrix(&Cel, 6, 6);
mallocmatrix(&NoN,6,6);
mallocmatrix(&Idev,6,6);
mallocmatrix(&mat1,6,6);
mallocvector(&strstrNor0,6);
mallocvector(&Str,6);
mallocvector(&Ntr,6);
mallocvector(&Stau,6);
mallocvector(&dstrs,6);
initialized = true;
}
cleartens4(Cel);
cleartens4(NoN);
cleartens4(Idev);
cleartens4(mat1);
cleartens2(strstrNor0);
cleartens2(Str);
cleartens2(Ntr);
cleartens2(Stau);
cleartens2(dstrs);
//deviatoric operator
for(i=0;i<3;i++){
for(j=0;j<3;j++){
Idev[i][j]=-1./3.;
}
Idev[i][i]+=1.;
Idev[i+3][i+3]=1.;
}
//reset tangent operators to zero
for(i=0;i<6;i++){
for(j=0;j<6;j++){
Calgo[i][j]=0.;
}
}
strseqtr = vonmises(strstr);
//normal to the yield surface
//dev(strstr, Str);
//normal to the residual stress
addtens2(1.0, strstr, -1.0, Rstrs_n, strstrNor0);
dev(strstrNor0, Str);
strseqtr0 = vonmises(strstrNor0);
if(strseqtr > 0.){
addtens2(1.5/strseqtr0, Str, 0., Str, Ntr);
produit_tensoriel(Ntr,Ntr,NoN); }
//elastic increment
Dp=DP_MIN;
j2hard (p_n+Dp, sy0, htype, hmod1, hmod2, hp0, hexp, &R, &dR,&ddR);
strseq= sy0 + R;
h = (3.*G) + dR;
kp = strseq + 3.*G*Dp - strseqtr;
if(kp>=0.){
//printf("elastic: !");
copyvect(strstr, strs, 6);
copyvect(pstrn_n, pstrn,6);
*p = p_n;
Dp=0.;
copymatr(Cel,Calgo,6,6);
for(i=0;i<6;i++){
for(j=0;j<6;j++){
for(k=0;k<6;k++){
dCalgo[i][j][k]=0.;
}
}
}
//***** By Wu Ling 06/04/2011************************************************************************************
//compute dpdE accumulate plastice strain p derivative w.r.t. strain increment
cleartens2(dpdE);
//******* By Wu Ling 06/04/2011**********************************************************************************
}
//Plastic correction: return mapping algorithm
else{
Dp = fmax(DP_MIN,-kp/h); //correction must be such that Dp >= DP_MIN
it=0;
res=1.;
res0=res;
while(res>tol && stay){
it++;
if(it>20){
printf("no convergence in return mapping algorithm, p: %lf\n",*p);
stay=false;
Dp=1.;
//return(1);
}
*p=p_n+Dp;
j2hard (*p, sy0, htype, hmod1, hmod2, hp0, hexp, &R, &dR,&ddR);
kp = 3.*G*Dp + R + sy0 - strseqtr;
res = fabs(kp)/G;
if(it==1||res<res0){
res0=res;
cp = - kp/(3.*G+dR);
tmp1 = 0.99*(DP_MIN-Dp);
if(cp<tmp1) {
//printf("new Dp below Dpmin : correction!\n");
cp=tmp1;
}
Dp+=cp;
}
else{
//printf("residual is not decreasing: correction!\n");
cp=cp/2.; Dp=Dp-cp;
res0=res;
}
//printf("j2plast iteration: %d, res: %.10e, p: %.10e\n",it,res,*p);
}
//printf("Number of iterations in radial return algorithm: %d\n",it);
strseq = sy0 + R;
addtens2(1., strstr,-2.*G*(Dp), Ntr, strs);
h=(3.*G) + dR;
addtens2(1.,pstrn_n, Dp, Ntr, pstrn);
h3 = h*h*h;
tmp1 = 4.*G*G/h;
tmp2 = 4.*G*G*(*p-p_n)/strseqtr0;
addtens4(1.,Cel,(-tmp1+tmp2),NoN, Calgo);
addtens4(1.,Calgo,-1.5*tmp2,Idev,Calgo);
addtens4(1.5,Idev,-1.,NoN,mat1);
for(i=0;i<6;i++){
for(j=0;j<6;j++){
for(k=0;k<6;k++){
dCalgo[i][j][k] = (8.*G3*ddR/h3)*NoN[i][j]*Ntr[k] - (16.*G3/(h*strseqtr0))*Ntr[i]*mat1[j][k] - mat1[i][j]*8.*G3*Ntr[k]*((1./(h*strseqtr0))-(Dp/(strseqtr*strseqtr0))) + (16.*G3*Dp/(strseqtr0*strseqtr0))*Ntr[i]*mat1[j][k];
}
}
}
//***** By Wu Ling 06/04/2011************************************************************************************
//compute dpdE accumulate plastice strain p derivative w.r.t. strain increment
addtens2((2.*G/h), Ntr, 0., Ntr, dpdE);
} //end case plastic correction
if(!stay) return(1); //no convergence
return(0);
}
//end elastoPlastic for large deformation
// EVALUATION OF THE HARDENING FUNCTION
void j2hard (double p, double sy0, int htype, double hmod1, double hmod2, double hp0, double hexp, double* R, double* dR, double* ddR){
double tmp;
//Perfect plasticity
if(hmod1 == 0){
*R = 0.;
*dR = 0.;
*ddR = 0.;
return;
}
switch(htype) {
//POWER LAW HARDENING
case H_POW:
if(p<DP_MIN){
//printf("Should not be here: p: %.10e\n",p);
*R=0.;
if(hexp<1.0){ *dR = 1e20;
*ddR = -1e20;
}
else{
*dR = hmod1;
*ddR = 0.;
}
}
else{
*R=hmod1*pow(p,hexp);
*dR = (*R)*hexp/p;
*ddR = (*dR)*(hexp-1.)/p;
}
break;
//POWER LAW HARDENING FOR DRUCKER PRAGER
case H_POW_DP:
if(p<DP_MIN){
//printf("Should not be here: p: %.10e\n",p);
*R=0.;
if(hexp<1.0){
*dR = 1e20;
*ddR = -1e20;
}
else{
*dR = hmod1;
*ddR = 0.;
}
}
else{
*R=hmod1*pow(p,hexp);
*dR = (*R)*hexp/p;
*ddR = (*dR)*(hexp-1.)/p;
}
break;
//EXPONENTIAL HARDENING
case H_EXP:
if(p<DP_MIN){
//printf("Should not be here: p: %.10e\n",p);
*R=0.;
*dR = hmod1*hexp;
*ddR = -hexp*(*dR);
}
else{
tmp = exp(-hexp*p);
*R = hmod1*(1.-tmp);
*dR = hmod1*hexp*tmp;
*ddR = -hexp*(*dR);
}
break;
//EXPONENTIAL HARDENING FOR DRUCKER PRAGER
case H_EXP_DP:
if(p<DP_MIN){
//printf("Should not be here: p: %.10e\n",p);
*R=0.;
*dR = hmod1*hexp;
*ddR = -hexp*(*dR);
}
else{
tmp = exp(-hexp*p);
*R = hmod1*(1.-tmp);
*dR = hmod1*hexp*tmp;
*ddR = -hexp*(*dR); }
break;
//SWIFT LAW
case H_SWIFT:
if(p<DP_MIN){
//printf("Should not be here: p: %.10e\n",p);
*R=0.;
*dR = sy0*hmod1*hexp;
*ddR = (*dR)*(hexp-1.)*hmod1;
}
else{
tmp = 1.+hmod1*p;
*R = sy0*pow(tmp,hexp) - sy0;
tmp = hmod1/tmp;
*dR = (*R)*hexp*tmp;
*ddR = (*dR)*(hexp-1.)*tmp;
}
break;
//SWIFT LAW FOR DRUCKER PRAGER
case H_SWIFT_DP:
if(p<DP_MIN){
//printf("Should not be here: p: %.10e\n",p);
*R=0.;
*dR = sy0*hmod1*hexp;
*ddR = (*dR)*(hexp-1.)*hmod1;
}
else{
tmp = 1.+hmod1*p;
*R = sy0*pow(tmp,hexp) - sy0;
tmp = hmod1/tmp;
*dR = (*R)*hexp*tmp;
*ddR = (*dR)*(hexp-1.)*tmp;
}
break;
//LINEAR-EXPONENTIAL HARDENING
case H_LINEXP:
if(p<DP_MIN){
//printf("Should not be here: p: %.10e\n",p);
*R = 0;
*dR = hmod1 + hmod2*hexp;
*ddR = -hmod2*hexp*hexp;
}
else{
tmp = exp(-hexp*p);
*R = hmod1*p + hmod2*(1.-tmp);
*dR = hmod1 + hmod2*hexp*tmp;
*ddR = - hmod2*hexp*hexp*tmp;
}
break;
//LINEAR-EXPONENTIAL HARDENING FOR DRUCKER PRAGER
case H_LINEXP_DP:
if(p<DP_MIN){
//printf("Should not be here: p: %.10e\n",p);
*R = 0;
*dR = hmod1 + hmod2*hexp;
*ddR = -hmod2*hexp*hexp;
}
else{
tmp = exp(-hexp*p);
*R = hmod1*p + hmod2*(1.-tmp);
*dR = hmod1 + hmod2*hexp*tmp;
*ddR = - hmod2*hexp*hexp*tmp;
} break;
//POWER LAW HARDENING EXTRAPOLATED AFTER 16% DEFO TO MIMIC DIGIMAT TO ABAQUS
case H_POW_EXTRAPOL:
if(p<DP_MIN){
//printf("Should not be here: p: %.10e\n",p);
*R=0.;
if(hmod1<1){
*dR = 1e20;
*ddR = -1e20;
}
else{
*dR = hmod1;
*ddR = 0.;
}
}
else if(p<hp0)
{
*R=hmod1*pow(p,hexp);
*dR = (*R)*hexp/p;
*ddR = (*dR)*(hexp-1.)/p;
}
else if(p<10.*hp0)
{
*R=hmod1*pow(hp0,hexp)+hmod2*(p-hp0);
*dR = hmod2;
*ddR = 0.;
}
else
{
*R=hmod1*pow(hp0,hexp)+hmod2*(10.*hp0-hp0)+hmod2*(p-10.*hp0)/1000.;
*dR = hmod2/1000.;
*ddR = 0.;
}
break;
//POWER LAW HARDENING EXTRAPOLATED AFTER 16% DEFO TO MIMIC DIGIMAT TO ABAQUS DRUCKER PRAGER
case H_POW_EXTRAPOL_DP:
if(p<DP_MIN){
//printf("Should not be here: p: %.10e\n",p);
*R=0.;
if(hmod1<1){
*dR = 1e20;
*ddR = -1e20;
}
else{
*dR = hmod1;
*ddR = 0.;
}
}
else if(p<hp0)
{
*R=hmod1*pow(p,hexp);
*dR = (*R)*hexp/p;
*ddR = (*dR)*(hexp-1.)/p;
}
else if(p<10.*hp0)
{
*R=hmod1*pow(hp0,hexp)+hmod2*(p-hp0);
*dR = hmod2;
*ddR = 0.;
}
else
{
*R=hmod1*pow(hp0,hexp)+hmod2*(10.*hp0-hp0)+hmod2*(p-10.*hp0)/1000.;
*dR = hmod2/1000.;
*ddR = 0.;
}
break;
//POWER EXPONENTIAL HARDENING
case H_POW_EXP:
if(p<DP_MIN){
*R = 0.0;
*dR = 0.0;
if(hmod2<1.0){
*ddR = 1e20;
}
else{
*ddR = 0.;
}
}
else{
tmp = exp(-hexp*p);
*R = hmod1*pow(p,hmod2)*(1.-tmp);
*dR = hmod2*hmod1*pow(p,hmod2-1.)*(1.-tmp) + hexp*hmod1*pow(p,hmod2)*tmp;
*ddR = hmod2*(hmod2-1.0)*hmod1*pow(p,hmod2-2.)*(1.-tmp) + 2.0*hmod2*hexp*hmod1*pow(p,hmod2-1.0)*tmp - hexp*hexp*hmod1*pow(p,hmod2)*tmp;
}
break;
//POWER EXPONENTIAL HARDENING FOR DRUCKER PRAGAER
case H_POW_EXP_DP:
if(p<DP_MIN){
*R = 0.0;
*dR = 0.0;
if(hmod2<1.0){
*ddR = 1e20;
}
else{
*ddR = 0.;
}
}
else{
tmp = exp(-hexp*p);
*R = hmod1*pow(p,hmod2)*(1.-tmp);
*dR = hmod2*hmod1*pow(p,hmod2-1.)*(1.-tmp) + hexp*hmod1*pow(p,hmod2)*tmp;
*ddR = hmod2*(hmod2-1.0)*hmod1*pow(p,hmod2-2.)*(1.-tmp) + 2.0*hmod2*hexp*hmod1*pow(p,hmod2-1.0)*tmp - hexp*hexp*hmod1*pow(p,hmod2)*tmp;
}
break;
default:
printf("Bad choice of hardening law in j2hard: %d\n",htype);
break;
}
//printf("p: %lf\t, R: %.10e, \t dR: %.10e \t, ddR: %.10e \n",p, *R,*dRdp,*ddRdp);
}
/***********************************************
* COMPUTATION OF ALGORITHMIC TANGENT OPERATOR
* FROM STATE VARIABLES
* ********************************************/
void algoJ2(double E, double nu, double sy0, int htype, double hmod1, double hmod2, double hp0, double hexp,
double* strs, double p, double dp, double** Calgo){
int i,j;
double seqtr;
double R,dR,ddR;
double mu;
static double *N=NULL;
static double *devstrs=NULL;
static double **Cel=NULL;
static double **mat1=NULL;
static double **mat2=NULL; static double ** NoN=NULL;
static double ** Idev=NULL;
double h, tmp1, tmp2;
if(N==NULL)
{
mallocvector(&devstrs,6);
mallocvector(&N,6);
mallocmatrix(&Cel,6,6);
mallocmatrix(&NoN,6,6);
mallocmatrix(&mat1,6,6);
mallocmatrix(&mat2,6,6);
mallocmatrix(&Idev,6,6);
}
cleartens2(devstrs);
cleartens2(N);
cleartens4(Cel);
cleartens4(NoN);
cleartens4(mat1);
cleartens4(mat2);
cleartens4(Idev);
//deviatoric operator
for(i=0;i<3;i++){
for(j=0;j<3;j++){
Idev[i][j]=-1./3.;
}
Idev[i][i]+=1.;
Idev[i+3][i+3]=1.;
}
//reset tangent operators to zero
mu = E/(2.*(1.+nu));
compute_Hooke(E,nu,Cel);
dev(strs,devstrs);
j2hard (p, sy0, htype, hmod1, hmod2, hp0, hexp, &R, &dR, &ddR);
addtens2(1.5/(sy0+R),devstrs,0.,devstrs,N);
seqtr = 3*mu*dp + sy0 + R;
produit_tensoriel(N,N,NoN);
h = 3.*mu + dR;
tmp1 = 4.*mu*mu/h;
tmp2 = 4.*mu*mu*dp/seqtr;
addtens4(tmp1,NoN,0,NoN,mat1);
addtens4(3.*tmp2/2.,Idev,-tmp2,NoN,mat2);
addtens4(1.,Cel,-1.,mat1,Calgo);
addtens4(1.,Calgo,-1.,mat2,Calgo);
//free(devstrs);
//free(N);
//freematrix(Cel,6);
//freematrix(mat1,6);
//freematrix(mat2,6);
//freematrix(Idev,6);
//freematrix(NoN,6);
}
//ELASTO-PLASTIC CONSTITUTIVE BOX for second moment mehtod
int constbox_ep(double E, double nu, double sy0, int htype, double hmod1, double hmod2, double hp0, double hexp, double *DE, double *dstrn, double *strs_n, double eq2strs_n, double *eq2strs, double dstrseq_tr, double *dstrseq_trdDE, double *strs, double* pstrn_n, double* pstrn, double p_n, double *p, FY2nd *FY2){
int i,j,it;
double strseqtr2; //Variables for the elastic prediction
double R, dR, ddR, strseq2, mu_el; //Variables for the correction phase
double Dp, temp1, temp2;
double cp=0.0;
double tol = 1e-12; //tolerance on the residual
double res, res0; //residual
R=0.;
dR=0.;
ddR=0.;
mu_el = E/(2.*(1.+nu));
static bool initialized=false;
// initialize only once!
static double **Cel=NULL;
static double **Idev=NULL;
static double *dstrstr=NULL;
static double *dstrstr_dev=NULL;
static double *Str=NULL;
static double *Ntr=NULL;
static double *strs_ndev=NULL;
static double *mat1=NULL;
bool stay=true;
if(!initialized)
{
mallocmatrix(&Cel, 6, 6);
mallocmatrix(&Idev,6,6);
mallocvector(&dstrstr,6);
mallocvector(&dstrstr_dev,6);
mallocvector(&Str,6);
mallocvector(&Ntr,6);
mallocvector(&strs_ndev,6);
mallocvector(&mat1,6);
initialized = true;
}
cleartens4(Cel);
cleartens4(Idev);
cleartens2(dstrstr);
cleartens2(dstrstr_dev);
cleartens2(Str);
cleartens2(Ntr);
cleartens2(strs_ndev);
cleartens2(mat1);
//deviatoric operator
for(i=0;i<3;i++){
for(j=0;j<3;j++){
Idev[i][j]=-1./3.;
}
Idev[i][i]+=1.;
Idev[i+3][i+3]=1.;
}
//elastic predictor
compute_Hooke(E,nu,Cel);
//compute trial incremental stress
contraction42(Cel, dstrn, dstrstr);
contraction42(Idev, dstrstr, dstrstr_dev);
contraction42(Idev, strs_n, strs_ndev);
if(dstrseq_tr > 0.){
addtens2(1.5/dstrseq_tr, dstrstr_dev, 0., dstrstr_dev, Ntr);
}
addtens2(1.0, strs_n, 1.0, dstrstr, Str); //1. elastic increment
Dp=DP_MIN;
j2hard (p_n+Dp, sy0, htype, hmod1, hmod2, hp0, hexp, &R, &dR,&ddR);
strseq2= (sy0 + R)*(sy0 + R);
//1.1 trial (eq_stress)^2
addtens2(1.0, dstrstr_dev, -2.0*mu_el*Dp, Ntr, mat1);
int scheme = 1;
//if(contraction22(strs_ndev,Str)<=0.0){
// scheme = 1;
//}
if(scheme==1){
strseqtr2 = (dstrseq_tr-3.0*mu_el*Dp)*(dstrseq_tr-3.0*mu_el*Dp);
}
else if(scheme==2){
strseqtr2 = (dstrseq_tr-3.0*mu_el*Dp)*(dstrseq_tr-3.0*mu_el*Dp) + 3.0*contraction22(strs_ndev,mat1);
}
else{
strseqtr2 = eq2strs_n + (dstrseq_tr-3.0*mu_el*Dp)*(dstrseq_tr-3.0*mu_el*Dp) + 3.0*contraction22(strs_ndev,mat1);
}
FY2->f = strseqtr2 - strseq2;
//1.2 check yielding condition
if(FY2->f<=0.){
*eq2strs = strseqtr2;
copyvect(Str, strs,6);
copyvect(pstrn_n, pstrn,6);
*p = p_n;
Dp=0.;
// clean structure FY2
FY2->f=0.0;
FY2->dfdp=0.0;
for(i=0;i<6;i++){
FY2->dfdstrn_m[i]=0.0;
FY2->dfdstrn_c[i]=0.0;
FY2->dpdstrn_m[i]=0.0;
FY2->dpdstrn_c[i]=0.0;
}
}
//2. Plastic correction: return mapping algorithm
else{
it=0;
res=1.0;
res0=res;
while(res>tol && stay){
it++;
if(it>20){
printf(":( no convergence in return mapping algorithm, p: %lf\n",*p);
stay=false;
Dp=1.;
//return(1);
}
*p=p_n+Dp;
j2hard (*p, sy0, htype, hmod1, hmod2, hp0, hexp, &R, &dR,&ddR);
strseq2= (sy0 + R)*(sy0 + R);
addtens2(1.0, dstrstr_dev, -2.0*mu_el*Dp, Ntr, mat1);
if(scheme==1){
strseqtr2 = (dstrseq_tr-3.0*mu_el*Dp)*(dstrseq_tr-3.0*mu_el*Dp);
}
else if(scheme==2){
strseqtr2 = (dstrseq_tr-3.0*mu_el*Dp)*(dstrseq_tr-3.0*mu_el*Dp) + 3.0*contraction22(strs_ndev,mat1);
}
else{
strseqtr2 = eq2strs_n + (dstrseq_tr-3.0*mu_el*Dp)*(dstrseq_tr-3.0*mu_el*Dp) + 3.0*contraction22(strs_ndev,mat1); }
FY2->f = strseqtr2 - strseq2;
res = fabs(FY2->f)/strseq2; //mu_el;
if(it==1||res<res0){
res0=res;
// 2.1 newton-raphson iteration , compute Dp
if(scheme==1 ){
FY2->dfdp = -6.0*mu_el*((dstrseq_tr-3.0*mu_el*Dp))-2.0*(sy0 + R)*dR;
}
else{
FY2->dfdp = -6.0*mu_el*((dstrseq_tr-3.0*mu_el*Dp) + contraction22(strs_ndev, Ntr))-2.0*(sy0 + R)*dR;
}
cp = -FY2->f/FY2->dfdp;
temp1 = 0.99*(DP_MIN-Dp);
if(cp<temp1) {
//printf("new Dp below Dpmin : correction!\n");
cp=temp1;
}
Dp+=cp;
}
else{
//printf("residual is not decreasing: correction!\n");
cp=cp/2.;
Dp=Dp-cp;
res0=res;
}
//printf("j2plast iteration: %d, res: %.10e, p: %.10e\n",it,res,*p);
}
// 2.2 After obtainning the value of Dp, compute the out put variables
*eq2strs = (sy0 + R)*(sy0 + R);
addtens2(1., Str,-2.*mu_el*Dp, Ntr, strs);
addtens2(1.,pstrn_n, Dp, Ntr, pstrn);
//2.3 compute the derivatives of f and p
if(scheme==1){
temp1 = 0.0;
temp2 = 2.0*(dstrseq_tr-3.0*mu_el*Dp);
}
else{
temp1 = 6.0*mu_el*(1.0-3.0*mu_el*Dp/dstrseq_tr);
temp2 = 2.0*(dstrseq_tr-3.0*mu_el*Dp)+ 9.0*mu_el*Dp*contraction22(strs_ndev, dstrstr_dev)/dstrseq_tr/dstrseq_tr;
}
for(i=0;i<6;i++){
FY2->dfdstrn_m[i] = temp1*strs_ndev[i];
FY2->dpdstrn_m[i] = -FY2->dfdstrn_m[i]/FY2->dfdp;
FY2->dfdstrn_c[i] = temp2*dstrseq_trdDE[i];
FY2->dpdstrn_c[i] = -FY2->dfdstrn_c[i]/FY2->dfdp;
}
} //end case plastic correction
if(!stay) return(1); //no convergence
return(0);
}
/////**********************************************///////
//**** ELASTO-VISCOPLASTIC CONSTITUTIVE BOX ************
/////******* 07.2016, By Ling Wu **************///////
/////**********************************************///////
int constbox_evp(double E, double nu, double sy0, int htype, double hmod1, double hmod2, double hp0,
double hexp, int viscmodel, double vmod, double vstressy, double vexp, double *dstrn, double *strs_n, double *strs, double* pstrn_n, double* pstrn, double p_n, double *p, double *twomu_iso, double* dtwomu_iso, double** Calgo, double*** dCalgo,
double* dpdE, bool residual, double dt)
{
int i,j,k,it,error;
double strseqtr, strseqtrNor0, trstrstr; //Variables for the elastic prediction
double R, dR, ddR, strseq, G; //Variables for the correction phase
double gv, hv, dgdf, dhvdSeq, dhvdp;
double kp;
double Dp;
double h, tmp1, tmp2;
double cp=0.;
double tol = 1.e-12; //tolerance on the residual
double res, res0; //residual
double mu_iso, k_iso;
double G3, h3;
double sq2;
double eq2strs_n;
sq2 = sqrt(2.);
R=0.;
dR=0.;
ddR=0.;
G = E/(2.*(1.+nu));
// k_iso = E/(3.*(1.-(2.*nu)));
G3 = G*G*G;
static bool initialized=false;
// initialize only once!
static double **Cel=NULL;
static double **NoN=NULL;
static double **Idev=NULL;
static double **mat1=NULL;
static double *strstr=NULL;
static double *strstrNor0=NULL;
static double *Str=NULL;
static double *Ntr=NULL;
static double *dstrs=NULL;
bool stay=true;
if(!initialized)
{
mallocmatrix(&Cel, 6, 6);
mallocmatrix(&NoN,6,6);
mallocmatrix(&Idev,6,6);
mallocmatrix(&mat1,6,6);
mallocvector(&strstr,6);
mallocvector(&strstrNor0,6);
mallocvector(&Str,6);
mallocvector(&Ntr,6);
mallocvector(&dstrs,6);
initialized = true;
}
cleartens4(Cel);
cleartens4(NoN);
cleartens4(Idev);
cleartens4(mat1);
cleartens2(strstr);
cleartens2(strstrNor0);
cleartens2(Str);
cleartens2(Ntr);
cleartens2(dstrs);
//deviatoric operator
for(i=0;i<3;i++){
for(j=0;j<3;j++){
Idev[i][j]=-1./3.; }
Idev[i][i]+=1.;
Idev[i+3][i+3]=1.;
}
//reset tangent operators to zero
for(i=0;i<6;i++){
for(j=0;j<6;j++){
Calgo[i][j]=0.;
}
}
//elastic predictor
compute_Hooke(E,nu,Cel);
//compute trial stress
contraction42(Cel, dstrn, dstrs);
addtens2(1., strs_n, 1., dstrs, strstrNor0);
addtens2(1., strs_n, 1., dstrs, strstr);
strseqtrNor0 = vonmises(strstrNor0);
// strseqtrNor0 = vonmises(dstrs);
strseqtr = vonmises(strstr);
//normal to the yield surface
dev(strstrNor0, Str);
// dev(dstrs, Str);
if(strseqtrNor0 > 0.){
addtens2(1.5/strseqtrNor0, Str, 0., Str, Ntr);
produit_tensoriel(Ntr,Ntr,NoN);
}
//elastic increment
Dp=DP_MIN;
j2hard (p_n+Dp, sy0, htype, hmod1, hmod2, hp0, hexp, &R, &dR,&ddR);
strseq= sy0 + R;
kp = strseq + 3.*G*Dp - strseqtr;
if(kp>=0.){
//printf("elastic: !");
copyvect(strstr, strs, 6);
copyvect(pstrn_n, pstrn,6);
*p = p_n;
Dp=0.;
copymatr(Cel,Calgo,6,6);
for(i=0;i<6;i++){
for(j=0;j<6;j++){
for(k=0;k<6;k++){
dCalgo[i][j][k]=0.;
}
}
}
//compute dpdE accumulate plastice strain p derivative w.r.t. strain increment
cleartens2(dpdE);
}
//ViscoPlastic correction: return mapping algorithm
else{
Dp = DP_MIN; //correction must be such that Dp >= DP_MIN
it=0;
res=1.;
res0=res;
while(res>tol && stay){
it++;
if(it>20){ printf("no convergence in return mapping algorithm, p: %lf\n",*p);
//return(1);
stay=false;
Dp=1.;
}
*p=p_n+Dp;
j2hard (*p, sy0, htype, hmod1, hmod2, hp0, hexp, &R, &dR,&ddR);
strseq = strseqtr-3.*G*Dp;
solve_gv(viscmodel, vmod, vstressy, vexp, G, strseq, sy0, dt, R, dR, ddR, &gv, &hv, &dgdf, &dhvdSeq, &dhvdp,eq2strs_n);
kp = Dp - gv*dt;
res = fabs(kp);
h=hv;
if(it==1||res<res0){
res0=res;
cp = - kp/(h*dgdf*dt);
tmp1 = 0.99*(DP_MIN-Dp);
if(cp<tmp1) {
//printf("new Dp below Dpmin : correction!\n");
cp=tmp1;
}
Dp+=cp;
}
else{
//printf("residual is not decreasing: correction!\n");
cp=cp/2.;
Dp=Dp-cp;
res0=res;
}
//printf("j2plast iteration: %d, res: %.10e, p: %.10e\n",it,res,*p);
}
//printf("Number of iterations in radial return algorithm: %d\n",it);
addtens2(1., strstr,-2.*G*(Dp), Ntr, strs);
addtens2(1.,pstrn_n, Dp, Ntr, pstrn);
h3 = h*h*h;
tmp1 = 4.*G*G/h;
tmp2 = 4.*G*G*(*p-p_n)/strseqtr;
addtens4(1.,Cel,(-tmp1+tmp2),NoN, Calgo);
addtens4(1.,Calgo,-1.5*tmp2,Idev,Calgo);
addtens4(1.5,Idev,-1.,NoN,mat1);
for(i=0;i<6;i++){
for(j=0;j<6;j++){
for(k=0;k<6;k++){
dCalgo[i][j][k] = (8.*G3/(h*h))*(dhvdSeq*(1.0-3.0*G/h) + dhvdp/h)*NoN[i][j]*Ntr[k] - (16.*G3/(h*strseqtr))*Ntr[i]*mat1[j][k] - mat1[i][j]*8.*G3*Ntr[k]*((1./(h*strseqtr))-(Dp/(strseqtr*strseqtr))) + (16.*G3*Dp/(strseqtr*strseqtr))*Ntr[i]*mat1[j][k];
}
}
}
//compute dpdE accumulate plastice strain p derivative w.r.t. strain increment
addtens2((2.*G/h), Ntr, 0, Ntr, dpdE);
} //end case plastic correction
//compute isoJ2 tangent operator and its derivative w.r.t. strain increment if(Dp<DP_MIN){ //that is, Dp=0
mu_iso=G;
for(i=0;i<6;i++){
dtwomu_iso[i]=0.;
}
}
else {
mu_iso = G*(1.-(3.*G/h));
for(i=0;i<6;i++){
dtwomu_iso[i] = 12.*pow(G/h,3.)*ddR*Ntr[i];
}
}
*twomu_iso = 2.*mu_iso;
// no free as initialized once
// freematrix(Cel, 6);
// freematrix(Idev,6);
// freematrix(NoN,6);
// freematrix(mat1,6);
if(!stay) return(1); //no convergence
return(0);
}
//ELASTO-VISCOPLASTIC CONSTITUTIVE BOX for second moment mehtod
int constbox_evp(double E, double nu, double sy0, int htype, double hmod1, double hmod2, double hp0, double hexp, int viscmodel, double vmod, double vstressy, double vexp, double *DE, double *dstrn, double *strs_n, double eq2strs_n, double *eq2strs, double dstrseq_tr, double *dstrseq_trdDE, double *strs, double* pstrn_n, double* pstrn, double p_n, double *p, FY2nd *FY2, double dt, ELcc *LCC){
int i,j,it;
double strseqtr; //Variables for the elastic prediction
double gv, hv, dgdf, dhvdSeq, dhvdp;
double R, dR, ddR, strseq, mu_el; //Variables for the correction phase
double kp;
double Dp, temp0, temp1, temp2;
double cp=0.0;
double tol = 1e-12; //tolerance on the residual
double res, res0; //residual
R=0.;
dR=0.;
ddR=0.;
mu_el = E/(2.*(1.+nu));
static bool initialized=false;
// initialize only once!
static double **Cel=NULL;
static double **Idev=NULL;
static double **I=NULL;
//static double **ImA=NULL;
static double *dstrstr=NULL;
static double *dstrstr_dev=NULL;
static double *Str=NULL;
static double *Ntr=NULL;
static double *strs_ndev=NULL;
//static double *Edstrn_i=NULL;
//static double *Edstrn_m=NULL;
static double *mat1=NULL;
//static double *Edstrstr=NULL;
//static double *Edstrstr_dev=NULL;
//static double *P;
bool stay=true;
if(!initialized)
{
mallocmatrix(&Cel, 6, 6);
mallocmatrix(&Idev,6,6);
mallocmatrix(&I,6,6);
//mallocmatrix(&ImA,6,6); mallocvector(&dstrstr,6);
mallocvector(&dstrstr_dev,6);
mallocvector(&Str,6);
mallocvector(&Ntr,6);
mallocvector(&strs_ndev,6);
//mallocvector(&Edstrn_i,6);
//mallocvector(&Edstrn_m,6);
mallocvector(&mat1,6);
//mallocvector(&Edstrstr,6);
//mallocvector(&Edstrstr_dev,6);
//mallocvector(&P,6);
initialized = true;
}
cleartens4(Cel);
cleartens4(Idev);
cleartens4(I);
//cleartens4(ImA);
cleartens2(dstrstr);
cleartens2(dstrstr_dev);
cleartens2(Str);
cleartens2(Ntr);
cleartens2(strs_ndev);
//cleartens2(Edstrn_i);
//cleartens2(Edstrn_m);
cleartens2(mat1);
//cleartens2(Edstrstr);
//cleartens2(Edstrstr_dev);
//cleartens2(P);
//deviatoric operator
for(i=0;i<3;i++){
for(j=0;j<3;j++){
Idev[i][j]=-1./3.;
}
Idev[i][i]+=1.;
Idev[i+3][i+3]=1.;
}
for(i=0;i<6;i++){
I[i][i] = 1.;
}
//elastic predictor
compute_Hooke(E,nu,Cel);
//compute trial incremental stress
contraction42(Cel, dstrn, dstrstr);
contraction42(Idev, dstrstr, dstrstr_dev);
contraction42(Idev, strs_n, strs_ndev);
//double v1= LCC->v1;
//double v0 = 1.0 - v1;
//contraction42(LCC->A, DE, Edstrn_i); //unloading strain of inclusion at tn
//addtens2 (1./v0, DE, -1.*v1/v0, Edstrn_i, Edstrn_m); //unloading strain of matrix at tn
//addtens4(1.0/v0, I, -1.0*v1/v0, LCC->A, ImA);
//compute trial incremental stress
//contraction42(Cel, dstrn, Edstrstr);
//contraction42(Idev, Edstrstr, Edstrstr_dev);
//contraction42(Idev, strs_n, strs_ndev);
//contraction42(ImA, strs_ndev,P);
//addtens2(1.0, strs_n, 1.0, Edstrstr, Str); addtens2(1.0, strs_n, 1.0, dstrstr, Str);
//1. elastic increment
Dp=DP_MIN;
j2hard (p_n+Dp, sy0, htype, hmod1, hmod2, hp0, hexp, &R, &dR,&ddR);
//1.1 trial (eq_stress)^2
strseqtr = dstrseq_tr;
int soft_inclusions=0;
if(soft_inclusions==1)
{
// for soft inclusion, a modification is applied-----------
temp0 = (dstrseq_tr*dstrseq_tr)+6.0*mu_el*contraction22(strs_ndev, dstrn)+ 1.5*contraction22(strs_ndev,strs_ndev);
strseqtr = fabs(pow(temp0,0.5)-3.0*mu_el*Dp);
}
if(dstrseq_tr > 0.){
addtens2(1.5/dstrseq_tr, dstrstr_dev, 0., strs_ndev, Ntr);
//old version
//addtens2(1.5/dstrseq_tr, dstrstr_dev, 1.5, strs_ndev, Ntr);
//addtens2(1.5/strseqtr, dstrstr_dev, 1.5/strseqtr, strs_ndev, Ntr);
//addtens2(1.5/strseqtr, dstrstr_dev, 1.5/strseqtr, strs_ndev, Ntr);
//----------------------------------------------------------------------
}
FY2->f = strseqtr - (sy0 + R);
//1.2 check yielding condition
if(FY2->f<=0.){
*eq2strs = strseqtr*strseqtr;
copyvect(Str, strs,6);
copyvect(pstrn_n, pstrn,6);
*p = p_n;
Dp=0.;
// clean structure FY2
FY2->f=0.0;
FY2->dfdp=0.0;
for(i=0;i<6;i++){
FY2->dfdstrn_m[i]=0.0;
FY2->dfdstrn_c[i]=0.0;
FY2->dpdstrn_m[i]=0.0;
FY2->dpdstrn_c[i]=0.0;
}
}
//2. Plastic correction: return mapping algorithm
else{
it=0;
res=1.0;
res0=res;
while(res>tol && stay){
it++;
if(it>20){
printf("no convergence in return mapping algorithm, p: %lf\n",*p);
stay=false;
Dp=1.;
//return(1);
}
*p=p_n+Dp;
j2hard (*p, sy0, htype, hmod1, hmod2, hp0, hexp, &R, &dR,&ddR);
strseq = strseqtr-3.0*mu_el*Dp;
solve_gv(viscmodel, vmod, vstressy, vexp, mu_el, strseq, sy0, dt, R, dR, ddR, &gv, &hv, &dgdf, &dhvdSeq, &dhvdp,eq2strs_n);
kp = Dp - gv*dt;
res = fabs(kp);
if(it==1||res<res0){
res0=res;
cp = - kp/(hv*dgdf*dt);
temp1 = 0.99*(DP_MIN-Dp);
if(cp<temp1) {
//printf("new Dp below Dpmin : correction!\n");
cp=temp1;
}
Dp+=cp;
}
else{
//printf("residual is not decreasing: correction!\n");
cp=cp/2.;
Dp=Dp-cp;
res0=res;
}
//printf("j2plast iteration: %d, res: %.10e, p: %.10e\n",it,res,*p);
}
// 2.2 After obtainning the value of Dp, compute the out put variables
*eq2strs = strseq*strseq;
addtens2(1., Str,-2.*mu_el*Dp, Ntr, strs);
addtens2(1.,pstrn_n, Dp, Ntr, pstrn);
//2.3 compute the derivatives of f and p
for(i=0;i<6;i++){
FY2->dpdstrn_m[i] = 0.0;
if(soft_inclusions==1)
{
// for soft inclusion, a modification is applied-----------
FY2->dpdstrn_m[i] = strs_ndev[i]*3.0*mu_el/(strseqtr)/hv;
//---------------------
}
FY2->dpdstrn_c[i] = dstrseq_trdDE[i]/hv;
}
} //end case plastic correction
if(!stay) return(1); //no convergence
return(0);
}
//****************************************************
/**********J2 viscoplasticity for large deformation***/
//******************************************************
int constbox_evp(double E, double nu, double sy0, int htype, double hmod1, double hmod2, double hp0,
double hexp, int viscmodel, double vmod, double vstressy, double vexp, double *strs_tr, double *Rstrs_n, double *strs, double* pstrn_n,
double* pstrn, double p_n, double *p, double** Calgo, double*** dCalgo, double* dpdE, bool residual, double dt)
{
int i,j,k,it,error;
double strseqtr,strseqtr0; //Variables for the elastic prediction
double R, dR, ddR, strseq, G; //Variables for the correction phase
double gv, hv, dgdf, dhvdSeq, dhvdp;
double kp;
double Dp;
double h, tmp1, tmp2;
double cp=0.0;
double tol = 1.e-12; //tolerance on the residual
double res, res0; //residual
double G3, h3;
double sq2;
double eq2strs_n;
bool stay=true;
sq2 = sqrt(2.);
R=0.;
dR=0.;
ddR=0.;
G = E/(2.*(1.+nu));
// k_iso = E/(3.*(1.-(2.*nu)));
G3 = G*G*G;
static bool initialized=false;
// initialize only once!
static double **Cel;
static double **NoN;
static double **Idev;
static double **mat1;
static double *Str;
static double *Ntr;
static double *strstrNor0;
if(!initialized)
{
mallocmatrix(&Cel, 6, 6);
mallocmatrix(&NoN,6,6);
mallocmatrix(&Idev,6,6);
mallocmatrix(&mat1,6,6);
mallocvector(&Str,6);
mallocvector(&Ntr,6);
mallocvector(&strstrNor0,6);
initialized = true;
}
cleartens4(Cel);
cleartens4(NoN);
cleartens4(Idev);
cleartens4(mat1);
cleartens2(Str);
cleartens2(Ntr);
cleartens2(strstrNor0);
//deviatoric operator
for(i=0;i<3;i++){
for(j=0;j<3;j++){
Idev[i][j]=-1./3.;
}
Idev[i][i]+=1.;
Idev[i+3][i+3]=1.;
}
//reset tangent operators to zero
for(i=0;i<6;i++){
for(j=0;j<6;j++){
Calgo[i][j]=0.;
}
}
//elastic predictor
compute_Hooke(E,nu,Cel);
strseqtr = vonmises(strs_tr);
//normal to the yield surface
//dev(strs_tr, Str);
//normal to the residual stress
addtens2(1.0, strs_tr, -1.0, Rstrs_n, strstrNor0);
dev(strstrNor0, Str);
strseqtr0 = vonmises(strstrNor0);
if(strseqtr > 0.){
addtens2(1.5/strseqtr0, Str, 0., Str, Ntr);
produit_tensoriel(Ntr,Ntr,NoN);
}
//elastic increment
Dp=DP_MIN;
j2hard (p_n+Dp, sy0, htype, hmod1, hmod2, hp0, hexp, &R, &dR,&ddR);
strseq= sy0 + R;
kp = strseq + 3.*G*Dp - strseqtr;
if(kp>=0.){
//printf("elastic: !");
copyvect(strs_tr, strs, 6);
copyvect(pstrn_n, pstrn,6);
*p = p_n;
Dp=0.;
copymatr(Cel,Calgo,6,6);
for(i=0;i<6;i++){
for(j=0;j<6;j++){
for(k=0;k<6;k++){
dCalgo[i][j][k]=0.;
}
}
}
//compute dpdE accumulate plastice strain p derivative w.r.t. strain increment
cleartens2(dpdE);
}
//ViscoPlastic correction: return mapping algorithm
else{
Dp = DP_MIN; //correction must be such that Dp >= DP_MIN
it=0;
res=1.;
res0=res;
while(res>tol && stay){
it++;
if(it>20){
printf("no convergence in return mapping algorithm, p: %lf\n",*p);
stay=false;
Dp=1.;
//return(1);
}
*p=p_n+Dp;
j2hard (*p, sy0, htype, hmod1, hmod2, hp0, hexp, &R, &dR,&ddR);
strseq = strseqtr-3.*G*Dp;
solve_gv(viscmodel, vmod, vstressy, vexp, G, strseq, sy0, dt, R, dR, ddR, &gv, &hv, &dgdf, &dhvdSeq, &dhvdp,eq2strs_n);
kp = Dp - gv*dt;
res = fabs(kp);
h=hv;
if(it==1||res<res0){
res0=res;
cp = - kp/(h*dgdf*dt);
tmp1 = 0.99*(DP_MIN-Dp);
if(cp<tmp1) {
//printf("new Dp below Dpmin : correction!\n");
cp=tmp1;
}
Dp+=cp;
}
else{
//printf("residual is not decreasing: correction!\n");
cp=cp/2.;
Dp=Dp-cp;
res0=res;
}
//printf("j2plast iteration: %d, res: %.10e, p: %.10e\n",it,res,*p);
}
//printf("Number of iterations in radial return algorithm: %d\n",it);
addtens2(1., strs_tr,-2.*G*(Dp), Ntr, strs);
addtens2(1.,pstrn_n, Dp, Ntr, pstrn);
h3 = h*h*h;
tmp1 = 4.*G*G/h;
tmp2 = 4.*G*G*(*p-p_n)/strseqtr;
addtens4(1.,Cel,(-tmp1+tmp2),NoN, Calgo);
addtens4(1.,Calgo,-1.5*tmp2,Idev,Calgo);
addtens4(1.5,Idev,-1.,NoN,mat1);
for(i=0;i<6;i++){
for(j=0;j<6;j++){
for(k=0;k<6;k++){
dCalgo[i][j][k] = (8.*G3/(h*h))*(dhvdSeq*(1.0-3.0*G/h) + dhvdp/h)*NoN[i][j]*Ntr[k] - (16.*G3/(h*strseqtr))*Ntr[i]*mat1[j][k] - mat1[i][j]*8.*G3*Ntr[k]*((1./(h*strseqtr))-(Dp/(strseqtr*strseqtr))) + (16.*G3*Dp/(strseqtr*strseqtr))*Ntr[i]*mat1[j][k];
}
}
}
//compute dpdE accumulate plastice strain p derivative w.r.t. strain increment
addtens2((2.*G/h), Ntr, 0, Ntr, dpdE);
} //end case plastic correction
// no free as initialized once
// freematrix(Cel, 6);
// freematrix(Idev,6);
// freematrix(NoN,6);
// freematrix(mat1,6);
if(!stay) return(1); //no convergence
return(0);
}
/***********************************************************
* SOLVE VISCOPLASTIC FUNCTION gv *
**********************************************************/
void solve_gv(int viscmodel, double vmod, double vstressy, double vexp, double mu, double strseq, double yldstrs, double dt, double R, double dR, double ddR, double *gv, double *hv, double *dgdf, double *dhvdSeq, double *dhvdp, double eq2strs_n){
double f;
double tmp1,tmp2, dgdR;
double ddgdSdR, ddgdR;
f = strseq - yldstrs - R;
if(f < 0.0){
printf("f < 0 Should be elastic, f = %.10e\n",f);
return;
}
switch(viscmodel){
case V_NORTON:
*gv = vmod*pow(f/yldstrs,vexp);
tmp1 = vmod*vexp/yldstrs*pow(f/yldstrs,vexp-1.0);
*dgdf = tmp1;
*hv = 1.0/(tmp1*dt)+3.0*mu + dR;
tmp2 = vmod*vexp*(vexp-1.0)/(yldstrs*yldstrs)*pow(f/yldstrs,vexp-2.0); //ddgvdSeq
*dhvdSeq = -1.0/(tmp1*tmp1*dt)*tmp2;
*dhvdp = 1.0/(tmp1*tmp1*dt)*tmp2*dR + ddR;
break;
//PERZYNA law with drag stress implemented in Abaqus
case V_PERZYNA:
*gv = vmod*pow(f/(yldstrs+R),vexp);
tmp1 = vmod*vexp/(yldstrs+R)*pow(f/(yldstrs+R),vexp-1.0);
*dgdf = tmp1;
dgdR = -tmp1-vmod*vexp/(yldstrs+R)*pow(f/(yldstrs+R),vexp);
*hv = 1.0/(tmp1*dt) + 3.0*mu - dgdR/tmp1*dR ;
tmp2 = vmod*vexp*(vexp-1.0)/((yldstrs+R)*(yldstrs+R))*pow(f/(yldstrs+R),vexp-2.0); //ddgvdSeq
ddgdSdR = -tmp2*(1.0 + vexp*f/(yldstrs+R)/(vexp-1.0));
ddgdR = tmp2*(1.0 + 2.0*vexp*f/(yldstrs+R)/(vexp-1.0) + (vexp+1.0)*f/(yldstrs+R)*f/(yldstrs+R)/(vexp-1.0));
*dhvdSeq = 1.0/(tmp1*tmp1)*tmp2*(-1.0/dt + dgdR*dR) - ddgdSdR*dR/tmp1;
*dhvdp = (-1.0/dt + dgdR*dR)/(tmp1*tmp1)*ddgdSdR*dR - (ddgdR*dR*dR+ dgdR*ddR)/tmp1;
break;
//Power law
case V_POWER:
*gv = vmod*pow(strseq/(yldstrs+R),vexp);
tmp1 = vmod*vexp/(yldstrs+R)*pow(strseq/(yldstrs+R),vexp-1.0);
*dgdf = tmp1;
dgdR = -vmod*vexp/(yldstrs+R)*pow(strseq/(yldstrs+R),vexp);
*hv = 1.0/(tmp1*dt) + 3.0*mu - dgdR/tmp1*dR ;
tmp2 = vmod*vexp*(vexp-1.0)/((yldstrs+R)*(yldstrs+R))*pow(strseq/(yldstrs+R),vexp-2.0); //ddgvdSeq
ddgdSdR = -tmp1*vexp/(yldstrs+R);
ddgdR = vmod*vexp*(vexp+1.0)/(yldstrs+R)/(yldstrs+R)*pow(strseq/(yldstrs+R),vexp);
*dhvdSeq = 1.0/(tmp1*tmp1)*tmp2*(-1.0/dt + dgdR*dR) - ddgdSdR*dR/tmp1;
*dhvdp = (-1.0/dt + dgdR*dR)/(tmp1*tmp1)*ddgdSdR*dR - (ddgdR*dR*dR+ dgdR*ddR)/tmp1;
break;
case V_PRANDTL:
double Beta;
Beta = 10.0;
*gv = vmod*pow(sinh(f/Beta),vexp);
tmp1 = vmod*vexp/Beta*pow(sinh(f/Beta),vexp-1.0)*cosh(f/Beta);
*dgdf = tmp1;
*hv = 1.0/(tmp1*dt)+3.0*mu + dR;
tmp2 = vmod*vexp/Beta/Beta*((vexp-1.0)*pow(sinh(f/Beta),vexp-2.0)*cosh(f/Beta)*cosh(f/Beta)+
pow(sinh(f/Beta),vexp)); //ddgvdSeq
*dhvdSeq = -1.0/(tmp1*tmp1*dt)*tmp2;
*dhvdp = 1.0/(tmp1*tmp1*dt)*tmp2*dR + ddR;
break;
//Power law with hardening: vmod is actually dot eps0/sy0^n
case V_POWER_NO_HARDENING:
*gv = vmod*pow(strseq/vstressy,vexp);
tmp1 = vmod/vstressy*vexp*pow(strseq/vstressy,vexp-1.0);
*dgdf = tmp1;
dgdR = 0.;
*hv = 1.0/(tmp1*dt) + 3.0*mu ;
tmp2 = vmod*vexp*(vexp-1.0)/vstressy/vstressy*pow(strseq/vstressy,vexp-2.0); //ddgvdSeq
ddgdSdR = 0.0;
ddgdR = 0.0;
if(vexp == 1.0){
*dhvdSeq = 0.0;
*dhvdp =0.0;
} else{
*dhvdSeq = -1.0/(dt*tmp1*tmp1)*tmp2;
*dhvdp = 0.0;
}
break;
default:
printf("Invalid viscosity law: %d\n", viscmodel);
break;
}
return;
}
#endif