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

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

学院：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、、G12；各单层对应的厚度；各单层对应的主方向夹角。然后就要计算每个单层板的二维刚度矩阵Q，具体公式如下：
；；；；
得到Q矩阵后，根据课本上讲到的得到。
然后根据z坐标的定义求出到，接下来，最重要的一步，根据下式计算A、B、D。

一、书上P110的几个问题可以归纳为以下几个类型。
（1）正交铺设5层对称层合板（T5-7）
数据文档
层数 5
层序数 厚度m E1（Pa） E2（Pa） v12 v21 G12（Pa） 角度（°）
1 1.00E-03 9.60E+10 2.40E+10 0.10 0.40 1.00E+10 0.00
2 1.00E-03 9.60E+10 2.40E+10 0.10 0.40 1.00E+10 90.00
3 1.00E-03 9.60E+10 2.40E+10 0.10 0.40 1.00E+10 0.00
4 1.00E-03 9.60E+10 2.40E+10 0.10 0.40 1.00E+10 90.00
5 1.00E-03 9.60E+10 2.40E+10 0.10 0.40 1.00E+10 0.00

结果文档
拉伸刚度A=
3.5000E+08 5.0000E+07 -4.3711E-01
5.0000E+07 2.7500E+08 -6.1196E+00
-4.3711E-01 -6.1196E+00 5.0000E+07
耦合刚度B=
-9.0523E-12 -7.2831E-13 -3.8514E-21
-7.2831E-13 -2.2631E-12 -7.0304E-20
-3.8514E-21 -9.1480E-20 -9.0417E-13
弯曲刚度D=
8.7917E+02 1.0417E+02 -4.7354E-07
1.0417E+02 4.2292E+02 -6.6296E-06
-4.7354E-07 -6.6296E-06 1.0417E+02
由此可以从课本上了解到的：
A16=A26=0；
D16=D26=0；
相吻合。
这里B显然是等于零的。

（2） 正交铺设6层反对称层合板（T5-8）
数据文档
层数 6
层序数 厚度m E1（Pa） E2（Pa） v12 v21 G12（Pa） 角度（°）
1 1.00E-03 9.60E+10 2.40E+10 0.10 0.40 1.00E+10 0.00
2 1.00E-03 9.60E+10 2.40E+10 0.10 0.40 1.00E+10 90.00
3 1.00E-03 9.60E+10 2.40E+10 0.10 0.40 1.00E+10 0.00
4 1.00E-03 9.60E+10 2.40E+10 0.10 0.40 1.00E+10 90.00
5 1.00E-03 9.60E+10 2.40E+10 0.10 0.40 1.00E+10 0.00
6 1.00E-03 9.60E+10 2.40E+10 0.10 0.40 1.00E+10 90.00

结果文档
拉伸刚度A=
3.7500E+08 6.0000E+07 -6.5567E-01
6.0000E+07 3.7500E+08 -9.1794E+00
-6.5567E-01 -9.1794E+00 6.0000E+07
耦合刚度B=
-1.1250E+05 1.8943E-10 -3.2784E-04
1.8943E-10 1.1250E+05 -4.5897E-03
-3.2784E-04 -4.5897E-03 1.8102E-10
弯曲刚度D=
1.1250E+03 1.8000E+02 -1.9670E-06
1.8000E+02 1.1250E+03 -2.7538E-05
-1.9670E-06 -2.7538E-05 1.8000E+02

由此可以和课本上了解到的：
A11=A22；
D11=D22；
A16=A26=D16=D26=0；
相吻合。

（3）5层对称角铺设层合板（T5-9）
数据文档
层数 5
层序数 厚度m E1（Pa） E2（Pa） v12 v21 G12（Pa） 角度（°）
1 1.00E-03 9.60E+10 2.40E+10 0.10 0.40 1.00E+10 30.00
2 1.00E-03 9.60E+10 2.40E+10 0.10 0.40 1.00E+10 -30.00
3 1.00E-03 9.60E+10 2.40E+10 0.10 0.40 1.00E+10 30.00
4 1.00E-03 9.60E+10 2.40E+10 0.10 0.40 1.00E+10 -30.00
5 1.00E-03 9.60E+10 2.40E+10 0.10 0.40 1.00E+10 30.00

结果文档
拉伸刚度A=
3.4531E+08 1.1094E+08 2.3274E+07
1.1094E+08 1.5781E+08 9.2015E+06
2.3274E+07 9.2015E+06 1.1094E+08
耦合刚度B=
-1.2975E-11 -1.4779E-12 1.4282E-12
-8.4022E-12 1.3998E-12 -6.5192E-13
1.4282E-12 -8.2956E-13 -8.4022E-12
弯曲刚度D=
7.1940E+02 2.3112E+02 1.4159E+02
2.3112E+02 3.2878E+02 5.5976E+01
1.4159E+02 5.5976E+01 2.3112E+02

