FEM_NHK_11.17
目录
bending_of_five_elements.m
%
% Bending of five elements: linear elastic
%
% Nodal coordinates (mm)
XYZ=[0 0 0 ;1 0 0 ;1 1 0 ;0 1 0 ;
0 0 2 ;1 0 2 ;1 1 2 ;0 1 2 ;
0 0 4 ;1 0 4 ;1 1 4 ;0 1 4 ;
0 0 6 ;1 0 6 ;1 1 6 ;0 1 6 ;
0 0 7 ;1 0 7 ;1 1 8 ;0 1 8 ;
0 0 10;1 0 10;1 1 10;0 1 10];
%
% Element connectivity
LE=[ 1 2 3 4 5 6 7 8;
5 6 7 8 9 10 11 12;
9 10 11 12 13 14 15 16;
13 14 15 16 17 18 19 20;
17 18 19 20 21 22 23 24];
%
% External forces [Node, DOF, Value]
EXTFORCE=[21 2 1.0E5; 22 2 1.0E5; 23 2 1.0E5; 24 2 1.0E5];
%
% Prescribed displacements [Node, DOF, Value]
SDISPT=[1 1 0;1 2 0;1 3 0;
2 1 0;2 2 0;2 3 0;
3 1 0;3 2 0;3 3 0;
4 1 0;4 2 0;4 3 0];
%
% Material properties
% MID:0(Linear elastic) PROP=[LAMBDA MNU]
%MID=0;
MID=31;
%
%for elastic(MID=0)
%E=2E11; NU=0.3; LAMBDA=E*NU/((1+NU)*(1-2*NU)); MU=E/(2*(1+NU));
%PROP=[LAMBDA MU];
%for elasto plastic mutiply(MID=31)
Young = 24000; nu=0.2; mu=Young/2/(1+nu); lambda=nu*Young/((1+nu)*(1-2*nu));
beta = 0; H = 1000; sY = 200*sqrt(3);
mp = [lambda mu beta H sY];
%
% Load increments [Start End Increment InitialFactor FinalFactor]
%TIMS=[0.0 1.0 1.0 0.0 1.0]';
TIMS=[0.0 1.0 0.1 0.1 0.2]';
%
% Set program parameters
ITRA=30; ATOL=1.0E5; NTOL=6; TOL=1E-6;
%
% Calling main function
NOUT = fopen('output.txt','w');
NLFEA(ITRA, TOL, ATOL, NTOL, TIMS, NOUT, MID, PROP, EXTFORCE, SDISPT, XYZ, LE);
fclose(NOUT);
cntelm2d.m
function [FORCE, STIFF] = cntelm2d(OMEGAN, OMEGAT, CFRI, ELXY, ELXYP, LTAN)
%***********************************************************************
% SEARCH CONTACT POINT AND RETURN STIFFNESS AND RESIDUAL FORCE
% IF CONTACTED FOR NORMAL CONTACT
%***********************************************************************
%
ZERO = 0.D0; ONE = 1.D0; EPS = 1.E-6; P05 = 0.05; FORCE=[]; STIFF=[];
XT = ELXY(:,3)-ELXY(:,2); XLEN = norm(XT);
if XLEN < EPS, return; end
XTP = ELXYP(:,3)-ELXYP(:,2); XLENP = norm(XTP);
%
% UNIT NORMAL AND TANGENTIAL VECTOR
XT = XT/XLEN;
XTP = XTP/XLENP;
XN = [-XT(2); XT(1)];
%
% NORMAL GAP FUNCTION Gn = (X_s - X_1).N
GAPN = (ELXY(:,1)-ELXY(:,2))'*XN;
%
% CHECK IMPENETRATION CONDITION
if (GAPN >= ZERO) || (GAPN <= -XLEN), return; end
%
% NATURAL COORDINATE AT CONTACT POINT
ALPHA = (ELXY(:,1) - ELXY(:,2))'*XT/XLEN;
ALPHA0 = ((ELXYP(:,1)-ELXYP(:,2))'*XTP)/XLENP;
%
% OUT OF SEGMENT
if (ALPHA > ONE+P05) || (ALPHA < -P05), return; end
%
% CONTACT OCCURS IN THIS SEGMENT
XLAMBN = -OMEGAN*GAPN;
XLAMBT = 0;
LFRIC = 1; if CFRI == 0, LFRIC = 0; end
if LFRIC
GAPT = (ALPHA - ALPHA0)*XLENP;
XLAMBT = -OMEGAT*GAPT;
FRTOL = XLAMBN*CFRI;
LSLIDE = 0;
if abs(XLAMBT) > FRTOL
LSLIDE = 1;
XLAMBT = -FRTOL*SIGN(ONE,GAPT);
end
end
%
% DEFINE VECTORS
NN = [XN; -(ONE-ALPHA)*XN; -ALPHA*XN];
TT = [XT; -(ONE-ALPHA)*XT; -ALPHA*XT];
PP = [ZERO; ZERO; -XN; XN];
QQ = [ZERO; ZERO; -XT; XT];
CN = NN - GAPN*QQ/XLEN;
CT = TT + GAPN*PP/XLEN;
%
% CONTACT FORCE
FORCE = XLAMBN*CN + XLAMBT*CT;
%
% FORM STIFFNESS
if LTAN
STIFF = OMEGAN*(CN*CN');
if LFRIC
TMP1 = -CFRI*OMEGAN*SIGN(ONE,GAPT);
TMP2 = -XLAMBT/XLEN;
if LSLIDE
STIFF = STIFF + TMP1*(CT*CN') + TMP2*(CN*PP'+PP*CN'-CT*QQ'-QQ*CT');
else
STIFF = STIFF + OMEGAT*(CT*CT') + TMP2*(CN*PP'+PP*CN'-CT*QQ'-QQ*CT');
end
end
end
end
cntelm3d.m
function [FORCE, STIFF] = cntelm3d(OMEGAN, ELXY, LTAN)
%***********************************************************************
% SEARCH CONTACT POINT AND RETURN STIFFNESS AND RESIDUAL FORCE
% IF CONTACTED FOR NORMAL CONTACT
%***********************************************************************
%
EPS=1.E-6; TL1=2; TL2=0.1; TL3=1.01; FORCE=[]; STIFF=[];
%
% COMPUTE TWO TANGENT VECTORS & APPROXIMATED CONTACT POINT AT THE CENTER
[T1, T2, XS] = CUTL(0,0,ELXY);
XN = cross(T1, T2); XN=XN/norm(XN);
XI = (XS'*T1)/(2*norm(T1)^2);
ETA = (XS'*T2)/(2*norm(T2)^2);
GN = XN'*XS;
XX=(ELXY(:,2)-ELXY(:,3)+ELXY(:,4)-ELXY(:,5))/4;
%
% IF NATURAL COORD. IS OUT OF BOUND, NO CONTACT IS ASSUMED
if((XI<-TL1)||(XI>TL1)||(ETA<-TL1)||(ETA>TL1)||GN>TL2), return; end
%
% FIND EXACT CONTACT POINT THROUGH NEWTON-RAPHSON METHOD
for ICOUNT=1:20
[T1, T2, XS] = CUTL(ETA,XI,ELXY);
A=[-T1'*T1, XS'*XX-T2'*T1; XS'*XX-T2'*T1, -T2'*T2];
B=[-XS'*T1; -XS'*T2];
DXI=A\B;
XI=XI+DXI(1); ETA=ETA+DXI(2);
if(norm(DXI)<EPS), break; end
end
%
% CHECK FOR REFERENCE COORDINATE WITHIN RANGE
if((XI<-TL3)||(XI>TL3)||(ETA<-TL3)||(ETA>TL3)), return; end
%
% NORMAL GAP FUNCTION
XN = cross(T1, T2); XN=XN/norm(XN);
GN = XN'*XS;
if GN>0, return; end
%
% CONTACT FORCE
FORCE = -OMEGAN*GN*XN;
%
% FORM STIFFNESS (NONFRICTION)
if LTAN, STIFF = OMEGAN*(XN*XN'); end
end
function [T1, T2, XS] = CUTL(ETA,XI,ELXY)
%***********************************************************************
% COMPUTE COORD. OF CENTEROID AND TWO TANGENT VECTORS
%***********************************************************************
XNODE=[0 -1 1 1 -1; 0 -1 -1 1 1];
T1 = zeros(3,1); T2 = zeros(3,1); XC = zeros(3,1); XS = zeros(3,1);
for J = 1:3
T1(J) = sum(XNODE(1,2:5).*(1+ETA*XNODE(2,2:5)).*ELXY(J,2:5)./4);
T2(J) = sum(XNODE(2,2:5).*(1+XI *XNODE(1,2:5)).*ELXY(J,2:5)./4);
XC(J) = sum((1+XI*XNODE(1,2:5)).*(1+ETA*XNODE(2,2:5)).*ELXY(J,2:5)./4);
XS(J) = ELXY(J,1) - XC(J);
end
end
combHard.m
%
% Linear combined isotropic/kinamtic hardening model
%
function [stress, alpha, ep]=combHard(mp,D,deps,stressN,alphaN,epN)
% Inputs:
% mp = [lambda, mu, beta, H, Y0];
% D = elastic stiffness matrix
% stressN = [s11, s22, s33, t12, t23, t13];
% alphaN = [a11, a22, a33, a12, a23, a13];
%
Iden = [1 1 1 0 0 0]';
two3 = 2/3; stwo3=sqrt(two3); %constants
mu=mp(2); beta=mp(3); H=mp(4); Y0=mp(5); %material properties
ftol = Y0*1E-6; %tolerance for yield
stresstr = stressN + D*deps; %trial stress
I1 = sum(stresstr(1:3)); %trace(sigmatr)
str = stresstr - I1*Iden/3; %deviatoric stress
eta = str - alphaN; %shifted stress
etat = sqrt(eta(1)^2 + eta(2)^2 + eta(3)^2 ...
+ 2*(eta(4)^2 + eta(5)^2 + eta(6)^2));%norm of eta
fyld = etat - stwo3*(Y0+(1-beta)*H*epN); %trial yield function
if fyld < ftol %yield test
stress = stresstr; alpha = alphaN; ep = epN;%trial states are final
return;
else
gamma = fyld/(2*mu + two3*H); %plastic consistency param
ep = epN + gamma*stwo3; %updated eff. plastic strain
end
N = eta/etat; %unit vector normal to f
stress = stresstr - 2*mu*gamma*N; %updated stress
alpha = alphaN + two3*beta*H*gamma*N; %updated back stress
combHard1D.m
%
% 1D Linear combined isotropic/kinamtic hardening model
%
function [stress, alpha, ep]=combHard1D(mp,deps,stressN,alphaN,epN)
% Inputs:
% mp = [E, beta, H, Y0];
% deps = strain increment
% stressN = stress at load step N
% alphaN = back stress at load step N
% epN = plastic strain at load step N
%
E=mp(1); beta=mp(2); H=mp(3); Y0=mp(4); %material properties
ftol = Y0*1E-6; %tolerance for yield
stresstr = stressN + E*deps; %trial stress
etatr = stresstr - alphaN; %trial shifted stress
fyld = abs(etatr) - (Y0+(1-beta)*H*epN); %trial yield function
if fyld < ftol %yield test
stress = stresstr; alpha = alphaN; ep = epN;%trial states are final
return;
else
dep = fyld/(E+H); %plastic strain increment
end
stress = stresstr - sign(etatr)*E*dep; %updated stress
alpha = alphaN + sign(etatr)*beta*H*dep; %updated back stress
ep = epN + dep; %updated plastic strain
return
combHardTan.m
%
% Tangent stiffness for linear combined isotropic/kinamtic hardening model
%
function [Dtan]=combHardTan(mp,D,deps,stressN,alphaN,epN)
% Inputs:
% mp = [lambda, mu, beta, H, Y0];
% D = elastic stiffness matrix
% stressN = [s11, s22, s33, t12, t23, t13];
% alphaN = [a11, a22, a33, a12, a23, a13];
%
Iden = [1 1 1 0 0 0]';
two3 = 2/3; stwo3=sqrt(two3); %constants
mu=mp(2); beta=mp(3); H=mp(4); Y0=mp(5); %material properties
ftol = Y0*1E-6; %tolerance for yield
stresstr = stressN + D*deps; %trial stress
I1 = sum(stresstr(1:3)); %trace(sigmatr)
str = stresstr - I1*Iden/3; %deviatoric stress
eta = str - alphaN; %shifted stress
etat = sqrt(eta(1)^2 + eta(2)^2 + eta(3)^2 ...
+ 2*(eta(4)^2 + eta(5)^2 + eta(6)^2));%norm of eta
fyld = etat - stwo3*(Y0+(1-beta)*H*epN); %trial yield function
if fyld < ftol %yield test
Dtan = D; return; %elastic
end
gamma = fyld/(2*mu + two3*H); %plastic consistency param
N = eta/etat; %unit vector normal to f
var1 = 4*mu^2/(2*mu+two3*H);
var2 = 4*mu^2*gamma/etat; %coefficients
Dtan = D - (var1-var2)*N*N' + var2*Iden*Iden'/3;%tangent stiffness
Dtan(1,1) = Dtan(1,1) - var2; %contr. from 4th-order I
Dtan(2,2) = Dtan(2,2) - var2;
Dtan(3,3) = Dtan(3,3) - var2;
Dtan(4,4) = Dtan(4,4) - .5*var2;
Dtan(5,5) = Dtan(5,5) - .5*var2;
Dtan(6,6) = Dtan(6,6) - .5*var2;
E4_19.m
%
% Example 4.19 - Shear deformation of a square (finite rotation)
%
Young = 24000; nu=0.2; mu=Young/2/(1+nu); lambda=nu*Young/((1+nu)*(1-2*nu));
beta = 0; H = 1000; sY = 200*sqrt(3);
mp = [lambda mu beta H sY];
Iden=[1 1 1 0 0 0]';
D=2*mu*eye(6) + lambda*Iden*Iden';
D(4,4) = mu; D(5,5) = mu; D(6,6) = mu;
L = zeros(3,3);
stressN=[0 0 0 0 0 0]';
deps=[0 0 0 0 0 0]';
alphaN = [0 0 0 0 0 0]';
epN=0;
stressRN=stressN; alphaRN=alphaN;epRN=epN;
for i=1:15
deps(4) = 0.004; L(1,2) = 0.024; L(2,1) = -0.02;
[stressRN, alphaRN] = rotatedStress(L, stressRN, alphaRN);
[stressR, alphaR, epR]=combHard(mp,D,deps,stressRN,alphaRN,epRN);
[stress, alpha, ep]=combHard(mp,D,deps,stressN,alphaN,epN);
X(i) = i*deps(4); Y1(i) = stress(4); Y2(i) = stressR(4);
stressN = stress; alphaN = alpha; epN = ep;
stressRN = stressR; alphaRN = alphaR; epRN = epR;
end
X = [0 X]; Y1=[0 Y1]; Y2=[0 Y2]; plot(X,Y1,X,Y2);
E4_22.m
%
% Example 4.22 - Shear deformation of a square
%
clear
Young = 24000; nu=0.2; mu=Young/2/(1+nu); lambda=nu*Young/((1+nu)*(1-2*nu));
beta = 0; H = 1000; sY = 200*sqrt(3);
mp = [lambda mu beta H sY];
Iden=[1 1 1 0 0 0]';
D=2*mu*eye(6) + lambda*Iden*Iden';
D(4,4) = mu; D(5,5) = mu; D(6,6) = mu;
Iden=[1 1 1]';
%DM=2*mu*eye(3) + lambda*Iden*Iden';
L = zeros(3,3);
Dtan = zeros(6,6);
stressN=[0 0 0 0 0 0]';
deps=[0 0 0 0 0 0]';
alphaN = [0 0 0 0 0 0]';
epN=0;
stressRN=stressN; alphaRN=alphaN;epRN=epN;
bMN=[1 1 1 0 0 0]';
alphaMN = [0 0 0]';
epMN=0;
for i=1:15
DM=2*mu*eye(3) + lambda*Iden*Iden';
deps(4) = 0.004; L(1,2) = 0.024; L(2,1) = -0.02;
[stressRN, alphaRN] = rotatedStress(L, stressRN, alphaRN);
[stressR, alphaR, epR]=combHard(mp,D,deps,stressRN,alphaRN,epRN);
[stress, alpha, ep]=combHard(mp,D,deps,stressN,alphaN,epN);
[stressM, bM, alphaM, epM]=mulPlast(mp,DM,L,bMN,alphaMN,epMN);
[Dtan]=mulPlastTan(mp,DM,L,bMN,alphaMN,epMN);
X(i)=i*deps(4);Y1(i)=stress(4);Y2(i)=stressR(4);Y3(i)=stressM(4);
stressN = stress; alphaN = alpha; epN = ep;
stressRN = stressR; alphaRN = alphaR; epRN = epR;
bMN=bM; alphaMN = alphaM; epMN = epM;
% csvwrite('De.csv',Dtan)
end
Y3
X = [0 X]; Y1=[0 Y1]; Y2=[0 Y2]; Y3 = [0 Y3]; plot(X,Y1,X,Y2,X,Y3);
ELAST3D.m
function ELAST3D(ETAN, UPDATE, LTAN, NE, NDOF, XYZ, LE)
%***********************************************************************
% MAIN PROGRAM COMPUTING GLOBAL STIFFNESS MATRIX RESIDUAL FORCE FOR
% PLASTIC MATERIAL MODELS
%***********************************************************************
%%
global DISPTD FORCE GKF SIGMA
%
% Integration points and weights (2-point integration)
XG=[-0.57735026918963D0, 0.57735026918963D0];
WGT=[1.00000000000000D0, 1.00000000000000D0];
%
% Index for history variables (each integration pt)
INTN=0;
%
%LOOP OVER ELEMENTS, THIS IS MAIN LOOP TO COMPUTE K AND F
for IE=1:NE
% Nodal coordinates and incremental displacements
ELXY=XYZ(LE(IE,:),:);
% Local to global mapping
IDOF=zeros(1,24);
for I=1:8
II=(I-1)*NDOF+1;
IDOF(II:II+2)=(LE(IE,I)-1)*NDOF+1:(LE(IE,I)-1)*NDOF+3;
end
DSP=DISPTD(IDOF);
DSP=reshape(DSP,NDOF,8);
%
%LOOP OVER INTEGRATION POINTS
for LX=1:2, for LY=1:2, for LZ=1:2
E1=XG(LX); E2=XG(LY); E3=XG(LZ);
INTN = INTN + 1;
%
% Determinant and shape function derivatives
[~, SHPD, DET] = SHAPEL([E1 E2 E3], ELXY);
FAC=WGT(LX)*WGT(LY)*WGT(LZ)*DET;
%
% Strain
DEPS=DSP*SHPD';
DDEPS=[DEPS(1,1) DEPS(2,2) DEPS(3,3) ...
DEPS(1,2)+DEPS(2,1) DEPS(2,3)+DEPS(3,2) DEPS(1,3)+DEPS(3,1)]';
%
% Stress
STRESS = ETAN*DDEPS;
%
% Update stress
if UPDATE
SIGMA(:,INTN)=STRESS;
continue;
end
%
% Add residual force and tangent stiffness matrix
BM=zeros(6,24);
for I=1:8
COL=(I-1)*3+1:(I-1)*3+3;
BM(:,COL)=[SHPD(1,I) 0 0;
0 SHPD(2,I) 0;
0 0 SHPD(3,I);
SHPD(2,I) SHPD(1,I) 0;
0 SHPD(3,I) SHPD(2,I);
SHPD(3,I) 0 SHPD(1,I)];
end
%
% Residual forces
FORCE(IDOF) = FORCE(IDOF) - FAC*BM'*STRESS;
%
% Tangent stiffness
if LTAN
EKF = BM'*ETAN*BM;
GKF(IDOF,IDOF)=GKF(IDOF,IDOF)+FAC*EKF;
end
end, end, end
end
end
HYPER3D.m
function HYPER3D(MID, PROP, UPDATE, LTAN, NE, NDOF, XYZ, LE)
%***********************************************************************
% MAIN PROGRAM COMPUTING GLOBAL STIFFNESS MATRIX RESIDUAL FORCE FOR
% PLASTIC MATERIAL MODELS
%***********************************************************************
%%
global DISPTD FORCE GKF SIGMA
%
% Integration points and weights
XG=[-0.57735026918963D0, 0.57735026918963D0];
WGT=[1.00000000000000D0, 1.00000000000000D0];
%
% Index for history variables (each integration pt)
INTN=0;
%
%LOOP OVER ELEMENTS, THIS IS MAIN LOOP TO COMPUTE K AND F
for IE=1:NE
% Nodal coordinates and incremental displacements
ELXY=XYZ(LE(IE,:),:);
% Local to global mapping
IDOF=zeros(1,24);
for I=1:8
II=(I-1)*NDOF+1;
IDOF(II:II+2)=(LE(IE,I)-1)*NDOF+1:(LE(IE,I)-1)*NDOF+3;
end
DSP=DISPTD(IDOF);
DSP=reshape(DSP,NDOF,8);
%
%LOOP OVER INTEGRATION POINTS
for LX=1:2, for LY=1:2, for LZ=1:2
E1=XG(LX); E2=XG(LY); E3=XG(LZ);
INTN = INTN + 1;
%
% Determinant and shape function derivatives
[~, SHPD, DET] = SHAPEL([E1 E2 E3], ELXY);
FAC=WGT(LX)*WGT(LY)*WGT(LZ)*DET;
%
% Deformation gradient
F=DSP*SHPD' + eye(3);
%
% Computer stress and tangent stiffness
[STRESS DTAN] = Mooney(F, PROP(1), PROP(2), PROP(3), LTAN);
%
% Update plastic variables
if UPDATE
SIGMA(:,INTN)=STRESS;
continue;
end
%
% Add residual force and tangent stiffness matrix
BN=zeros(6,24);
BG=zeros(9,24);
for I=1:8
COL=(I-1)*3+1:(I-1)*3+3;
BN(:,COL)=[SHPD(1,I)*F(1,1) SHPD(1,I)*F(2,1) SHPD(1,I)*F(3,1);
SHPD(2,I)*F(1,2) SHPD(2,I)*F(2,2) SHPD(2,I)*F(3,2);
SHPD(3,I)*F(1,3) SHPD(3,I)*F(2,3) SHPD(3,I)*F(3,3);
SHPD(1,I)*F(1,2)+SHPD(2,I)*F(1,1) SHPD(1,I)*F(2,2)+SHPD(2,I)*F(2,1) SHPD(1,I)*F(3,2)+SHPD(2,I)*F(3,1);
SHPD(2,I)*F(1,3)+SHPD(3,I)*F(1,2) SHPD(2,I)*F(2,3)+SHPD(3,I)*F(2,2) SHPD(2,I)*F(3,3)+SHPD(3,I)*F(3,2);
SHPD(1,I)*F(1,3)+SHPD(3,I)*F(1,1) SHPD(1,I)*F(2,3)+SHPD(3,I)*F(2,1) SHPD(1,I)*F(3,3)+SHPD(3,I)*F(3,1)];
%
BG(:,COL)=[SHPD(1,I) 0 0;
SHPD(2,I) 0 0;
SHPD(3,I) 0 0;
0 SHPD(1,I) 0;
0 SHPD(2,I) 0;
0 SHPD(3,I) 0;
0 0 SHPD(1,I);
0 0 SHPD(2,I);
0 0 SHPD(3,I)];
end
%
% Residual forces
FORCE(IDOF) = FORCE(IDOF) - FAC*BN'*STRESS;
%
% Tangent stiffness
if LTAN
SIG=[STRESS(1) STRESS(4) STRESS(6);
STRESS(4) STRESS(2) STRESS(5);
STRESS(6) STRESS(5) STRESS(3)];
SHEAD=zeros(9);
SHEAD(1:3,1:3)=SIG;
SHEAD(4:6,4:6)=SIG;
SHEAD(7:9,7:9)=SIG;
%
EKF = BN'*DTAN*BN + BG'*SHEAD*BG;
GKF(IDOF,IDOF)=GKF(IDOF,IDOF)+FAC*EKF;
end
end; end; end;
end
end
Mooney.m
%
% Calculate 2nd PK stress and stiffness for Mooney-Rivlin material
%
function [Stress D] = Mooney(F, A10, A01, K, ltan)
% Inputs:
% F = Deformation gradient [3x3]
% A10, A01, K = Material constants
% ltan = 0 Calculate stress alone
% 1 Calculate stress and material stiffness
% Outputs:
% Stress = 2nd PK stress [S11, S22, S33, S12, S23, S13];
% D = Material stiffness [6x6]
%
X12 = 1/2; X13 = 1/3; X23 = 2/3; X43 = 4/3; X53 = 5/3; X89 = 8/9;
%
C = F'*F;
C1=C(1,1); C2=C(2,2); C3=C(3,3); C4=C(1,2); C5=C(2,3); C6=C(1,3);
I1 = C1+C2+C3;
I2 = C1*C2+C1*C3+C2*C3-C4^2-C5^2-C6^2;
I3 = det(C);
J3 = sqrt(I3);
J3M1 = J3 - 1.0D+00;
%
I1E = 2*[1 1 1 0 0 0]';
I2E = 2*[C2+C3, C3+C1, C1+C2, -C4, -C5, -C6]';
I3E = 2*[C2*C3-C5^2, C3*C1-C6^2, C1*C2-C4^2, ...
C5*C6-C3*C4, C6*C4-C1*C5, C4*C5-C2*C6]';
%
W1 = I3^(-X13); W2 = X13*I1*I3^(-X43); W3 = I3^(-X23);
W4 = X23*I2*I3^(-X53); W5 = X12*I3^(-X12);
%
J1E = W1*I1E - W2*I3E;
J2E = W3*I2E - W4*I3E;
J3E = W5*I3E;
%
Stress = A10*J1E + A01*J2E + K*J3M1*J3E;
%
% Material stiffness
%
D=zeros(6);
if ltan == 1
%
I2EE = [0 4 4 0 0 0; 4 0 4 0 0 0; 4 4 0 0 0 0;
0 0 0 -2 0 0; 0 0 0 0 -2 0; 0 0 0 0 0 -2];
I3EE = [ 0 4*C3 4*C2 0 -4*C5 0;
4*C3 0 4*C1 0 0 -4*C6;
4*C2 4*C1 0 -4*C4 0 0;
0 0 -4*C4 -2*C3 2*C6 2*C5;
-4*C5 0 0 2*C6 -2*C1 2*C4;
0 -4*C6 0 2*C5 2*C4 -2*C2];
%
W1 = X23*I3^(-X12); W2 = X89*I1*I3^(-X43); W3 = X13*I1*I3^(-X43);
W4 = X43*I3^(-X12); W5 = X89*I2*I3^(-X53); W6 = I3^(-X23);
W7 = X23*I2*I3^(-X53); W8 = I3^(-X12); W9 = X12*I3^(-X12);
%
J1EE = -W1*(J1E*J3E' + J3E*J1E') + W2*(J3E*J3E') - W3*I3EE;
J2EE = -W4*(J2E*J3E' + J3E*J2E') + W5*(J3E*J3E') + W6*I2EE - W7*I3EE;
J3EE = -W8*(J3E*J3E') + W9*I3EE;
%
D = A10*J1EE + A01*J2EE + K*(J3E*J3E') + K*J3M1*J3EE;
end
return;
mulPlast.m
%
% Multiplicative plasticity with linear combined hardening
%
function [stress, b, alpha, ep]=mulPlast(mp,D,L,b,alpha,ep)
%mp = [lambda, mu, beta, H, Y0];
%D = elasticity matrix b/w prin stress & log prin stretch (3x3)
%L = [dui/dxj] velocity gradient
%b = elastic left C-G deformation vector (6x1)
%alpha = principal back stress (3x1)
%ep = effective plastic strain
%
EPS=1E-12;
Iden = [1 1 1]'; two3 = 2/3; stwo3=sqrt(two3); %constants
mu=mp(2); beta=mp(3); H=mp(4); Y0=mp(5); %material properties
ftol = Y0*1E-6; %tolerance for yield
R = inv(eye(3)-L); %inc. deformation gradient
bm=[b(1) b(4) b(6);b(4) b(2) b(5);b(6) b(5) b(3)];
bm = R*bm*R'; %trial elastic left C-G
b=[bm(1,1) bm(2,2) bm(3,3) bm(1,2) bm(2,3) bm(1,3)]';
[~,P]=eig(bm); %eigenvalues
eigen=sort(real([P(1,1) P(2,2) P(3,3)]))'; %principal stretch
%
% Duplicated eigenvalues
TMP=-1;
for I=1:2
if abs(eigen(1)-eigen(3)) < EPS
eigen(I)=eigen(I)+TMP*EPS;
TMP=-TMP;
end
end
if abs(eigen(1)-eigen(2)) < EPS; eigen(2) = eigen(2) + EPS; end;
if abs(eigen(2)-eigen(3)) < EPS; eigen(2) = eigen(2) + EPS; end;
%
% EIGENVECTOR MATRIX N*N' = M(6,*)
M=zeros(6,3); %eigenvector matrices
for K=1:3
KB=1+mod(K,3);
KC=1+mod(KB,3);
EA=eigen(K);
EB=eigen(KB);
EC=eigen(KC);
D1=EB-EA;
D2=EC-EA;
DA = 1 / (D1 * D2);
M(1,K)=((b(1)-EB)*(b(1)-EC)+b(4)*b(4)+b(6)*b(6))*DA;
M(2,K)=((b(2)-EB)*(b(2)-EC)+b(4)*b(4)+b(5)*b(5))*DA;
M(3,K)=((b(3)-EB)*(b(3)-EC)+b(5)*b(5)+b(6)*b(6))*DA;
M(4,K)=(b(4)*(b(1)-EB+b(2)-EC)+b(5)*b(6))*DA;
M(5,K)=(b(5)*(b(2)-EB+b(3)-EC)+b(4)*b(6))*DA;
M(6,K)=(b(6)*(b(3)-EB+b(1)-EC)+b(4)*b(5))*DA;
end
%
eigen=sort(real([P(1,1) P(2,2) P(3,3)]))'; %principal stretch
deps = 0.5*log(eigen); %logarithmic
sigtr = D*deps; %trial principal stress
eta = sigtr - alpha - sum(sigtr)*Iden/3; %shifted stress
etat = norm(eta); %norm of eta
fyld = etat - stwo3*(Y0+(1-beta)*H*ep); %trial yield function
if fyld < ftol %yield test
sig = sigtr; %trial states are final
stress = M*sig; %stress (6x1)
else
gamma = fyld/(2*mu + two3*H); %plastic consistency param
ep = ep + gamma*stwo3; %updated eff. plastic strain
N = eta/etat; %unit vector normal to f
deps = deps - gamma*N; %updated elastic strain
sig = sigtr - 2*mu*gamma*N; %updated stress
alpha = alpha + two3*beta*H*gamma*N; %updated back stress
stress = M*sig; %stress (6x1)
b = M*exp(2*deps); %updated elastic left C-G
end
mulPlastTan.m
%
% Tangent stiffness of multiplicative plasticity with linear hardening
%
function [Dtan] = mulPlastTan(mp, D, L, b, alpha, ep)
%
EPS = 1E-12;
Iden = [1 1 1]'; two3 = 2/3; stwo3 = sqrt(two3); %constants
mu = mp(2); beta = mp(3); H = mp(4); Y0 = mp(5); %material properties
ftol = Y0 * 1E-6; %tolerance for yield
R = inv(eye(3) - L); %inc. deformation gradient
bm = [b(1) b(4) b(6); b(4) b(2) b(5); b(6) b(5) b(3)];
bm = R * bm * R'; %trial elastic left C-G
b = [bm(1, 1) bm(2, 2) bm(3, 3) bm(1, 2) bm(2, 3) bm(1, 3)]';
[~, P] = eig(bm); %eigenvalues
eigen = sort(real([P(1, 1) P(2, 2) P(3, 3)]))'; %principal stretch
%
% Duplicated eigenvalues
TMP = -1;
for I = 1:2
if abs(eigen(1) - eigen(3)) < EPS
eigen(I) = eigen(I) + TMP * EPS;
TMP = -TMP;
end
end
if abs(eigen(1) - eigen(2)) < EPS; eigen(2) = eigen(2) + EPS; end;
if abs(eigen(2) - eigen(3)) < EPS; eigen(2) = eigen(2) + EPS; end;
%
% EIGENVECTOR MATRIX N*N = M(6,*)
M = zeros(6, 3); %eigenvector matrices
for K = 1:3
KB = 1 + mod(K, 3);
KC = 1 + mod(KB, 3);
EA = eigen(K);
EB = eigen(KB);
EC = eigen(KC);
D1 = EB - EA;
D2 = EC - EA;
DA = 1 / (D1 * D2);
M(1, K) = ((b(1) - EB) * (b(1) - EC) + b(4) * b(4) + b(6) * b(6)) * DA;
M(2, K) = ((b(2) - EB) * (b(2) - EC) + b(4) * b(4) + b(5) * b(5)) * DA;
M(3, K) = ((b(3) - EB) * (b(3) - EC) + b(5) * b(5) + b(6) * b(6)) * DA;
M(4, K) = (b(4) * (b(1) - EB + b(2) - EC) + b(5) * b(6)) * DA;
M(5, K) = (b(5) * (b(2) - EB + b(3) - EC) + b(4) * b(6)) * DA;
M(6, K) = (b(6) * (b(3) - EB + b(1) - EC) + b(4) * b(5)) * DA;
end
%
deps = 0.5 * log(eigen); %logarithmic
sigtr = D * deps; %trial principal stress
sig = sigtr;
eta = sigtr - alpha - sum(sigtr) * Iden / 3; %shifted stress
etat = norm(eta); %norm of eta
fyld = etat - stwo3 * (Y0 + (1 - beta) * H * ep); %trial yield function
if fyld >= ftol %yield test
gamma = fyld / (2 * mu + two3 * H); %plastic consistency param
N = eta / etat; %unit vector normal to f
sig = sigtr - 2 * mu * gamma * N; %updated stress
var1 = 4 * mu^2 / (2 * mu + two3 * H);
var2 = 4 * mu^2 * gamma / etat; %coefficients
D = D - (var1 - var2) * (N * N') + var2 * (Iden * Iden') / 3; %tangent stiffness
D(1, 1) = D(1, 1) - var2; %contr. from 4th-order I
D(2, 2) = D(2, 2) - var2;
D(3, 3) = D(3, 3) - var2;
end
J1 = sum(eigen);
J3 = eigen(1) * eigen(2) * eigen(3);
I2 = [1 1 1 0 0 0]';
I4 = eye(6); I4(4, 4) = .5; I4(5, 5) = .5; I4(6, 6) = .5;
Ibb = [0, b(4)^2 - b(1) * b(2), b(6)^2 - b(1) * b(3), 0, b(4) * b(6) - b(1) * b(5), 0;
b(4) * b(4) - b(1) * b(2), 0, b(5)^2 - b(2) * b(3), 0, 0, b(4) * b(5) - b(2) * b(6);
b(6)^2 - b(1) * b(3), b(5)^2 - b(2) * b(3), 0, b(5) * b(6) - b(3) * b(4), 0, 0;
0, 0, b(5) * b(6) - b(3) * b(4), (b(1) * b(2) - b(4)^2) / 2, (b(2) * b(6) - b(4) * b(5)) / 2, (b(1) * b(5) - b(4) * b(6)) / 2;
b(4) * b(6) - b(1) * b(5), 0, 0, (b(2) * b(6) - b(4) * b(5)) / 2, (b(2) * b(3) - b(5)^2) / 2, (b(3) * b(4) - b(5) * b(6)) / 2;
0, b(4) * b(5) - b(2) * b(6), 0, (b(1) * b(5) - b(4) * b(6)) / 2, (b(3) * b(4) - b(5) * b(6)) / 2, (b(1) * b(3) - b(6)^2) / 2];
%
d1 = 1 / ((eigen(2) - eigen(1)) * (eigen(3) - eigen(1)));
d2 = 1 / ((eigen(3) - eigen(2)) * (eigen(1) - eigen(2)));
d3 = 1 / ((eigen(1) - eigen(3)) * (eigen(2) - eigen(3)));
t11 = -J3 * d1 / eigen(1); t12 = -J3 * d2 / eigen(2); t13 = -J3 * d3 / eigen(3);
t21 = d1 * eigen(1); t22 = d2 * eigen(2); t23 = d3 * eigen(3);
t31 = t21 * (J1 - 4 * eigen(1)); t32 = t22 * (J1 - 4 * eigen(2)); t33 = t23 * (J1 - 4 * eigen(3));
%
CT1 = d1 * Ibb + t11 * (I4 - (I2 - M(:, 1)) * (I2 - M(:, 1))') + t21 * (b * M(:, 1)' + M(:, 1) * b') + t31 * M(:, 1) * M(:, 1)';
CT2 = d2 * Ibb + t12 * (I4 - (I2 - M(:, 2)) * (I2 - M(:, 2))') + t22 * (b * M(:, 2)' + M(:, 2) * b') + t32 * M(:, 2) * M(:, 2)';
CT3 = d3 * Ibb + t13 * (I4 - (I2 - M(:, 3)) * (I2 - M(:, 3))') + t23 * (b * M(:, 3)' + M(:, 3) * b') + t33 * M(:, 3) * M(:, 3)';
%
Dtan = M * D * M' + 2 * (CT1 * sig(1) + CT2 * sig(2) + CT3 * sig(3));
NLFEA.m
function NLFEA(ITRA,TOL,ATOL,NTOL,TIMS,NOUT,MID,PROP,EXTFORCE,SDISPT,XYZ,LE)
%***********************************************************************
% MAIN PROGRAM FOR HYPERELASTIC/ELASTOPLASTIC ANALYSIS
%***********************************************************************
%%
global DISPDD DISPTD FORCE GKF %Global variables
%
[NUMNP, NDOF] = size(XYZ); % Analysis parameters
NE = size(LE,1);
NEQ = NDOF*NUMNP;
%
DISPTD=zeros(NEQ,1); DISPDD=zeros(NEQ,1); % Nodal displacement & increment
if MID >= 0, ETAN=PLSET(PROP, MID, NE); end % Initialize material properties
%
ITGZONE(XYZ, LE, NOUT); % Check element connectivity
%
% Load increments [Start End Increment InitialLoad FinalLoad]
NLOAD=size(TIMS,2);
ILOAD=1; % First load increment
TIMEF=TIMS(1,ILOAD); % Starting time
TIMEI=TIMS(2,ILOAD); % Ending time
DELTA=TIMS(3,ILOAD); % Time increment
CUR1=TIMS(4,ILOAD); % Starting load factor
CUR2=TIMS(5,ILOAD); % Ending load factor
DELTA0 = DELTA; % Saved time increment
TIME = TIMEF; % Starting time
TDELTA = TIMEI - TIMEF; % Time interval for load step
ITOL = 1; % Bisection level
TARY=zeros(NTOL,1); % Time stamps for bisections
%
% Load increment loop
%----------------------------------------------------------------------
ISTEP = -1; FLAG10 = 1;
while(FLAG10 == 1) % Solution has been converged
FLAG10 = 0; FLAG11 = 1; FLAG20 = 1;
%
CDISP = DISPTD; % Store converged displacement
%
if(ITOL==1) % No bisection
DELTA = DELTA0;
TARY(ITOL) = TIME + DELTA;
else % Recover previous bisection
ITOL = ITOL-1; % Reduce the bisection level
DELTA = TARY(ITOL)-TARY(ITOL+1); % New time increment
TARY(ITOL+1) = 0; % Empty converged bisection level
ISTEP = ISTEP - 1; % Decrease load increment
end
TIME0 = TIME; % Save the current time
%
% Update stresses and history variables
UPDATE=true; LTAN=false;
if MID ==0, ELAST3D(ETAN, UPDATE, LTAN, NE, NDOF, XYZ, LE);
elseif MID > 0, PLAST3D(MID, PROP, ETAN, UPDATE, LTAN, NE, NDOF, XYZ, LE);
elseif MID < 0, HYPER3D(PROP, UPDATE, LTAN, NE, NDOF, XYZ, LE);
else fprintf(NOUT,'\t\t *** Wrong material ID ***\n'); return;
end
%
% Print results
if(ISTEP>=0), PROUT(NOUT, TIME, NUMNP, NE, NDOF); end
%
TIME = TIME + DELTA; % Increase time
ISTEP = ISTEP + 1;
%
% Check time and control bisection
while(FLAG11 == 1) % Bisection loop start
FLAG11 = 0;
if ((TIME-TIMEI)>1E-10) % Time passed the end time
if ((TIMEI+DELTA-TIME)>1E-10) % One more at the end time
DELTA=TIMEI+DELTA-TIME; % Time increment to the end
DELTA0=DELTA; % Saved time increment
TIME=TIMEI; % Current time is the end
else
ILOAD=ILOAD+1; % Progress to next load step
if(ILOAD>NLOAD) % Finished final load step
FLAG10 = 0; % Stop the program
break;
else % Next load step
TIME=TIME-DELTA;
DELTA=TIMS(3,ILOAD);
DELTA0=DELTA;
TIME = TIME + DELTA;
TIMEF = TIMS(1,ILOAD);
TIMEI = TIMS(2,ILOAD);
TDELTA = TIMEI - TIMEF;
CUR1 = TIMS(4,ILOAD);
CUR2 = TIMS(5,ILOAD);
end
end
end
%
% Load factor and prescribed displacements
FACTOR = CUR1 + (TIME-TIMEF)/TDELTA*(CUR2-CUR1);
SDISP = DELTA*SDISPT(:,3)/TDELTA*(CUR2-CUR1);
%
% Start convergence iteration
%------------------------------------------------------------------
ITER = 0;
DISPDD = zeros(NEQ,1);
while(FLAG20 == 1)
FLAG20 = 0;
ITER = ITER + 1;
%
% Initialize global stiffness K and residual vector F
GKF = sparse(NEQ,NEQ);
FORCE = sparse(NEQ,1);
%
% Assemble K and F
UPDATE=false; LTAN=true;
if MID ==0, ELAST3D(ETAN, UPDATE, LTAN, NE, NDOF, XYZ, LE);
elseif MID > 0, PLAST3D(MID, PROP, ETAN, UPDATE, LTAN, NE, NDOF, XYZ, LE);
elseif MID < 0, HYPER3D(PROP, UPDATE, LTAN, NE, NDOF, XYZ, LE);
end
%
% Increase external force
if size(EXTFORCE,1)>0
LOC = NDOF*(EXTFORCE(:,1)-1)+EXTFORCE(:,2);
FORCE(LOC) = FORCE(LOC) + FACTOR*EXTFORCE(:,3);
end
%
% Prescribed displacement BC
NDISP=size(SDISPT,1);
if NDISP~=0
FIXEDDOF=NDOF*(SDISPT(:,1)-1)+SDISPT(:,2);
GKF(FIXEDDOF,:)=zeros(NDISP,NEQ);
GKF(FIXEDDOF,FIXEDDOF)=PROP(1)*eye(NDISP);
%
FORCE(FIXEDDOF)=0;
if ITER==1, FORCE(FIXEDDOF) = PROP(1)*SDISP(:); end
end
% Check convergence
if(ITER>1)
FIXEDDOF=NDOF*(SDISPT(:,1)-1)+SDISPT(:,2);
ALLDOF=1:NEQ;
FREEDOF=setdiff(ALLDOF,FIXEDDOF);
RESN=max(abs(FORCE(FREEDOF)));
OUTPUT(1, ITER, RESN, TIME, DELTA)
%
if(RESN<TOL)
FLAG10 = 1;
break;
end
%
if ((RESN>ATOL)||(ITER>=ITRA)) % Start bisection
ITOL = ITOL + 1;
if(ITOL<=NTOL)
DELTA = 0.5*DELTA;
TIME = TIME0 + DELTA;
TARY(ITOL) = TIME;
DISPTD=CDISP;
fprintf(1,'Not converged. Bisecting load increment %3d\n',ITOL);
else
fprintf(1,'\t\t *** Max No. of bisection ***\n'); return;
end
FLAG11 = 1;
FLAG20 = 1;
break;
end
end
%
% Solve the system equation
if(FLAG11 == 0)
SOLN = GKF\FORCE;
DISPDD = DISPDD + SOLN;
DISPTD = DISPTD + SOLN;
FLAG20 = 1;
else
FLAG20 = 0;
end
if(FLAG10 == 1), break; end
end %20 Convergence iteration
end %11 Bisection
end %10 Load increment
%
% Successful end of program
fprintf(1,'\t\t *** Successful end of program ***\n');
fprintf(NOUT,'\t\t *** Successful end of program ***\n');
end
function OUTPUT(FLG, ITER, RESN, TIME, DELTA)
%*************************************************************************
% Print convergence iteration history
%*************************************************************************
%%
if FLG == 1
if ITER>2
fprintf(1,'%27d %14.5e \n',ITER,full(RESN));
else
fprintf(1,'\n \t Time Time step Iter \t Residual \n');
fprintf(1,'%10.5f %10.3e %5d %14.5e \n',TIME,DELTA,ITER,full(RESN));
end
end
end
function PROUT(NOUT, TIME, NUMNP, NE, NDOF)
%*************************************************************************
% Print converged displacements and stresses
%*************************************************************************
%%
global SIGMA DISPTD
%
fprintf(NOUT,'\r\n\r\nTIME = %11.3e\r\n\r\nNodal Displacements\r\n',TIME);
fprintf(NOUT,'\r\n Node U1 U2 U3');
for I=1:NUMNP
II=NDOF*(I-1);
fprintf(NOUT,'\r\n%5d %11.3e %11.3e %11.3e',I,DISPTD(II+1:II+3));
end
fprintf(NOUT,'\r\n\r\nElement Stress\r\n');
fprintf(NOUT,'\r\n S11 S22 S33 S12 S23 S13');
for I=1:NE
fprintf(NOUT,'\r\nElement %5d',I);
II=(I-1)*8;
fprintf(NOUT,'\r\n%11.3e %11.3e %11.3e %11.3e %11.3e %11.3e',SIGMA(1:6,II+1:II+8));
end
fprintf(NOUT,'\r\n\r\n');
end
function ETAN=PLSET(PROP, MID, NE)
%**********************************************************************
% Initialize history variables and elastic stiffness matrix
% XQ : 1-6 = Back stress alpha, 7 = Effective plastic strain
% SIGMA : Stress for rate-form plasticity
% : Left Cauchy-Green tensor XB for multiplicative plasticity
% ETAN : Elastic stiffness matrix
%**********************************************************************
%%
global SIGMA XQ
%
LAM=PROP(1);
MU=PROP(2);
%
N = 8*NE;
%
if MID > 30
SIGMA=zeros(12,N);
XQ=zeros(4,N);
SIGMA(7:9,:)=1;
ETAN=[LAM+2*MU LAM LAM ;
LAM LAM+2*MU LAM ;
LAM LAM LAM+2*MU];
else
SIGMA=zeros(6,N);
XQ=zeros(7,N);
ETAN=[LAM+2*MU LAM LAM 0 0 0;
LAM LAM+2*MU LAM 0 0 0;
LAM LAM LAM+2*MU 0 0 0;
0 0 0 MU 0 0;
0 0 0 0 MU 0;
0 0 0 0 0 MU];
end
end
function VOLUME = ITGZONE(XYZ, LE, NOUT)
%*************************************************************************
% Check element connectivity and calculate volume
%*************************************************************************
%%
EPS=1E-7;
NE = size(LE,1);
VOLUME=0;
for I=1:NE
ELXY=XYZ(LE(I,:),:);
[~, ~, DET] = SHAPEL([0 0 0], ELXY);
DVOL = 8*DET;
if DVOL < EPS
fprintf(NOUT,'\n??? Negative Jacobian ???\nElement connectivity\n');
fprintf(NOUT,'%5d',LE(I,:));
fprintf(NOUT,'\nNodal Coordinates\n');
fprintf(NOUT,'%10.3e %10.3e %10.3e\n',ELXY');
error('Negative Jacobian');
end
VOLUME = VOLUME + DVOL;
end
end
ONE_ELEMENT.m
%
% One element example
%
% Nodal coordinates
XYZ=[0 0 0;1 0 0;1 1 0;0 1 0;0 0 1;1 0 1;1 1 1;0 1 1]*0.01;
%
% Element connectivity
LE=[1 2 3 4 5 6 7 8];
%
% External forces [Node, DOF, Value]
EXTFORCE=[5 3 10.0E3; 6 3 10.0E3; 7 3 10.0E3; 8 3 10.0E3];
%EXTFORCE=[];
%
% Prescribed displacements [Node, DOF, Value]
SDISPT=[1 1 0;1 2 0;1 3 0;2 2 0;2 3 0;3 3 0;4 1 0;4 3 0];
%SDISPT=[1 1 0;1 2 0;1 3 0;2 2 0;2 3 0;3 3 0;4 1 0;4 3 0;5 3 0.5;6 3 0.5;7 3 0.5;8 3 0.5];
%
% Load increments [Start End Increment InitialFactor FinalFactor]
TIMS=[0.0 0.5 0.1 0.0 0.5; 0.5 1.0 0.1 0.5 1.0]';
%TIMS=[0.0 0.8 0.4 0.0 0.8; 0.8 1.1 0.1 0.8 1.1]';
%TIMS=[0.0 1.0 1.0 0.0 1.0]';
%
% Contact elements
LEC=[];
%
% Material properties
% MID:0(Linear elastic)
% PROP=[LAMBDA NU]
%MID=0;
%PROP=[1.1538E6, 7.6923E5];
% 1(Linear hardening) 2(Saturated hardening)
% 31(Linear hardening finite deformation)
% 32(Saturated hardening finite deformation)
% PROP=[LAMDA MU BETA H Y0 CE1 FTOL]
%MID=1;
%PROP=[110.747E3 80.1938E3 0.0 10.E3 50.0E3 1.0 0.0001];
% -N(Hyperelasticity with N parameters)
% PROP=[A10 A01 D]
MID=31;
PROP=[110.747E9 80.1938E9 0.0 1.E9 4.0E8];
%MID=-2;
%PROP=[80 20 10000];
%
% Set program parameters
ITRA=70; ATOL=1.0E5; NTOL=6; TOL=1E-6;
%clc
%
% Calling main function
NOUT = fopen('output.txt','w');
NLFEA(ITRA, TOL, ATOL, NTOL, TIMS, NOUT, MID, PROP, EXTFORCE, SDISPT, XYZ, LE);
fclose(NOUT);
PLAST3D.m
function PLAST3D(MID, PROP, ETAN, UPDATE, LTAN, NE, NDOF, XYZ, LE)
%***********************************************************************
% MAIN PROGRAM COMPUTING GLOBAL STIFFNESS MATRIX RESIDUAL FORCE FOR
% PLASTIC MATERIAL MODELS
%***********************************************************************
%%
global DISPDD DISPTD FORCE GKF SIGMA XQ
%
% Integration points and weights (2-point integration)
XG=[-0.57735026918963D0, 0.57735026918963D0];
WGT=[1.00000000000000D0, 1.00000000000000D0];
%
% Index for history variables (each integration pt)
INTN=0;
%
%LOOP OVER ELEMENTS, THIS IS MAIN LOOP TO COMPUTE K AND F
for IE=1:NE
% Nodal coordinates and incremental displacements
ELXY=XYZ(LE(IE,:),:);
% Local to global mapping
IDOF=zeros(1,24);
for I=1:8
II=(I-1)*NDOF+1;
IDOF(II:II+2)=(LE(IE,I)-1)*NDOF+1:(LE(IE,I)-1)*NDOF+3;
end
DSP=DISPTD(IDOF);
DSPD=DISPDD(IDOF);
DSP=reshape(DSP,NDOF,8);
DSPD=reshape(DSPD,NDOF,8);
%
%LOOP OVER INTEGRATION POINTS
for LX=1:2, for LY=1:2, for LZ=1:2
E1=XG(LX); E2=XG(LY); E3=XG(LZ);
INTN = INTN + 1;
%
% Determinant and shape function derivatives
[~, SHPD, DET] = SHAPEL([E1 E2 E3], ELXY);
FAC=WGT(LX)*WGT(LY)*WGT(LZ)*DET;
%
% Previous converged history variables
if MID > 30
NALPHA=3;
STRESSN=SIGMA(7:12,INTN);
else
NALPHA=6;
STRESSN=SIGMA(1:6,INTN);
end
ALPHAN=XQ(1:NALPHA,INTN);
EPN=XQ(NALPHA+1,INTN);
%
% Strain increment
if MID == 2 || MID == 31
F=DSP*SHPD' + eye(3);
SHPD=inv(F)'*SHPD;
end
DEPS=DSPD*SHPD';
DDEPS=[DEPS(1,1) DEPS(2,2) DEPS(3,3) ...
DEPS(1,2)+DEPS(2,1) DEPS(2,3)+DEPS(3,2) DEPS(1,3)+DEPS(3,1)]';
%
% Computer stress, back stress & effective plastic strain
if MID == 1
% Infinitesimal plasticity
[STRESS, ALPHA, EP]=combHard(PROP,ETAN,DDEPS,STRESSN,ALPHAN,EPN);
elseif MID == 2
% Plasticity with finite rotation
FAC=FAC*det(F);
[STRESSN, ALPHAN] = rotatedStress(DEPS, STRESSN, ALPHAN);
[STRESS, ALPHA, EP]=combHard(PROP,ETAN,DDEPS,STRESSN,ALPHAN,EPN);
elseif MID == 31
[STRESS, B, ALPHA, EP]=mulPlast(PROP,ETAN,DEPS,STRESSN,ALPHAN,EPN);
end
%
% Update plastic variables
if UPDATE
SIGMA(1:6,INTN)=STRESS;
XQ(:,INTN)= [ALPHA; EP];
if MID > 30
SIGMA(7:12,INTN)=B;
end
continue;
end
%
% Add residual force and tangent stiffness matrix
BM=zeros(6,24);
BG=zeros(9,24);
for I=1:8
COL=(I-1)*3+1:(I-1)*3+3;
BM(:,COL)=[SHPD(1,I) 0 0;
0 SHPD(2,I) 0;
0 0 SHPD(3,I);
SHPD(2,I) SHPD(1,I) 0;
0 SHPD(3,I) SHPD(2,I);
SHPD(3,I) 0 SHPD(1,I)];
%
BG(:,COL)=[SHPD(1,I) 0 0;
SHPD(2,I) 0 0;
SHPD(3,I) 0 0;
0 SHPD(1,I) 0;
0 SHPD(2,I) 0;
0 SHPD(3,I) 0;
0 0 SHPD(1,I);
0 0 SHPD(2,I);
0 0 SHPD(3,I)];
end
%
% Residual forces
FORCE(IDOF) = FORCE(IDOF) - FAC*BM'*STRESS;
%
% Tangent stiffness
if LTAN
if MID == 1
DTAN=combHardTan(PROP,ETAN,DDEPS,STRESSN,ALPHAN,EPN);
EKF = BM'*DTAN*BM;
elseif MID == 2
DTAN=combHardTan(PROP,ETAN,DDEPS,STRESSN,ALPHAN,EPN);
CTAN=[-STRESS(1) STRESS(1) STRESS(1) -STRESS(4) 0 -STRESS(6);
STRESS(2) -STRESS(2) STRESS(2) -STRESS(4) -STRESS(5) 0;
STRESS(3) STRESS(3) -STRESS(3) 0 -STRESS(5) -STRESS(6);
-STRESS(4) -STRESS(4) 0 -0.5*(STRESS(1)+STRESS(2)) -0.5*STRESS(6) -0.5*STRESS(5);
0 -STRESS(5) -STRESS(5) -0.5*STRESS(6) -0.5*(STRESS(2)+STRESS(3)) -0.5*STRESS(4);
-STRESS(6) 0 -STRESS(6) -0.5*STRESS(5) -0.5*STRESS(4) -0.5*(STRESS(1)+STRESS(3))];
SIG=[STRESS(1) STRESS(4) STRESS(6);
STRESS(4) STRESS(2) STRESS(5);
STRESS(6) STRESS(5) STRESS(3)];
SHEAD=kron(eye(3),SIG);
EKF = BM'*(DTAN+CTAN)*BM + BG'*SHEAD*BG;
elseif MID == 31
DTAN=mulPlastTan(PROP,ETAN,DEPS,STRESSN,ALPHAN,EPN);
SIG=[STRESS(1) STRESS(4) STRESS(6);
STRESS(4) STRESS(2) STRESS(5);
STRESS(6) STRESS(5) STRESS(3)];
SHEAD=kron(eye(3),SIG);
%
EKF = BM'*DTAN*BM + BG'*SHEAD*BG;
end
GKF(IDOF,IDOF)=GKF(IDOF,IDOF)+FAC*EKF;
end
end, end, end
end
end
rotatedStress.m
%
% Rotate stress and back stress to the mid-configuration
%
function [stress, alpha] = rotatedStress(L, S, A)
%L = [dui/dxj]
str=[S(1) S(4) S(6);S(4) S(2) S(5);S(6) S(5) S(3)];
alp=[A(1) A(4) A(6);A(4) A(2) A(5);A(6) A(5) A(3)];
R = L*inv(eye(3) - 0.5*L);
W = .5*(R-R');
R = eye(3) + inv(eye(3) - 0.5*W)*W;
str = R*str*R';
alp = R*alp*R';
stress=[str(1,1) str(2,2) str(3,3) str(1,2) str(2,3) str(1,3)]';
alpha =[alp(1,1) alp(2,2) alp(3,3) alp(1,2) alp(2,3) alp(1,3)]';
return;
SHAPEL.m
function [SF, GDSF, DET] = SHAPEL(XI, ELXY)
%*************************************************************************
% Compute shape function, derivatives, and determinant of hexahedron element
%*************************************************************************
%%
XNODE=[-1 1 1 -1 -1 1 1 -1;
-1 -1 1 1 -1 -1 1 1;
-1 -1 -1 -1 1 1 1 1];
QUAR = 0.125;
SF=zeros(8,1);
DSF=zeros(3,8);
for I=1:8
XP = XNODE(1,I);
YP = XNODE(2,I);
ZP = XNODE(3,I);
%
XI0 = [1+XI(1)*XP 1+XI(2)*YP 1+XI(3)*ZP];
%
SF(I) = QUAR*XI0(1)*XI0(2)*XI0(3);
DSF(1,I) = QUAR*XP*XI0(2)*XI0(3);
DSF(2,I) = QUAR*YP*XI0(1)*XI0(3);
DSF(3,I) = QUAR*ZP*XI0(1)*XI0(2);
end
GJ = DSF*ELXY;
DET = det(GJ);
GJINV=inv(GJ);
GDSF=GJINV*DSF;
end
TWO_ELEMENT.m
%
% Two-element example page 131
%
clear
clc
% Nodal coordinates
XYZ=[0 0 0; 1 0 0; 1 1 0; 0 1 0;
0 0 1; 1 0 1; 1 1 1; 0 1 1;
0 0 2; 1 0 2; 1 1 2; 0 1 2]*0.01;
%
% Element connectivity
LE=[1 2 3 4 5 6 7 8;
5 6 7 8 9 10 11 12];
%
% External forces [Node, DOF, Value]
EXTFORCE=[9 3 10.0E3; 10 3 10.0E3; 11 3 10.0E3; 12 3 10.0E3];
%
% Prescribed displacements [Node, DOF, Value]
SDISPT=[1 1 0;1 2 0;1 3 0;2 2 0;2 3 0;3 3 0;4 1 0;4 3 0];
%
% Load increments [Start End Increment InitialFactor FinalFactor]
TIMS=[0.0 0.8 0.4 0.0 0.8; 0.8 1.1 0.1 0.8 1.1]';
%
% Material properties PROP=[LAMDA MU BETA H Y0]
MID=31;
PROP=[110.747E9 80.1938E9 0.0 1.E9 4.0E8];
%
% Set program parameters
ITRA=70; ATOL=1.0E5; NTOL=6; TOL=1E-6;
%
% Calling main function
NOUT = fopen('output.txt','w');
NLFEA(ITRA, TOL, ATOL, NTOL, TIMS, NOUT, MID, PROP, EXTFORCE, SDISPT, XYZ, LE);
fclose(NOUT);
浙公网安备 33010602011771号