复合材料力学上机编程作业(计算层合板刚度)要点

复合材料力学上机编程作业

学院：School of Civil Engineering专业：Engineering Mechanics

小组成员信息：James Wilson（2012031890015）、Tau Young（2012031890011）复合材料力学学了五个星期，这是这门课的第一次编程作业。我和杨涛结成一个小组，我用的是Fortran编制的程序，Tau Young用的是matlab编制。其中的算例以我的Fortran计算结果为准。Matlab作为可视化界面有其独到之处，在附录2中将会有所展示。

作业的内容是层合板的刚度的计算和验算，包括拉伸刚度A、弯曲刚度D以及耦合刚度B。

首先要给定层合板的各个参数，具体有：层合板的层数N；各单层的弹性常数

E1、E2、υ21、G12；各单层

对应的厚度；各单层对应的主方向夹角θ。然后就要计算每个单层板的二维刚度矩阵Q，具体公式如下：

υ12=υ21E2

E1；Q11=E11-υ12υ21；Q22=E21-υ12υ21；Q12=υ12E1； 1-υ12υ21Q66=G12

得到Q矩阵后，根据课本上讲到的Q=(T-1)TQ(T-1)得到Q。

然后根据z坐标的定义求出z0到zn，接下来，最重要的一步，根据下式计算A、B、D。

n⎧⎪Aij=∑(Qij)k(zk-zk-1)

k=1⎪⎪1n22⎨Bij=∑(Qij)k(zk-zk-1) 2k=1⎪⎪1n33⎪Dij=∑(Qij)k(zk-zk-1)3k=1⎩

一、书上P110的几个问题可以归纳为以下几个类型。

第 1 页共 1 页

（4）6层反对称角铺设层合板（T5-10）第 2 页共 2 页

（5）我还想验证一个书上的例题，在课本P114。三层层合板，外层厚度t1，内层10t1，正交铺设比m=0.2,。

二、验证Verchery的论文里给出的数值算例。

这里一直到Table5的数据都是从Verchery的论文中截获。

Verchery论文中的18层序列，第（21）式【laminates without bending-extension coupling】的排列有两种材料，一种是Boron-Epoxy，一种是Glass-Epoxy。而且都给出了最终的计算结果Q，A*，D*。

下面是我的Fortran计算数据文档和结果文档。

（1）Boron-Epoxy材料。

第 3 页共 3 页

（2）Glass-Epoxy材料。

第 4 页共 4 页

（3）当然我也验证了第（22）【laminates with equal elastic properties in bending and extension】、（23）【quasi-homogeneous laminates】的排序，材料是Boron-Epoxy，下面给出计算的结果。

从下面的两个结果表中可以知道，（22）排序的确是C=0，（23）的排序的确是B=0且C=0。验证成功。

第 5 页共 5 页

附件1：计算所用的程序代码。

PROGRAM COMPOSITE !Coded by James Wilson IMPLICIT NONE

REAL(8)::A(3,3),B(3,3),D(3,3),MC(5),TEMP,ROT(3,3) !A拉伸刚度；B耦合刚度；D弯曲刚度； !MC读入材料常数；ROT旋转矩阵

REAL(8)::TOTAL_TH,HALF_TH !总厚度；半厚度

REAL(8),ALLOCATABLE::Q(:,:,:),AL(:),T(:),Z(:),Z1(:),Z2(:),Z3(:) !Q每层板相应刚度；AL转角；T每层板的厚度；Z坐标量 INTEGER(4)::N,I,J,K,SEQ,L

!____IJK循环变量；N板的层数；SEQ序数

CHARACTER(50)::CHR(8),TEMPC,filename1,filename2 !CHR、TEMPC:character variables

WRITE(*,*)"Please insert the INP file name(a.txt for example):" READ(*,*)filename1 OPEN(8,file=filename1)!Open data file !Read in data

READ(8,*)TEMPC,N

ALLOCATE(Q(3,3,N),AL(N),T(N),Z(N+1),Z1(N),Z2(N),Z3(N)) READ(8,*)CHR(1:8) DO I=1,N

READ(8,*)SEQ,T(I),MC(1:5),AL(I)

Q(:,:,I)=0!Calculate stiffness of each layer for the principal axis TEMP=1./(1-

MC(3)*MC(4)) Q(1,1,I)=MC(1)*TEMP Q(2,2,I)=MC(2)*TEMP Q(3,3,I)=MC(5)

Q(1,2,I)=MC(4)*MC(2)*TEMP Q(2,1,I)=Q(1,2,I)

AL(I)=AL(I)*3.1415926535898/180

ROT(1,1)=(cos(AL(I)))**2!Work out Rot Matrix ROT(2,2)=ROT(1,1)

ROT(3,3)=cos(2*AL(I)) ROT(2,1)=1-ROT(1,1) ROT(1,2)=ROT(2,1)

ROT(3,1)=0.5*sin(2*AL(I)) ROT(3,2)=-ROT(3,1) ROT(1,3)=-2*ROT(3,1)

ROT(2,3)=-2*ROT(3,2)

Q(:,:,I)=MATMUL(MATMUL(ROT,Q(:,:,I)),TRANSPOSE(ROT)) ENDDO

TOTAL_TH=sum(T) HALF_TH=TOTAL_TH/2 Z(1)=-HALF_TH !Work out Z DO

I=1,N

Z(I+1)=Z(I)+T(I) ENDDO

!calculate tensor A、B and D DO K=1,N

Z1(K)=(Z(K+1)-Z(K))

Z2(K)=(Z(K+1)-Z(K))*(Z(K+1)+Z(K))/2 Z3(K)=(Z(K+1)**3-Z(K)**3)/3 ENDDO

A=0;B=0;D=0

WRITE(*,*)"Please insert the OUP file name(b.txt for example):"

READ(*,*)filename2 OPEN(9,file=filename2)

!Write in stiffness tensor for each single ply DO K=1,N

WRITE(9,100)K

100 FORMAT("The stiffness of number",1X,I2,2X,"ply is:") DO I=1,3