复合材料的MATLAB计算M文件

%MOCMexp01.m %MECHANICS OF COMPOSITE MATERIALS EX1 %已知：层合板[θ/-θ]s,其单层板性能指标为E1=140GPa,E2=10GPa,G12=5GPa,u12=0.3,tp=0.125mm, %求弹性系数（弯曲和拉伸两组）Ex,Gxy,uxy,mx与theta的关系（θ为0~90°，Δθ=5°） %对结果进行列表，与单层板结果进行对比论证 %清屏，清除原有变量，程序开始 echo off; clc; clear; format short; %………………………………………………………………………………………………………… %1、预处理，定义已知材料性能指标及相关参数 E1=140; E2=10; G12=5; u12=0.3; tk=[0.125,0.125,0.125,0.125]; %对应1-4层单层板的厚度，单位mm zk=[-0.1875,-0.0625,0.0625,0.187]; %对应1-4层单层板的中面坐标，单位mm t=0.5; %对应层合板的总厚度，单位mm %由单层板性能指标求刚度系数Q Q=zeros(3,3); u21=u12*E2/E1; Q(1,1)=E1/(1-u12*u21); Q(2,2)=E2/(1-u12*u21); Q(3,3)=G12; Q(1,2)=u12*Q(2,2); Q(2,1)=Q(1,2); Q(1,3)=0;Q(3,1)=0;Q(2,3)=0;Q(3,2)=0; %此处单位为GPa %………………………………………………………………………………………………………… %2、计算，引入变量θ（0~90°）,以Δθ=5°的步进循环求值 for theta=0:5:90 c=theta/5+1; %c为计数变量 %转换成变换刚度矩阵Q（θ） x=[theta -theta]; m=cosd(x);n=sind(x);%以角度制计算，对7.0版本使用，若为6.0请改为弧度制 for k=1:2 T{k}=[m(k)^2,n(k)^2,2*m(k)*n(k);n(k)^2,m(k)^2,-2*m(k)*n(k);-m(k)*n(k),m(k)*n(k),(m(k)^2-n(k)^2)]; Qk{k}=(T{k}'*Q^-1*T{k})^-1; %T'为矩阵T的转置，Q^-1为矩阵Q的求逆 end; %由对称性 Qk{3}=Qk{2};Qk{4}=Qk{1}; %定义细胞元组矩阵T,Qk，方便计算矩阵A,D %此处单位为GPa %计算拉伸刚度A A=zeros(3,3); %赋初值0 for i=1:3, for j=1:3, for k=1:4, A(i,j)=A(i,j)+tk(k)*Qk{k}(i,j); end; end; end; %此处单位为KN/mm %计算弯曲刚度D D=zeros(3,3); for i=1:3, for j=1:3, for k=1:4, D(i,j)=D(i,j)+(tk(k)*zk(k)^2+tk(k)^3/12)*Qk{k}(i,j); end; end; end; %此处单位为(KN/mm)-1 %求单层板的弹性系数 Ex0(c)=(m(1)^4/E1+n(1)^4/E2+m(1)^2*n(1)^2*(1/G12-2*u12/E1))^-1; Gxy0(c)=(m(1)^2*n(1)^2*(4/E1+4/E2+8*u12/E1)+(m(1)^2-n(1)^2)^2*1/G12)^-1; uxy0(c)=Ex0(c)*((m(1)^4+n(1)^4)*u12/E1-m(1)^2*n(1)^2*(1/E1+1/E2-1/G12)); mx0(c)=Ex0(c)*(m(1)^3*n(1)*(1/G12-2*u12/E1-2/E1)-m(1)*n(1)^3*(1/G12-2*u12/E1-2/E2)); %求层合板的面内及弯曲弹性系数 a=A^-1;d=D^-1; %对矩阵A，D求逆 Ex1(c)=1/(t*a(1,1));Ex2(c)=12/(t^3*d(1,1)); %Ey1(c)=1/(t*a(2,2));Ey2(c)=12/(t^3*d(2,2)); Gxy1(c)=1/(t*a(3,3));Gxy2(c)=12/(t^3*d(3,3)); uxy1(c)=-a(1,2)/a(1,1);uxy2(c)=-d(1,2)/d(1,1); %uyx1(c)=-a(1,2)/a(2,2);uyx2(c)=-d(1,2)/d(2,2); mx1(c)=-a(1,3)/a(1,1);mx2(c)=-d(1,3)/d(1,1); %my1(c)=-a(2,3)/a(2,2);my2(c)=-d(2,3)/d(1,1); end; %………………………………………………………………………………………………………… %3、后处理，绘图，结果可视化分析 %先画出单层板的弹性系数图 figure('NumberTitle','off','Name','单层板弹性系数'); theta=[0:5:90]; subplot(2,2,1); plot(theta,Ex0,'g+:'); xlabel('层角θ(0-90°)'); ylabel('Ex/GPa'); title('弹性模量Ex随θ的变化'); axis auto; subplot(2,2,2); plot(theta,Gxy0,'g+:'); xlabel('层角θ(0-90°)'); ylabel('Gxy/GPa'); title('剪切模量Gxy随θ的变化'); axis auto; subplot(2,2,3); plot(theta,uxy0,'g+:'); xlabel('层角θ(0-90°)'); ylabel('μ'); title('柏松比随θ的变化'); axis auto; subplot(2,2,4); plot(theta,mx0,'g+:'); xlabel('层角θ(0-90°)'); ylabel('m'); title('剪切耦合系数随θ的变化'); axis auto; %画出层合板的弹性系数图 %此处，红色为面内弹性系数，蓝色为弯曲弹性系数 figure('NumberTitle','off','Name','层合板弹性系数'); theta=[0:5:90]; subplot(2,2,1); plot(theta,Ex1,'r+:'); hold on; plot(theta,Ex2,'b+:'); xlabel('层角θ(0-90°)'); ylabel('Ex/GPa'); title('弹性模量Ex随θ的变化'); axis auto; hold off; legend('面内','弯曲',1); %曲线标注 subplot(2,2,2); plot(theta,Gxy1,'r+:'); hold on; plot(theta,Gxy2,'b+:'); xlabel('层角θ(0-90°)'); ylabel('Gxy/GPa'); title('剪切模量Gxy随θ的变化'); axis auto; hold off; subplot(2,2,3); plot(theta,uxy1,'r+:'); hold on; plot(theta,uxy2,'b+:'); xlabel('层角θ(0-90°)'); ylabel('μ'); title('柏松比随θ的变化'); axis auto; hold off; subplot(2,2,4); plot(theta,mx1,'r+:'); hold on; plot(theta,mx2,'b+:'); xlabel('层角θ(0-90°)'); ylabel('m'); title('剪切耦合系数随θ的变化'); axis auto; hold off; %将层合板结果输出列表显示 temp=[theta' Ex1' Ex2' Gxy1' Gxy2' uxy1' uxy2' mx1' mx2']; disp(' θ Ex1 Ex2 Gxy1 Gxy2 uxy1 uxy2 mx1 mx2'); disp(temp); %对比分析 disp('由上表及对比两图可知：'); disp('1、该层合板在面内具有正交各向异性，即没有拉剪耦合效应，但存在弯扭耦合。'); disp('2、与单层板相比，层合板随θ的增加剪切模量先逐渐增大然后减小，相对于单层板有更高的抗剪切能力。'); %清除变量，回收内存，结束程序 clear; pack;