matlab codes for GINCA (geometry-inspired numerical convex analysis) - Home / 主页 - 娄 燕山
Due to my fault that the matlab codes were not properly published with the article <Y Lou, C Zhang, P Wu, JW Yoon, 2024. New geometry-inspired numerical convex analysis method for yield functions under isotropic and anisotropic hardenings. International Journal of Solids and Structures, 112582>. The matlab codes are shared below:
Appendix A: GINCA code to compute the convex domain of the one-parameter yield function (Drucker, Cazacu-Barlat2004 and Cazacu2018) based on Fig. 1
function [] = myfunction( )
clc
clear
tol=1.0e-6;
angled=180:-180:0; %angled = 0 for upper limit of c, while angled = 180 for lower limit of c
angle=angled*pi/180;
n=length(angled);
d_ifa=0.1;
delta_ifa=0.0001;
for i=1:n
radius_u=6;
radius_l=0;
radius_m=(radius_u+radius_l)/2;
dr=radius_u-radius_l;
while dr>tol
c=radius_m*cos(angle(i));
param=[c];
convex=1; %1 is for convex and 0 for concave
for l=0:d_ifa:360
angld_A=l/180*pi;
s11_A_try=cos(angld_A);
s22_A_try=sin(angld_A);
f_A=Yldfunction(s11_A_try,s22_A_try,param);
s11_A=s11_A_try/f_A;
s22_A=s22_A_try/f_A;
angld_B=l/180*pi+delta_ifa;
s11_B_try=cos(angld_B);
s22_B_try=sin(angld_B);
f_B=Yldfunction(s11_B_try,s22_B_try,param);
s11_B=s11_B_try/f_B;
s22_B=s22_B_try/f_B;
s11_C=(s11_A+s11_B)/2;
s22_C=(s22_A+s22_B)/2;
f_C=Yldfunction(s11_C,s22_C,param);
if f_C>1
convex=0;
break
end
end
if convex==1
radius_l=radius_m;
else
radius_u=radius_m;
end
radius_m=(radius_u+radius_l)/2;
dr=(radius_u-radius_l)/2;
end
c_limit(i)=radius_m*cos(angle(i))
end
end
function y=Yldfunction(s11,s22,param)
c=param(1);
I1=s11+s22;
s1=s11-I1/3;
s2=s22-I1/3;
s3=-I1/3;
J2=-s1*s2-s2*s3-s1*s3;
J3=s1*s2*s3;
y=(J2^3-c*J3^2)^(1/6); %Drucker
% y=(J2^(3/2)-c*J3)^(1/3); %Cazacu-Barlat2004
% y=(J2*(J2^3-c*J3^2))^(1/8); %Cazacu2018
end
Appendix B: Program for calculation of the convex domain of the newly proposed yield function for differential-anisotropic hardening under biaxial loading
function [] = myfunction( )
clc
clear
tol = 1.0e-6;
angled=0:1:360;
angle=angled*pi/180;
n=length(angled)
delta_ifa=0.001/180*pi;
d_ifa=0.1;
for i=1:n
radius_u=5;
radius_l= 0;
radius_m=(radius_u+radius_l)/2;
dr=radius_u-radius_l;
while dr>tol
c=radius_m*cos(angle(i));
d=radius_m*sin(angle(i));
param=[c d];
convex=1; %1 is for convex
for l=0:d_ifa:360
angld_A=l/180*pi;
s11_A_try=cos(angld_A);
s22_A_try=sin(angld_A);
f_A=Yldfunction(s11_A_try,s22_A_try,param);
s11_A=s11_A_try/f_A;
s22_A=s22_A_try/f_A;
angld_B=l/180*pi+delta_ifa;
s11_B_try=cos(angld_B);
s22_B_try=sin(angld_B);
f_B=Yldfunction(s11_B_try,s22_B_try,param);
s11_B=s11_B_try/f_B;
s22_B=s22_B_try/f_B;
s11_C=(s11_A+s11_B)/2;
s22_C=(s22_A+s22_B)/2;
f_C=Yldfunction(s11_C,s22_C,param);
if f_C>1
convex=0;
break
end
end
if convex==1
radius_l=radius_m;
else
radius_u=radius_m;
end
radius_m=(radius_u+radius_l)/2;
dr=(radius_u-radius_l)/2;
end
c_limit(i)=radius_m*cos(angle(i));
d_limit(i)=radius_m*sin(angle(i));
end
plot(c_limit,d_limit,'b-','linewidth',5)
end
function y=Yldfunction(s11,s22,param)
c=param(1);
d=param(2);
I1=s11+s22;
s1=s11-I1/3;
s2=s22-I1/3;
s3=-I1/3;
J2=-s1*s2-s2*s3-s1*s3;
J3=s1*s2*s3;
y=((J2^3-c*J3^2)^(1/2)-d*J3)^(1/3);
end
Appendix C: Program for calculation of the convex domain of the newly proposed yield function for differential-anisotropic hardening on the p-plane
function [] = myfunction( )
clc
clear
tol=1.0e-6;
angled=0:1:360;
angle=angled*pi/180;
n=length(angled)
delta_ifa=0.001/180*pi;
d_ifa=0.1;
for i=1:n
radius_u =5;
radius_l = 0;
radius_m=(radius_u+radius_l)/2;
dr=radius_u-radius_l;
while dr>tol
c=radius_m*cos(angle(i));
d=radius_m*sin(angle(i));
param=[c d];
convex=1; %1 is for convex
for l=0:d_ifa:360
angld_A=l/180*pi;
s1_A_try=2/3*cos(angld_A);
s2_A_try=2/3*cos(2/3*pi-angld_A);
s3_A_try=2/3*cos(4/3*pi-angld_A);
f_A=Yldfunction(s1_A_try,s2_A_try,s3_A_try,param);
s1_A=s1_A_try/f_A;
s2_A=s2_A_try/f_A;
s3_A=s3_A_try/f_A;
angld_B=l/180*pi+delta_ifa;
s1_B_try=2/3*cos(angld_B);
s2_B_try=2/3*cos(2/3*pi-angld_B);
s3_B_try=2/3*cos(4/3*pi-angld_B);
f_B=Yldfunction(s1_B_try,s2_B_try,s3_B_try,param);
s1_B=s1_B_try/f_B;
s2_B=s2_B_try/f_B;
s3_B=s3_B_try/f_B;
s1_C=(s1_A+s1_B)/2;
s2_C=(s2_A+s2_B)/2;
s3_C=(s3_A+s3_B)/2;
f_C=Yldfunction(s1_C,s2_C,s3_C,param);
if f_C>1
convex=0;
break
end
end
if convex==1
radius_l=radius_m;
else
radius_u=radius_m;
end
radius_m=(radius_u+radius_l)/2;
dr=(radius_u-radius_l)/2;
end
c_limit(i)=radius_m*cos(angle(i));
d_limit(i)=radius_m*sin(angle(i));
end
plot(c_limit,d_limit,'b-','linewidth',5)
end
function y=Yldfunction(s1,s2,s3,param)
c=param(1);
d=param(2);
J2=-s1*s2-s2*s3-s1*s3;
J3=s1*s2*s3;
y=((J2^3-c*J3^2)^(1/2)-d*J3)^(1/3);
end
Appendix D: Program for calculating the convex domain of the Hu2017 normalized stress invariant-based yield function
function [] = myfunction( )
clc
clear
tol=1.0e-6;
inter=1;
ifad=0:inter:360;
ifa=ifad*pi/180;
bitad=0:inter:180;
bita=bitad*pi/180;
n_ifa=length(ifa)
n_bita=length(bita)
delta_ifa=0.001/180*pi;
d_ifa=1;
for i=1:n_ifa
for j=1:n_bita
radius_u=50;
radius_l=0;
radius_m=(radius_u+radius_l)/2;
dr=radius_u-radius_l;
while dr>tol
a=radius_m*sin(bita(j))*cos(ifa(i));
b=radius_m*sin(bita(j))*sin(ifa(i));
c=radius_m*cos(bita(j));
param=[a b c];
convex=1; %1 is for convex
for l=0:d_ifa:360
angld_A=l/180*pi;
s1_A_try=2/3*cos(angld_A);
s2_A_try=2/3*cos(2/3*pi-angld_A);
s3_A_try=2/3*cos(4/3*pi-angld_A);
f_A=Yldfunction(s1_A_try,s2_A_try,s3_A_try,param);
s1_A=s1_A_try/f_A;
s2_A=s2_A_try/f_A;
s3_A=s3_A_try/f_A;
angld_B=l/180*pi+delta_ifa;
s1_B_try=2/3*cos(angld_B);
s2_B_try=2/3*cos(2/3*pi-angld_B);
s3_B_try=2/3*cos(4/3*pi-angld_B);
f_B=Yldfunction(s1_B_try,s2_B_try,s3_B_try,param);
s1_B=s1_B_try/f_B;
s2_B=s2_B_try/f_B;
s3_B=s3_B_try/f_B;
s1_C=(s1_A+s1_B)/2;
s2_C=(s2_A+s2_B)/2;
s3_C=(s3_A+s3_B)/2;
f_C=Yldfunction(s1_C,s2_C,s3_C,param);
if f_C>1
convex=0;
break
end
end
if convex==1
radius_l=radius_m;
else
radius_u=radius_m;
end
radius_m=(radius_u+radius_l)/2;
dr=(radius_u-radius_l)/2;
end
a_limit(i,j)=radius_m*sin(bita(j))*cos(ifa(i));
b_limit(i,j)=radius_m*sin(bita(j))*sin(ifa(i));
c_limit(i,j)=radius_m*cos(bita(j));
end
end
mesh(a_limit,b_limit,c_limit);
end
function y=Yldfunction(s1,s2,s3,param)
a=param(1);
b=param(2);
c=param(3);
J2=-s1*s2-s2*s3-s1*s3;
J3=s1*s2*s3;
y=J2^(1/2)-a*J3/J2-b*J3^2/J2^(5/2)-c*J3^3/J2^4;
end
-
2024
05-16
-
2024
05-15
-
2024
05-15
-
2024
05-15
-
2024
05-15
-
2024
05-15