由此可以和课本上了解到的：
A16、A26相对要小；
D16、D26相对要小；
相吻合。
这里B显然是等于零的。

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

数据文档
层数 6
层序数 厚度m E1（Pa） E2（Pa） v12 v21 G12（Pa） 角度（°）
1 1.00E-03 9.60E+10 2.40E+10 0.10 0.40 1.00E+10 30.00
2 1.00E-03 9.60E+10 2.40E+10 0.10 0.40 1.00E+10 -30.00
3 1.00E-03 9.60E+10 2.40E+10 0.10 0.40 1.00E+10 30.00
4 1.00E-03 9.60E+10 2.40E+10 0.10 0.40 1.00E+10 -30.00
5 1.00E-03 9.60E+10 2.40E+10 0.10 0.40 1.00E+10 30.00
6 1.00E-03 9.60E+10 2.40E+10 0.10 0.40 1.00E+10 -30.00

结果文档
拉伸刚度A=
4.1437E+08 1.3313E+08 -1.3479E-09
1.3313E+08 1.8938E+08 3.0013E-10
-1.3479E-09 3.4470E-10 1.3313E+08
耦合刚度B=
1.2989E-11 -9.3543E-12 -6.9823E+04
7.4607E-12 -2.7711E-12 -2.7605E+04
-6.9823E+04 -2.7605E+04 7.4607E-12
弯曲刚度D=
1.2431E+03 3.9938E+02 2.9643E-14
3.9938E+02 5.6813E+02 9.5583E-15
2.9643E-14 7.4281E-15 3.9938E+02

由此可以和课本上了解到的：
A16=A26=D16=D26=0;
B11=B22=B12=B21=B66=0；
相吻合。

（5）我还想验证一个书上的例题，在课本P114。三层层合板，外层厚度t1，内层10t1，正交铺设比m=0.2,。玻璃/环氧单层板性能：E1=5.4E10Pa，E2=1.8E10Pa，v21=0.25，G12=8.8E9Pa。
数据文档
层数 3
层序数 厚度m E1（Pa） E2（Pa） v12 v21 G12（Pa） 角度（°）
1 1.00E+00 5.40E+10 1.80E+10 0.083 0.250 8.80E+09 0.00
2 1.00E+01 5.40E+10 1.80E+10 0.083 0.250 8.80E+09 90.00
3 1.00E+00 5.40E+10 1.80E+10 0.083 0.250 8.80E+09 0.00

结果文档
The NORMALISED extension stiffness A* equals:
2.4509E+10 4.5954E+09 -1.3893E+02
4.5954E+09 4.9017E+10 -1.2002E+03
-1.3893E+02 -1.2002E+03 8.8000E+09
The NORMALISED coupling tensor C equals:
-9.3609E+09 1.4411E-05 -4.2450E+01
1.4411E-05 9.3609E+09 -3.6673E+02
-4.2450E+01 -3.6673E+02 1.6107E-05
The NORMALISED bending stiffness D* equals:
3.3869E+10 4.5954E+09 -9.6477E+01
4.5954E+09 3.9656E+10 -8.3347E+02
-9.6477E+01 -8.3347E+02 8.8000E+09

由于课本上只是分析Nx的荷载，给出了A*的数值和课本上计算的结果一致。

二、验证Verchery的论文里给出的数值算例。
这里一直到Table5的数据都是从Verchery的论文中截获。

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

（1）Boron-Epoxy材料。

（Boron-Epoxy）数据文档
层数 18
层序数 厚度m E1（Pa） E2（Pa） v12 v21 G12（Pa） 角度（°）
1 1.00E-03 2.04E+11 1.85E+10 0.021 0.230 5.59E+09 0.00
2 1.00E-03 2.04E+11 1.85E+10 0.021 0.230 5.59E+09 0.00
3 1.00E-03 2.04E+11 1.85E+10 0.021 0.230 5.59E+09 60.00
4 1.00E-03 2.04E+11 1.85E+10 0.021 0.230 5.59E+09 60.00
5 1.00E-03 2.04E+11 1.85E+10 0.021 0.230 5.59E+09 -60.00
6 1.00E-03 2.04E+11 1.85E+10 0.021 0.230 5.59E+09 -60.00
7 1.00E-03 2.04E+11 1.85E+10 0.021 0.230 5.59E+09 -60.00
8 1.00E-03 2.04E+11 1.85E+10 0.021 0.230 5.59E+09 60.00
9 1.00E-03 2.04E+11 1.85E+10 0.021 0.230 5.59E+09 0.00
10 1.00E-03 2.04E+11 1.85E+10 0.021 0.230 5.59E+09 -60.00
11 1.00E-03 2.04E+11 1.85E+10 0.021 0.230 5.59E+09 -60.00
12 1.00E-03 2.04E+11 1.85E+10 0.021 0.230 5.59E+09 60.00
13 1.00E-03 2.04E+11 1.85E+10 0.021 0.230 5.59E+09 60.00
14 1.00E-03 2.04E+11 1.85E+10 0.021 0.230 5.59E+09 0.00
15 1.00E-03 2.04E+11 1.85E+10 0.021 0.230 5.59E+09 0.00
16 1.00E-03 2.04E+11 1.85E+10 0.021 0.230 5.59E+09 0.00
17 1.00E-03 2.04E+11 1.85E+10 0.021 0.230 5.59E+09 60.00
18 1.00E-03 2.04E+11 1.85E+10 0.021 0.230 5.59E+09 -60.00

