L组Project专业课程设计改.doc

辽宁工程技术大学力学和工程学院
Project汇报书

题 目
工程中板壳强度、刚度计算和稳定 性分析（弹塑性力学）
班 级
工程力学09—1班
姓 名
李慧飞 李世全 李天舒 梁霜
指导老师
李利萍
成 绩

辽宁工程技术大学
力学和工程学院 制

project任务书

Project题目：工程中板壳强度、刚度计算和稳定性分析。
Project关键内容：经过FORTRAN语言对正方形薄板进行了连续体平面问题有限元分析；经过板壳理论基础对负重轮盘进行分析。并对初步板壳理论认识有了深入了解。

学生姓名： 李慧飞 李世全 李天舒 梁霜

指导老师签字：
年 月 日

摘 要
早期薄板断裂问题是依据Kirehhoff假设进行分析。1961年M.L.Williams利用特征值展开方法具体推导了弯曲时直线裂纹尖端应力表示式。1965年件G.C.Sih将Hilbert边值问题推广至板内，求得了直线裂纹尖端应力强度因子。应用复变函数理论和保角映射方法，苏联学者处理了板内含有多种曲线型缺点弯扭问题，这些结果已列入附录中我们也曾用复变函数方法计算唇型和翼型裂纹应力强度因子。1971年W.K.Wilson和Thompson首次用有限元法分析薄板弯曲断裂问题。
和薄板断裂分析相仿，最初薄壳断裂分析亦普遍应用Kirehhoff理论，基于扁壳理论，E.S.Folias研究含裂纹球壳问题应用叠加原理将球壳含裂纹问题分解为无裂纹问题和裂纹面受分布力问题。前者满足基础方程和远离裂纹面处边界条件，以后者满足基础方程并和前者解叠加后在裂纹面上满足内力自由条件经过Fourier积分变换，首先化为对偶积分, 因为核函数求积困难，故深入转换成奇异积分方程，求得含裂纹球壳裂尖应力表示式首项和强度因子。
为克服上述缺点, 多年来已经有较多研究者采取Reissner精化理论进行研究。J.K.Knowles等、R.J.Hartranft和G.C.Sih分别导得Reissner奇异性，在文中得出对称工型问题裂纹尖端应力、应变场展开式，处理了有限尺寸复合型弯曲断裂问题。
综观自1961年至今五十余年时间内, 含裂纹板壳断裂分析已取得很大进展, 不过还有部分问题, 如板壳弹塑性断裂分析, 复合材料层合壳断裂和实用工程计算方法等方面需要深入研究。
工程上大量采取薄板壳型结构，它们在压力作用下，会在内部应力远小于材料屈服极限应力（见材料力学性能）时，忽然产生垂直于压力方向位移而降低承载能力，甚至发生破坏，这种现象称为失稳、皱损或屈曲。板壳失稳是由侧向位移引发，所以失稳属于刚度问题。因为研究板壳失稳问题就形成了板壳稳定性理论。
关键词：板壳；强度；刚度；稳定性

Abstract
Early sheet fracture problem is analyzed according to the Kirehhoff hypothesis. In 1961 M.L.W illiams using eigenvalue method is deduced in detail on bending line of the crack tip stress expression. In 1965 a G.C.S ih will Hilbert boundary value problem is to promote in the plate, get the linear stress intensity factor of the crack tip. Application of complex function theory and conformal mapping method, the Soviet union scholars solved plate contains various curve type defect bending and twisting problem, these results have been listed in the appendix, we also used the methods of complex function calculation lip type and profile crack stress intensity factor. In 1971, for the first time W.K.W ilson and Thompson by finite element method analysis sheet bending fracture problems.
And sheet fracture analysis are similar, initial shell fracture analysis is also widely used Kirehhoff theory, based on the shallow shell theory, E.S.F olias research including crack of spherical shell applications of the superposition principle will be spherical shell with crack problem is decomposed into no crack and the crack surface by distribution force problem. The former meet basic equation and away from the crack surface in the boundary conditions, which meet the basic equation and the former solution stack in the crack surface after meet internal force free condition through the Fourier integral transform, first into the dual integral, due to the nuclear function quadrature difficulties, so further converted into a singular integral equation, get including crack spherical shell of the crack tip stress expression first term and intensity factor.
In order to overcome the above shortcomings, in recent years many researchers using Reissner refine theory research. J.K.K nowles, R.J.H artranft and G.C.S ih respectively to guide Reissner singularity, in the work that symmetrical type problems crack tip stress and strain field expansion, solve the hybrid finite size bending fracture problems.
Synoptic since 1961 twenty years time, including crack plate and shell of fracture analysis has made great progress, but there are still some problems, such as plate and shell elastic-plastic fracture analysis, composite laminated shell fracture and practical engineering calculation method for further study.
Works on a large number of thin lamella-type structure, and they are under pressure, in the internal stress is far less than the yield stress limit of the material (see the mechanical properties of the material), and suddenly generate a displacement perpendicular to the direction of the pressure and reduce the carrying capacity, or even destroyed, a phenomenon known as instability, wrinkled loss or buckling. Plates and Shells instability is caused by the lateral displacement, and thus instability belongs stiffness problem. As the study of plate and shell instability problems on the plate and shell stability theory.
Key words: shell;strength;stiffness;stability

