 宏观力学层合板强度分析MATLAB实现

最新推荐文章于 2024-01-03 10:51:36 发布

仰望星空1.0

最新推荐文章于 2024-01-03 10:51:36 发布

阅读量286

点赞数

文章标签： matlab 开发语言

本文链接：https://blog.csdn.net/2301_76701508/article/details/134292112

版权

n=input('please input the layer number of compusite materials:');

t=zeros(1,n);

theta=zeros(1,n);

e1=input('please input e1(Gpa): ')*10^9;

e2=input('please input e2(Gpa): ')*10^9;

g12=input('please input g12(Gpa): ')*10^9;

v21=input('please input v21: ');

F=zeros(5,1);

F(1,1)=input('please input Xt(unit/Mpa):')*10^6;

F(2,1)=input('please input Xc(unit/Mpa):')*10^6;

F(3,1)=input('please input Yt(unit/Mpa):')*10^6;

F(4,1)=input('please input Yc(unit/Mpa):')*10^6;

F(5,1)=input('please input S(unit/Mpa):')*10^6;

for i=1:n

disp(['please input elasetic constant from bottom to the top(The direction of the elastic spindle is' ...

' specified as 1 direction),while now you should input layer ' ,num2str(i)])

theta(1,i)=input('please input theta(°）: ');

t(1,i)=input('please input the thickness(m):');

theta(1,i)=theta(1,i)/180*pi;

end

N=zeros(1,3);M=zeros(1,3);

for i=1:3

N(i)=input(['please input N',num2str(i),'(force at per width). While axis 1 is the fiber direction(unit/N):']);

end

for i=1:3

M(i)=input(['please input M',num2str(i),'(force at per width). While axis 1 is the fiber direction(unit/N*m):']);

end

FAM=[N';M'];

scale0=zeros(1,n);

[q,s]=calculateq(e1,e2,g12,v21);

[q0,s0,T]=calculateq0(n,q,s,theta);

ddd=0;

epsilonx=zeros(1,3*n);

format long

for i=1:3*n

[Q0,S0,A,B,D,z0]=calculateQtotal(q0,t,n);

[scale,deslayer,desnum,layer,sc,epsilon,kappa]=calculatescale(q0,z0,T,F,n,t,S0,FAM);

[q0,desd]=recalculateq0(t,q,q0,T,F,scale,epsilon,kappa,z0,desnum,deslayer,layer);

scale0(i)=scale;

DF=scale0(i)*FAM;

if desd(1)==0

break

else if i>1&&scale0(i)*2<scale0(i-1)

break

end

epsilonx(2*i-1)=scale0(i)*epsilon(1);

[Q0,S0,A,B,D,z0]=calculateQtotal(q0,t,n);

ttt=scale0(i)*S0*FAM;

epsilonx(2*i)=ttt(1);

ddd=ddd+1;

disp(['第',num2str(i),'次总破坏层数为',num2str(desnum)])

for j=1:desnum

disp(['第',num2str(deslayer(j)),'层'])

disp(['该层破坏的方向为',num2str(desd(j)),'方向'])

end

disp('破坏时的外力矩阵为')

disp(num2str(DF,'%e'))

end

scale1=scale0;

clear scale0;

for i=1:ddd

scale0(2*i-1)=scale1(i);

scale0(2*i)=scale1(i);

epsilonx0(2*i-1)=epsilonx(2*i-1);

epsilonx0(2*i)=epsilonx(2*i);

end

for i=2*ddd:-1:2

scale0(i)=scale0(i-1);

epsilonx0(i)=epsilonx0(i-1);

end

scale0(1)=0;

epsilonx0(1)=0;

plot(epsilonx0,scale0)

xlabel('\epsilon');

ylabel('Fx/N','Rotation',90);

function [q,s]=calculateq(e1,e2,g12,v21)

s=zeros(3,3);

q=zeros(3,3);

s(1,1)=1/e1;

s(2,2)=1/e2;

s(1,2)=-v21/e1;

s(2,1)=s(1,2);

s(3,3)=1/g12;

q(1,1)=s(2,2)/(s(1,1)*s(2,2)-s(1,2)^2);

q(2,2)=s(1,1)/(s(1,1)*s(2,2)-s(1,2)^2);

q(1,2)=-s(1,2)/(s(1,1)*s(2,2)-s(1,2)^2);

q(2,1)=q(1,2);

q(3,3)=1/s(3,3);

end

function [q0,s0,T]=calculateq0(n,q,s,theta)