（Boron-Epoxy）结果文档
The stiffness of number 1 ply is:
2.0499E+11 4.2757E+09 0.0000E+00
4.2757E+09 1.8590E+10 0.0000E+00
0.0000E+00 0.0000E+00 5.5900E+09
The coupling stiffness B equals:
-5.6581E-10 2.1876E-10 -1.4447E-10
2.1876E-10 6.5563E-10 -5.0525E-10
-1.4447E-10 -4.5370E-10 2.5395E-10
The NORMALISED extension stiffness A* equals:
8.7706E+10 2.8359E+10 -6.5585E-07
2.8359E+10 8.7706E+10 -2.9973E-06
-6.5585E-07 -3.0072E-06 2.9674E+10
The NORMALISED coupling tensor C equals:
-2.1719E+10 4.4599E+09 -1.4444E+09
4.4599E+09 1.2799E+10 -4.5344E+09
-1.4444E+09 -4.5344E+09 4.4599E+09
The NORMALISED bending stiffness D* equals:
1.0943E+11 2.3899E+10 1.4444E+09
2.3899E+10 7.4907E+10 4.5344E+09
1.4444E+09 4.5344E+09 2.5214E+10

这里的结果显然是B=0，而且得到的Q，A*，D*与论文上的数据一致，

（2） Glass-Epoxy材料。

（Glass-Epoxy）数据文档
层数 18
层序数 厚度m E1（Pa） E2（Pa） v12 v21 G12（Pa） 角度（°）
1 1.00E-03 3.86E+10 8.27E+09 0.056 0.260 4.14E+09 0.00
2 1.00E-03 3.86E+10 8.27E+09 0.056 0.260 4.14E+09 0.00
3 1.00E-03 3.86E+10 8.27E+09 0.056 0.260 4.14E+09 60.00
4 1.00E-03 3.86E+10 8.27E+09 0.056 0.260 4.14E+09 60.00
5 1.00E-03 3.86E+10 8.27E+09 0.056 0.260 4.14E+09 -60.00
6 1.00E-03 3.86E+10 8.27E+09 0.056 0.260 4.14E+09 -60.00
7 1.00E-03 3.86E+10 8.27E+09 0.056 0.260 4.14E+09 -60.00
8 1.00E-03 3.86E+10 8.27E+09 0.056 0.260 4.14E+09 60.00
9 1.00E-03 3.86E+10 8.27E+09 0.056 0.260 4.14E+09 0.00
10 1.00E-03 3.86E+10 8.27E+09 0.056 0.260 4.14E+09 -60.00
11 1.00E-03 3.86E+10 8.27E+09 0.056 0.260 4.14E+09 -60.00
12 1.00E-03 3.86E+10 8.27E+09 0.056 0.260 4.14E+09 60.00
13 1.00E-03 3.86E+10 8.27E+09 0.056 0.260 4.14E+09 60.00
14 1.00E-03 3.86E+10 8.27E+09 0.056 0.260 4.14E+09 0.00
15 1.00E-03 3.86E+10 8.27E+09 0.056 0.260 4.14E+09 0.00
16 1.00E-03 3.86E+10 8.27E+09 0.056 0.260 4.14E+09 0.00
17 1.00E-03 3.86E+10 8.27E+09 0.056 0.260 4.14E+09 60.00
18 1.00E-03 3.86E+10 8.27E+09 0.056 0.260 4.14E+09 -60.00