目 录
1．绪论···························································1
2．板壳力学基础知识···············································2
2.1几项假设·······················································2
2.2 板壳应力、应变和应力和应变关系···························2
2.3 薄板内力····················································2
2.4 薄壳内力····················································3
2.5 薄板边界条件················································3
2.6 薄壳边界条件················································3
2.7 板壳应力计算公式············································3
3．工程中实例分析·················································5
3.1 连续体平面问题有限元分析····································5
3.1.1 理论依据和分析··············································5
3.1.2程序原理及实现···············································6
3.1.3 程序原理如框图············································6
3.1.4 算例原始数据和程序分析·····································10
3.1.5算例结果····················································21
3.2 负重轮轮盘强度和刚度板壳理论分析···························23
4.结论··························································27
参考文件 ·······················································28
工程中板壳强度、刚度计算和稳定性分析
1．绪论
有限单元法作为一门课程在现实工程中应用已经十分广泛，我组先对板壳力学基础知识进行简单学习，然后结合简单正方形薄板为例用软件进行计算和分析，最终在现实生活中发觉实例对实例进行计算并和理论进行比较得出结论。本文第一个例子利用计算机，结合FORTRAN语言和有限单元法理论知识，对平面应力问题薄板进行内力值计算和对应结构分析，并能够将由计算机计算实际结果和由用弹性力学理论知识计算出来解答相比较。第二个例子对负重轮轮盘进行了模型简化，利用板壳理论建立并求解了变形微分方程，得出轮盘内各点应变和应力分量，并深入对轮盘进行强度、刚度校核和设计，其结果和有限元分析结果靠近。从而有效对板壳强度，刚度稳定性进行了一定研究。

2. 板壳力学基础知识
2.1 几项假设
1.板壳是均匀、连续，而且是各向同性；
2.板壳是线弹性；
3.板壳变形是微小；
4.直法线假设，即认为板壳变形前垂直于中面法线线段在变形后仍保持为直线，并垂直于变形后中面，且其长度不变。
5.法向应力很小，能够忽略；
6.板中面没有变形。
2.2 板壳应力、应变和应力和应变关系
薄板壳内任一点沿z方向位移wA和坐标z无关，仅是坐标x、y函数，横向剪应变yz和zx应为零，
几何方程
（1-1）

物理方程

（1-2） 或 （1-3）

2.3 薄板内力

（1-4）

2.4 薄壳内力

（1-5）

2.5 薄板边界条件
1.简支边 ， （1-6）
2.固定边 ， （1-7）
3.自由边 ， （1-8）
2.6 薄壳边界条件
1.简支边 ， ，， （1-9）
2.固定边 ， ，， （1-10）
3.自由边 ，，， （1-11）
2.7 板壳应力计算公式
1.薄板
（1-12）
其中惯性矩，静面矩
2.薄壳
（1-13）
注意，在板壳弯曲问题中，数值上最大是法向应力、和切向应力，所以是关键应力，横向剪应力、数值较小，是次要应力，通常说来，无须对它们进行计算。

3．工程中实例分析
3.1 连续体平面问题有限元分析
正方形薄板四面受均匀载荷作用，该结构在边界上受正向分布压力，P=1kN/m，同时在沿对角线y轴上受一对集中压力，载荷为2kN，板厚t=1，泊松比ν=0，见下图：
Y

2kN
P=1kN/m
P
X
o

P P

2kN

图2.1正方形薄板四面受均匀载荷作用
Figure 2.1 square sheet around by the uniform load role
3.1.1 理论依据和分析
此问题，为弹性力学里平面应力问题，在板内部，四处全部有σz=0,τyz=0,τxz=0;σx=f1(x,y),σy=f2(x,y),τxy=f3(x,y)，应力含有这种性质问题，称为平面应力问题。弹性薄板在工程中应用很广泛，对于部分简单情况，如等厚、单跨、无大孔口，外形规则（如矩形，圆形等）薄板，已经有部分解答和表格可资利用盘。因为连续平板连续性，仅需要取其在第一象限四分之一部分研究计算，然后做出部分辅助线将平板分成若干部分，在为每个部分选择分子单元，采取此模型化为4个全等三角形单元，利用其对称性，四分之一边界约束，荷载可等效图所表示。

