复合材料力学上机编程作业(计算层合板刚度).docx
《复合材料力学上机编程作业(计算层合板刚度).docx》由会员分享，可在线阅读，更多相关《复合材料力学上机编程作业(计算层合板刚度).docx（12页珍藏版）》请在冰点文库上搜索。
复合材料力学上机编程作业
学院：
SchoolofCivilEngineering专业：
EngineeringMechanics
小组成员信息：
JamesWilson（201203189001）、TauYoung（201203189001）
复合材料力学学了五个星期，这是这门课的第一次编程作业。
我和杨涛结成一个小组，我用的是Fortran编制的程序，TauYoung用的是matlab编制。
其中的算例以我的Fortran计算结果为准。
Matlab作为可视化界面有其独到之处，在附录2中将会有所展示。
作业的内容是层合板的刚度的计算和验算，包括拉伸刚度A、弯曲刚度D以及耦合刚度B。
首先要给定层合板的各个参数，具体有：
层合板的层数N；各单层的弹性常数E1、E2、、G12；各单层对应的厚度；各单层对应的主方向夹角。
然后就要计算每个单层板的二维刚度矩阵Q，具体公式如下：
；；；；
得到Q矩阵后，根据课本上讲到的得到。
然后根据z坐标的定义求出到，接下来，最重要的一步，根据下式计算A、B、D。
一、书上P110的几个问题可以归纳为以下几个类型。
（1）正交铺设5层对称层合板（T5-7）
数据文档
层数5
层序数厚度mE1（Pa）E2（Pa）v12v21G12（Pa）角度（°）
11.00E-039.60E+102.40E+100.100.401.00E+100.00
21.00E-039.60E+102.40E+100.100.401.00E+1090.00
31.00E-039.60E+102.40E+100.100.401.00E+100.00
41.00E-039.60E+102.40E+100.100.401.00E+1090.00
51.00E-039.60E+102.40E+100.100.401.00E+100.00
（2）正交铺设6层反对称层合板（T5-8）
数据文档
层数6
层序数厚度mE1（Pa）E2（Pa）v12v21G12（Pa）角度（°）
11.00E-039.60E+102.40E+100.100.401.00E+100.00
21.00E-039.60E+102.40E+100.100.401.00E+1090.00
31.00E-039.60E+102.40E+100.100.401.00E+100.00
41.00E-039.60E+102.40E+100.100.401.00E+1090.00
51.00E-039.60E+102.40E+100.100.401.00E+100.00
61.00E-039.60E+102.40E+100.100.401.00E+1090.00
（3）5层对称角铺设层合板（T5-9）
数据文档
层数5
层序数厚度mE1（Pa）E2（Pa）v12v21G12（Pa）角度（°）
11.00E-039.60E+102.40E+100.100.401.00E+1030.00
21.00E-039.60E+102.40E+100.100.401.00E+10-30.00
31.00E-039.60E+102.40E+100.100.401.00E+1030.00
41.00E-039.60E+102.40E+100.100.401.00E+10-30.00
51.00E-039.60E+102.40E+100.100.401.00E+1030.00
（4）6层反对称角铺设层合板（T5-10）
数据文档
层数6
层序数厚度mE1（Pa）E2（Pa）v12v21G12（Pa）角度（°）
11.00E-039.60E+102.40E+100.100.401.00E+1030.00
21.00E-039.60E+102.40E+100.100.401.00E+10-30.00
31.00E-039.60E+102.40E+100.100.401.00E+1030.00
41.00E-039.60E+102.40E+100.100.401.00E+10-30.00
51.00E-039.60E+102.40E+100.100.401.00E+1030.00
61.00E-039.60E+102.40E+100.100.401.00E+10-30.00
（5）我还想验证一个书上的例题，在课本P114。
三层层合板，外层厚度t1，内层10t1，正交铺设比m=0.2,。
玻璃/环氧单层板性能：
E1=5.4E10Pa，E2=1.8E10Pa，v21=0.25，G12=8.8E9Pa。
数据文档
层数3
层序数厚度mE1（Pa）E2（Pa）v12v21G12（Pa）角度（°）
11.00E+005.40E+101.80E+100.0830.2508.80E+090.00
21.00E+015.40E+101.80E+100.0830.2508.80E+0990.00
31.00E+005.40E+101.80E+100.0830.2508.80E+090.00
二、验证Verchery的论文里给出的数值算例。
这里一直到Table5的数据都是从Verchery的论文中截获。
Verchery论文中的18层序列，第（21）式【laminateswithoutbending-extensioncoupling】的排列有两种材料，一种是Boron-Epoxy，一种是Glass-Epoxy。
而且都给出了最终的计算结果Q，A*，D*。
下面是我的Fortran计算数据文档和结果文档。
（1）Boron-Epoxy材料。
（Boron-Epoxy）数据文档
层数18
层序数厚度mE1（Pa）E2（Pa）v12v21G12（Pa）角度（°）
11.00E-032.04E+111.85E+100.0210.2305.59E+090.00
21.00E-032.04E+111.85E+100.0210.2305.59E+090.00
31.00E-032.04E+111.85E+100.0210.2305.59E+0960.00
41.00E-032.04E+111.85E+100.0210.2305.59E+0960.00
51.00E-032.04E+111.85E+100.0210.2305.59E+09-60.00
61.00E-032.04E+111.85E+100.0210.2305.59E+09-60.00
71.00E-032.04E+111.85E+100.0210.2305.59E+09-60.00
81.00E-032.04E+111.85E+100.0210.2305.59E+0960.00
91.00E-032.04E+111.85E+100.0210.2305.59E+090.00
101.00E-032.04E+111.85E+100.0210.2305.59E+09-60.00
111.00E-032.04E+111.85E+100.0210.2305.59E+09-60.00
121.00E-032.04E+111.85E+100.0210.2305.59E+0960.00
131.00E-032.04E+111.85E+100.0210.2305.59E+0960.00
141.00E-032.04E+111.85E+100.0210.2305.59E+090.00
151.00E-032.04E+111.85E+100.0210.2305.59E+090.00
161.00E-032.04E+111.85E+100.0210.2305.59E+090.00
171.00E-032.04E+111.85E+100.0210.2305.59E+0960.00
181.00E-032.04E+111.85E+100.0210.2305.59E+09-60.00
（2）Glass-Epoxy材料。
（Glass-Epoxy）数据文档
层数18
层序数厚度mE1（Pa）E2（Pa）v12v21G12（Pa）角度（°）
11.00E-033.86E+108.27E+090.0560.2604.14E+090.00
21.00E-033.86E+108.27E+090.0560.2604.14E+090.00
31.00E-033.86E+108.27E+090.0560.2604.14E+0960.00
41.00E-033.86E+108.27E+090.0560.2604.14E+0960.00
51.00E-033.86E+108.27E+090.0560.2604.14E+09-60.00
61.00E-033.86E+108.27E+090.0560.2604.14E+09-60.00
71.00E-033.86E+108.27E+090.0560.2604.14E+09-60.00
81.00E-033.86E+108.27E+090.0560.2604.14E+0960.00
91.00E-033.86E+108.27E+090.0560.2604.14E+090.00
101.00E-033.86E+108.27E+090.0560.2604.14E+09-60.00
111.00E-033.86E+108.27E+090.0560.2604.14E+09-60.00
121.00E-033.86E+108.27E+090.0560.2604.14E+0960.00
131.00E-033.86E+108.27E+090.0560.2604.14E+0960.00
141.00E-033.86E+108.27E+090.0560.2604.14E+090.00
151.00E-033.86E+108.27E+090.0560.2604.14E+090.00
161.00E-033.86E+108.27E+090.0560.2604.14E+090.00
171.00E-033.86E+108.27E+090.0560.2604.14E+0960.00
181.00E-033.86E+108.27E+090.0560.2604.14E+09-60.00
（3）当然我也验证了第（22）【laminateswithequalelasticpropertiesinbendingandextension】、（23）【quasi-homogeneouslaminates】的排序，材料是Boron-Epoxy，下面给出计算的结果。
从下面的两个结果表中可以知道，（22）排序的确是C=0，（23）的排序的确是B=0且C=0。
验证成功。
（A）第（22）排序。
（B）第（23）排序。
附件1：
计算所用的程序代码。
PROGRAMCOMPOSITE
!
CodedbyJamesWilson
IMPLICITNONE
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:
charactervariables
WRITE（*,*）"PleaseinserttheINPfilename（a.txtforexample）:
"
READ（*,*）filename1
OPEN（8,file=filename1）!
Opendatafile
!
Readindata
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）
DOI=1,N
READ（8,*）SEQ,T（I）,MC（1:
5）,AL（I）
Q（:
:
I）=0!
Calculatestiffnessofeachlayerfortheprincipalaxis
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!
WorkoutRotMatrix
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
!
WorkoutZ
DOI=1,N
Z（I+1）=Z（I）+T（I）
ENDDO
!
calculatetensorA、BandD
DOK=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（*,*）"PleaseinserttheOUPfilename（b.txtforexample）:
"
READ（*,*）filename2
OPEN（9,file=filename2）
!
Writeinstiffnesstensorforeachsingleply
DOK=1,N
WRITE（9,100）K
100FORMAT（"Thestiffnessofnumber",1X,I2,2X,"plyis:
"）
DOI=1,3
WRITE（9,200）Q（I,:
K）
200FORMAT（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
!
Outputtheactualstiffnessofthelaminate
WRITE（9,"（/）"）;WRITE（9,"（/）"）
WRITE（9,*）"TheACTUALstiffnesstensorofthelaminate:
"
WRITE（9,"（/）"）
WRITE（9,*）"TheextensionstiffnessAequals:
"
DOI=1,3
WRITE（9,200）A（I,1:
3）
ENDDO
WRITE（9,"（/）"）
WRITE（9,*）"ThecouplingstiffnessBequals:
"
DOI=1,3
WRITE（9,200）B（I,1:
3）
ENDDO
WRITE（9,"（/）"）
WRITE（9,*）"ThebendingstiffnessDequals:
"
DOI=1,3
WRITE（9,200）D（I,1:
3）
ENDDO
!
Normalisedtensoroutput
WRITE（9,"（/）"）;WRITE（9,"（/）"）
WRITE（9,*）"TheNORMALISEDstiffnesstensorofthelaminate:
"
WRITE（9,"（/）"）
WRITE（9,*）"TheNORMALISEDextensionstiffnessA*equals:
"
DOI=1,3
WRITE（9,200）A（I,1:
3）/TOTAL_TH
ENDDO
WRITE（9,"（/）"）
WRITE（9,*）"TheNORMALISEDcouplingtensorCequals:
"
DOI=1,3
WRITE（9,200）A（I,1:
3）/TOTAL_TH-12*D（I,1:
3）/TOTAL_TH**3
ENDDO
WRITE（9,"（/）"）
WRITE（9,*）"TheNORMALISEDbendingstiffnessD*equals:
"
DOI=1,3
WRITE（9,200）12*D（I,1:
3）/TOTAL_TH**3
ENDDO
WRITE（*,*）"OUTPUTsuccessfully,pleasepressanykeytoendprogram!
"
READ（*,*）
ENDPROGRAMCOMPOSITE
附2
杨涛同学的MATLAB（GUI）计算程序。
主要程序：
（编了个小界面，程序略长，删掉一些程序自带解释语句，添加了一些对关键语句的解释。
）界面是：
作的一个算例如下：
该算例结果与组内同伴JamesWilson同学基本一致，其余算例结果也基本一致，仅仅在趋近于零时有略微差异，在此不赘于。
后边附上源代码：
functionvarargout=composit_plate（varargin）
gui_Singleton=1;
gui_State=struct（'gui_Name',mfilename,...
'gui_Singleton',gui_Singleton,...
'gui_OpeningFcn',@composit_plate_OpeningFcn,...
'gui_OutputFcn',@composit_plate_OutputFcn,...
'gui_LayoutFcn',[],...
'gui_Callback',[]）;
ifnargin&&ischar（varargin{1}）
gui_State.gui_Callback=str2func（varargin{1}）;
end
ifnargout
[varargout{1:
nargout}]=gui_mainfcn（gui_State,varargin{:
}）;
else
gui_mainfcn（gui_State,varargin{:
}）;
end
functioncomposit_plate_OpeningFcn（hObject,eventdata,handles,varargin）
handles.output=hObject;
guidata（hObject,handles）;
ha=axes（'units','normalized','position',[0011]）;%嵌入坐标，为嵌入背景图片准备
uistack（ha,'down'）%作为背景
II=imread（'武汉大学.jpg'）;%读入图片信息
image（II）
colormaphsv
set（ha,'handlevisib
升级会员 