（Glass-Epoxy）结果文档
The stiffness of number 1 ply is:
3.9170E+10 2.1820E+09 0.0000E+00
2.1820E+09 8.3922E+09 0.0000E+00
0.0000E+00 0.0000E+00 4.1400E+09
The coupling stiffness B equals:
-3.2222E-10 8.5606E-11 -3.2021E-11
8.5606E-11 2.1697E-10 -8.1172E-11
-3.4829E-11 -8.1172E-11 -2.0606E-13 The The NORMALISED extension stiffness A* equals:
2.0451E+10 5.5118E+09 -1.4616E-07
5.5118E+09 2.0451E+10 -4.3469E-07
-1.4491E-07 -4.3469E-07 7.4698E+09
The NORMALISED coupling tensor C equals:
-3.4665E+09 6.1663E+08 -2.8000E+08
6.1663E+08 2.2332E+09 -7.0721E+08
-2.8000E+08 -7.0721E+08 6.1663E+08
The NORMALISED bending stiffness D* equals:
2.3918E+10 4.8952E+09 2.8000E+08
4.8952E+09 1.8218E+10 7.0721E+08
2.8000E+08 7.0721E+08 6.8532E+09

这里的结果显然是B=0，而且得到的Q，A*，D*仍然与论文上的数据一致。

（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。验证成功。

第 14 页 共 14 页

（A）第（22）排序。
（Boron-Epoxy）结果文档
The ACTUAL stiffness tensor of the laminate:
The extension stiffness A equals:
1.5787E+09 5.1047E+08 -1.9256E-08
5.1047E+08 1.5787E+09 -3.9050E-08
-1.9256E-08 -5.4129E-08 5.3412E+08

The coupling stiffness B equals:
-4.0463E+06 8.3088E+05 4.8750E+05
8.3088E+05 2.3845E+06 1.5303E+06
4.8750E+05 1.5303E+06 8.3088E+05

The bending stiffness D equals:
4.2625E+04 1.3783E+04 -7.3708E-13
1.3783E+04 4.2625E+04 -4.1602E-12
-7.3708E-13 -3.8027E-12 1.4421E+04

The NORMALISED stiffness tensor of the laminate:
The NORMALISED extension stiffness A* equals:
8.7706E+10 2.8359E+10 -1.0698E-06
2.8359E+10 8.7706E+10 -2.1694E-06
-1.0698E-06 -3.0072E-06 2.9674E+10

The NORMALISED coupling tensor C equals:
1.4275E-05 5.7798E-06 4.4685E-07
5.7798E-06 8.9034E-06 6.3907E-06
4.4685E-07 4.8174E-06 2.6114E-06

The NORMALISED bending stiffness D* equals:
8.7706E+10 2.8359E+10 -1.5166E-06
2.8359E+10 8.7706E+10 -8.5601E-06
-1.5166E-06 -7.8246E-06 2.9674E+10

（B）第（23）排序。
（Boron-Epoxy）结果文档
The ACTUAL stiffness tensor of the laminate:
The extension stiffness A equals:
1.5787E+09 5.1047E+08 -2.2981E-08
5.1047E+08 1.5787E+09 -5.3951E-08
-2.2981E-08 -5.4129E-08 5.3412E+08

The coupling stiffness B equals:
-1.0897E-09 2.1876E-10 -1.7357E-10
2.1876E-10 4.2280E-10 -4.4705E-10
-1.7357E-10 -4.5370E-10 3.7036E-10

The bending stiffness D equals:
4.2625E+04 1.3783E+04 -7.3708E-13
1.3783E+04 4.2625E+04 -2.3412E-12
-7.3708E-13 -2.8932E-12 1.4421E+04

The NORMALISED stiffness tensor of the laminate:
The NORMALISED extension stiffness A* equals:
8.7706E+10 2.8359E+10 -1.2767E-06
2.8359E+10 8.7706E+10 -2.9973E-06
-1.2767E-06 -3.0072E-06 2.9674E+10

The NORMALISED coupling tensor C equals:
-6.9290E-07 5.3458E-06 2.3989E-07
5.3458E-06 1.0632E-05 1.8201E-06
2.3989E-07 2.9460E-06 1.0096E-05

The NORMALISED bending stiffness D* equals:
8.7706E+10 2.8359E+10 -1.5166E-06
2.8359E+10 8.7706E+10 -4.8174E-06
-1.5166E-06 -5.9532E-06 2.9674E+10

附件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
WRITE(9,200)Q(I,:,K)
200 FORMAT(ES12.4,6X,ES12.4,6X,ES12.4)
ENDDO
WRITE(9,"(/)")
A=A+Q(:,:,K)*Z1(K)
B=B+Q(:,:,K)*Z2(K)
D=D+Q(:,:,K)*Z3(K)
ENDDO
!Output the actual stiffness of the laminate
WRITE(9,"(/)");WRITE(9,"(/)")
WRITE(9,*)"The ACTUAL stiffness tensor of the laminate:"
WRITE(9,"(/)")
WRITE(9,*)"The extension stiffness A equals:"
DO I=1,3
WRITE(9,200)A(I,1:3)
ENDDO
WRITE(9,"(/)")
WRITE(9,*)"The c