1kN/m

图2.2 三角形单元内力
Figure 2.2 triangle unit internal force
3.1.2 程序原理及实现
用FORTRAN程序实现。有节点信息文件NODE.IN和单元信息文件ELEMENT.IN,经过计算分析后输出一个通常性文件DATA.OUT。模型基础信息由文件BASIC.IN生成。该程序特点以下：问题类型：可用于计算弹性力学平面应力问题和平面应变问题。单元类型：采取常应变三角形单元。位移模式：用线性位移模式。载荷类型：节点载荷，非节点载荷应先转换为等效节点载荷。材料性质：弹性体由单一均匀材料组成。约束方法：为“0”位移固定约束，为确保无刚体位移，弹性体最少应有三个自由度独立约束。方程求解：针对半带宽刚度方程GUASS消去法。输出文件：由手工生成节点信息文件NODE.IN，和单元信息文件ELEMENT,IN。结果文件：输出通常结果文件DATA.OUT。
3.1.3 程序原理如框图
开始

输入数据（子程序READ_IN）
BASIC.IN(基础信息文件)
NODE.IN(节点信息文件)
ELEMEENT.IN(单元信息文件)
↓
形成单元刚度矩阵（子程序FORM_KE）
↓
以半带存放方法形成整体刚度矩阵
（BAND K）
↓
形成节点载荷向量（子程序FORM_P）
↓
处理边界条件（子程序DO_BC）
↓
求解方程取得节点位移（子程序SOLVE）
↓
计算单元及节点应力（子程序）
↓
输出文件
DATA.OUT
↓
结束
说明：
（1）关键变量：
ID：问题类型码：，ID=1时平面应力问题，ID=2时平面应变问题
N_NODE： 节点个数
N_LOAD： 节点荷载个数
N_DOF： 自由度，N_DOF=N_NODE*2（平面问题）
N_ELE： 单元个数
N_BAND： 矩阵半带宽
N_BC： 有约束节点个数
PE： 弹性模量
PR： 泊松比
PT： 厚度
LJK_ELE(I，3)：单元节点编号数组，LJK_ELE(I,1),LJK_ELE(I,2),LJK_ELE(I,3)分别放单元I三个节点整体编号。X（N_NODE）,Y(N_NODE):节点坐标数组，X(I),Y(I)分别存放节点Ix,y坐标值。
LJK_U(N_BC,3)：节点载荷数组，P_LJK(I,1)表示第I个作用有节点载荷节点编号，P_LJK(I,2),P_LJK(I,3)分别为该节点沿x,y方向节点载荷数值。
AK(N_DOF,N_BAND): 整体刚度矩阵
AKE(6,6): 单元刚度矩阵
BB(3,6): 位移应变转换矩阵（三节点几何矩阵）
DD(3,3): 弹性矩阵
SS(3,6): 应力矩阵
RESULT_N(N_NOF):节点荷载数组，存放节点荷载向量，解方程后该矩阵存放节点位移
DISP_E(6): 单元节点位移向量
STS_ELE(N_ELE,3): 单元应力分量
STS_ND(N_NODE,3): 节点应力分量
（2） 子程序说明：
READ_IN: 读入数据 BAND_K: 形成半带宽整体刚度矩阵
FORM_FE: 计算单元刚度矩阵 FORM_P: 计算节点载荷
CAL_AREA: 计算单元面积 DO_BC: 处理边界条件
CLA_DD: 计算单元弹性矩阵 SOLVE: 计算节点位移
CLA_BB: 计算单元位移......应变关系矩阵
CAL_STS: 计算单元和节点应力
（3） 文件处理：源程序文件：chengxu,for 程序读入数据文件：
BASIC.IN,NODE.IN,ELEMENT.IN（需要手工生成）
程序输出数据文件：
DATA.OUT
(4) 数据文件格式：
需读入模型、基础信息文件BASIC.IN格式以下表：
栏目
格式说明
实际需输入数据
基础模型数据
第1行，每两个数之间用“，”号
隔开
问题类型，单元个数，节点个数，有约束节点数，有载荷节点数
材料性质
第2行，每两个数之间用“，”号
隔开
弹性模量，泊松比，单元厚度
节点约束信息
在材料性质输入行以后另起行，每两个数之间用“，”号隔开
LJK_U(N_BC,3)
位移约束节点编号，该节点x方向约束代码，该节点y方向代码
节点载荷信息
在节点约束信息输入行以后另起行，每两个数之间用“，”号隔开
P_IJK(N_LOAD,3)
载荷作用节点编号，该节点x方向载荷
该节点y方向载荷，……
需读入节点信息文件NODE.IN格式以下表
栏目
格式说明
实际需输入数据
节点信息
每行为一个节点信息（每行三个数，每两个数之间用空格或“，”分开
LJK_U(N_BC,3)
节点号，该节点x坐标，该节点y方向坐标
需读入单元信息文件ELEMENT.IN格式以下表
栏目
格式说明
实际需输入数据
单元信息
每行为一个单元信息（每行有14个整型数4个为单元节点编号，对于3节点编号，第4个节点编号和第3个节点编号相同，后10个数无用，可输入“0”，每两个整型数之间用最少一个空格分开
NE_ANSYS(N_ELE,14)
单元节点号1（空格）
单元节点号2（空格）
单元节点号3（空格）
单元节点号4（空格）
0（空格）0（空格）0（空格）0（空格）0（空格）0（空格）0（空格）0（空格）0（空格）0
输出结果文件DATA.OUT格式以下表
栏目
实际需输入数据
节点位移
I RESULT_N(2*I_1) RESULT_N(2*I)
节点号 x方向位移 y方向位移
单元应力三个分量
IE STE_ELE(IE,1) STE_ELE(IE,2) STE_ELE(IE,3)
单元号 x方向应力 y方向应力 剪切应力
节点应力三个分量
I STS_ND(I,1) STS_ND(I,2) STS_ND(I,3)
节点号 x方向应力 y方向应力 剪切应力
3.1.4 算例原始数据和程序分析
（1） 模型基础信息文件BASIC.IN数据为
1,4,6,5,3
1,0,1.
1，1,0,2,1,0,4,1,1,5,0,1,6,0,1
1，-0.5，-1.5,3，-1，-1,6，-0.5，-0.5
（2） 手工准备节点信息文件 NODE.IN数据为
1 0.0 2.0
2 0.0 1.0
3 1.0 1.0
4 0. 0.
5 1.0 0.
6 2.0 0.
（3） 手工准备单元信息文件ELEMENT.IN数据为
1 2 3 3 0 0 0 0 1 1 1 1 0 1
2 4 5 5 0 0 0 0 1 1 1 1 0 2
5 3 2 2 0 0 0 0 1 1 1 1 0 3
3 5 6 6 0 0 0 0 1 1 1 1 0 4
（4） 源程序文件chengxu,for为：
PROGRAM FEM2D
DIMENSION IJK_ELE(500,3),X(500),Y(500),IJK_U(50,3),P_IJK(50,3),
&RESULT_N(500),AK(500,100)
DIMENSION STS_ELE(500,100),STS_ND(500,3)
OPEN(4,FILE='BASIC.IN')
OPEN(5,FILE='NODE.IN')
OPEN(6,FILE='ELEMENT.IN')
OPEN(8,FILE='DATA.OUT')
OPEN(9,FILE='FOR_POST.DAT')
READ(4,*)ID,N_ELE,N_NODE,N_BC,N_LOAD
IF(ID.EQ.1)WRITE(8,20)
IF(ID.EQ.2)WRITE(8,25)
20 FORMAT(/5X,'=========PLANE STRESS PROBLEM=========')
25 FORMAT(/5X,'=========PLANE STRAIN
PROBLEM=========')
CALL READ_IN(ID,N_ELE,N_NODE,N_BC,N_BAND,N_LOAD,PE,PR,PT,
& IJK_ELE,X,Y,IJK_U,P_IJK)
CALL BAND_K(N_DOF,N_BAND,N_ELE,IE,N_NODE,
& IJK_ELE,X,Y,PE,PR,PT,AK)
CALL FROM_P(N_ELE,N_NODE,N_LOAD,N_DOF,IJK_ELE,X,Y,P_IJK,
& RESULT_N)
CALL DO_BC(N_BC,N_BAND,N_DOF,IJK_U,AK,RESULT_N)
CALL SOLVE(N_NODE,N_DOF,N_BAND,AK,RESULT_N)
CALL CAL_STS(N_ELE,N_NODE,N_DOF,PE,PR,IJK_ELE,X,Y,RESULT_N
& STS_ELE,STS_ND)
C to putout a data file
WRITE(9,70)REAL(N_NODE),REAL(N_ELE)
70 FORMAT(2f9.4)
WRITE(9,71)(X(I),Y(I),RESULT_N(2*I-1),RESULT_N(2*I),
&STS_ND(I,1),STS_ND(I,2),STS_ND(I,3),I=1,N_NODE)
71 FORMAT(7F9.4)
WRITE(9,72)(REAL(IJK_ELE(I,1)),REAL(IJK_ELE(I,2)).
&REAL(IJK_ELE(I,3))REAL(IJK_ELE(I,3)),
&STS_ELE(I,1),STS_ELE(I,2),STS_ELE(I,3),I=1,N_ELE)
72 FORMAT(7F9.4)
C
CLOSE(4)
CLOSE(5)
CLOSE(8)
CLOSE(9)
END
C
C to get the original data in order to model the problem
SUBROUTINE READ_IN(ID,N_ELE,N_NODE,N_BC,N_BAND,N_LOAD,PE,PR,
&PT,IJK_ELE,X,Y,IJK_U,P_IJK)
DIMENSION IJK_ELE(500,3),X(N_NODE),Y(N_NODE),IJK_U(N_BC,3),
& P_IJK(N_LOAD,3),NE,ANSYS(N_ELE,14)
REAL ND_ANSYS(N_NODE,3)
READ(4,*)PE,PR,PT
READ(4,*)((IJK_U(I,J),J=1,3),I=1,N_BC)
READ(5,*)((ND_ANSYS(I,J),J=1,3),I=1,N_ELE)
DO 10 I=1,N_NODE
X(I)=ND_ANSYS(I,2)
Y(I)=ND_ANSYS(I,3)
10 CONTINUE
DO 11 I=1,N_ELE
DO 11 J=1,3
IJK_ELE(I,J)=NE_ANSYS(I,J)
11 CONTINUE
N_BAND=0
DO 20 IE=1,N_ELE
DO 20 I=1,3
DO 20 J=1,3
IW=IABS(IJKK_ELE(IE,I)-IJK_ELE(IE,J))
IF(N_BAND.LT.IW)N_BAND=IW
20 CONTINUE
N=BAND=(N_BAND+1)*2
IF(ID.EQ.1) THEN
ELSE
PE=PE/(1.0-PR*PR)
PR=PR/(1.0-PR)
END IF
RETURN
END
C
C to form the stiffness matrix of element
SUBROUTINE FORM_KE(IE,N_NODE,N_ELE,IJK_ELE,X,Y,PE,PR,PT,AKE)
DIMENSION IJK_ELE(500,3),X(N_NODE),Y(N_NODE),BB(3,6),DD(3,3),
& AKE(6,6),SS(6,6)
CALL CAL_DD(PE,PR,DD)
CALL CAL_BB(IE,N_NODE,N_ELE,IJK_ELE,X,Y,AE,BB)
DO 10 I=1,3
DO 10 J=1,6
SS(I,J)=SS(I,J)+DD(I,K)*BB(K,J)
DO 20 I=1,6
DO 20 J=1,6
AKE(I,J)=0.0
DO 20 K=1,3
20 AKE(I,J)=AKE(I,J)+SS(K,I)*BB(K,J)*AE*PT
RETURN
END
C to form banded global stiffness matrix
SUBROUTINE BAND_K(N_DOF,N_BAND,N_ELE,IE,N_NODE,IJK_ELE,X,Y,PE,
& PR,PT,AK)
DIMENSION
IJK_ELE(500,3),X(N_NODE),Y(N_NODE),AKE(6,6),AK(500,100)
N_DOF=2*N_NODE
DO 40 I=1,N_DOF
DO 40 J=1,N_BAND
40 AK(I,J)=0
DO 50 IE=1,N_ELE
CALL FORM_KE(IE,N_NODE,N_ELE,IJK_ELE,X,Y,PE,PR,PT,AKE)
DO 50 I=1,3
DO 50 II=1,2
IH=2*(I=1)+II
DO 50 J=1,3
DO 50 JJ=1,2
IL=2*(J-1)+JJ
IZL=2*(IJK_ELE(IE,J)-1)+JJ
IDL=IZL-IDH+1
IF(IDL.LE.0) THEN
ELSE
AK(IDH,IDL)=AK(IDH,IDL)+AKE(IH,IL)
END IF
50 CONTINUE
RETURN
END
C
C to calculate the area of element
SUBROUTINE CAL_AREA(IE,N_NODE,IJK_ELE,X,Y,AE)
DIMENSION IJK_ELE(500,3),X(N_NODE),Y(N_NODE)
I=IJK_ELE(IE,1)
J=IJK_ELE(IE,2)
K=IJK_ELE(IE,3)
XIJ=X(