% 计算各层刚度矩阵q0(3,3,n)和柔度矩阵s0(3,3,n)(转角后）

q0=zeros(3,3,n);

s0=zeros(3,3,n);

T=zeros(3,3,n);

for i=1:n

T(1,1,i)=cos(theta(1,i))^2;

T(1,2,i)=sin(theta(1,i))^2;

T(1,3,i)=2*sin(theta(1,i))*cos(theta(1,i));

T(2,1,i)=sin(theta(1,i))^2;

T(2,2,i)=cos(theta(1,i))^2;

T(2,3,i)=-2*sin(theta(1,i))*cos(theta(1,i));

T(3,1,i)=-sin(theta(1,i))*cos(theta(1,i));

T(3,2,i)=sin(theta(1,i))*cos(theta(1,i));

T(3,3,i)=cos(theta(1,i))^2-sin(theta(1,i))^2;

qt=(T(:,:,i)^-1)*q*((T(:,:,i)^-1)');

st=T(:,:,i)'*s*T(:,:,i);

for m=1:3

for n=1:3

q0(m,n,i)=qt(m,n);

s0(m,n,i)=st(m,n);

end

function [Q0,S0,A,B,D,z0]=calculateQtotal(q0,t,n)

A=zeros(3,3);

B=zeros(3,3);

D=zeros(3,3);

t0=0;

for i=1:n

t0=t0+t(i);

end

z0=t0/2;

tz=0;

for k=1:n

tz=tz+t(k);

for i=1:3

for j=1:3

A(i,j)=A(i,j)+q0(i,j,k)*t(k);

B(i,j)=B(i,j)+0.5*q0(i,j,k)*((tz-z0)^2-(tz-t(k)-z0)^2);

D(i,j)=D(i,j)+1/3*q0(i,j,k)*((tz-z0)^3-(tz-t(k)-z0)^3);

end

Q0=[A,B;B,D];

S0=Q0^(-1);

% S0=zeros(6,6);

% S0(1:3,1:3)=A^-1;

z0=t0/2;

end

function [scale,deslayer,desnum,layer,sc,epsilon,kappa]=calculatescale(q0,z0,T,F,n,t,S0,FAM)

%判断最易破坏层，初步分析可确定破坏平面为两层之间的连接面。输出最易破坏层的放大因子scale（与加载力相比，放大因子）及层数desnum、deslayer

syms scalebelow scaleup;

tz=0;

sigma=zeros(3,1,n);

sc=zeros(2*n,1);

EAK0=S0*FAM;

epsilon=zeros(3,1);

kappa=zeros(3,1);

for i=1:3

epsilon(i,1)=EAK0(i);

kappa(i,1)=EAK0(i+3);

end

for i=1:n

% if tz<=0

q0tr=q0(:,:,i);

sigmatrbelow=scalebelow*q0tr*(epsilon+(tz-z0)*kappa);%计算i层下表面sigma矩阵

tz=tz+t(i);

sigmatrup=scaleup*q0tr*(epsilon+(tz-z0)*kappa);%计算i层上表面sigma矩阵

sigmabelow0=T(:,:,i)*sigmatrbelow;

sigmaup0=T(:,:,i)*sigmatrup;

if q0tr==zeros(3,3)

continue

end

%Hill-蔡强度理论判断

fbelow=sigmabelow0(1)^2/(F(1)*F(2))-sigmabelow0(1)*sigmabelow0(2)/(F(1)*F(2))+sigmabelow0(2)^2/(F(3)*F(4))+(F(2)-F(1))/(F(1)*F(2))*sigmabelow0(1)+(F(4)-F(3))/(F(3)*F(4))*sigmabelow0(2)+sigmabelow0(3)^2/F(5)^2-1;

sc(2*i-1)=double(max(solve(fbelow==0,scalebelow)));%计算下表面的应力系数

fup=sigmaup0(1)^2/(F(1)*F(2))-sigmaup0(1)*sigmaup0(2)/(F(1)*F(2))+sigmaup0(2)^2/(F(3)*F(4))+(F(2)-F(1))/(F(1)*F(2))*sigmaup0(1)+(F(4)-F(3))/(F(3)*F(4))*sigmaup0(2)+sigmaup0(3)^2/F(5)^2-1;

sc(2*i)=double(max(solve(fup==0,scaleup)));%计算上表面的应力系数

end

scale=min(sc);

sc1=(sc-min(sc))/(min(sc)+1*10^-10);

desnum=0;

j=1;

for i=1:n

if sc1(2*i-1)<0.01||sc1(2*i)<0.01

desnum=desnum+1;

deslayer(j)=i;

j=j+1;

end