(新)弹性力学平面程序设计——精品.docVIP

弹性力学平面程序设计 概述 序编制说明 本程序用于计算弹性力学平面应力和应变两类问题。这两类问题用类型码ISP区别： ISP=1——平面应力问题 ISP=2——平面应变问题 本程序适用于弹性力学平面问题在给定载荷，如结点荷载、自重作用下的静力计算问题。若有其它非结点载荷作用，则应事先换算成等效结点荷载。 所计算的平面弹性体至少应由三个独立的支承，以保证弹性体的几何性质的不变性。支承方式可以在任何一结点的水平支承和竖直支承。 弹性体由单一的弹性材料组成。 本程序计算方法为有限元位移法。组集整体刚度矩阵时采用直接刚度法和后处理法。 由于整体刚度矩阵K的非零元素分布在以主对角线为中心的斜带形区域，为带形矩阵，根据带形矩阵的特点和对称性，采用半带存储方式。 解方程求解结点位移采用等带宽半带存储高斯消元法。 计算模型及计算方法 1．计算模型、杆端力及杆端位移 采用三角形单元进行分析。由于每个结点有两个位移分量，因此每个单元共有六个位移分量。以这六个结点位移分量作为基本未知量。单元中结点力和结点位移列阵为 2．单元刚度矩阵的建立 在每个单元上，结点力用结点位移表示，可建立单元刚度方程 确定单元刚度矩阵采用步骤： 应用式（8-9）或式（8-10）求弹性矩阵D； 应用式（8-26）求应变转换矩阵B； 依下式求单元刚度矩阵 3．整体刚度矩阵的建立 同杆件结构位移法一样，由单元刚度矩阵按直接刚度法集成。 4．支承条件的引入 首先将结点力向量F换成结点载荷向量P，结构的整体刚度方程变为 K△=P 根据支承条件，把K和P中与支承有杆的元素进行修改。 5．方程的求解及应力计算 采用等带宽半带存贮高斯消元法解方程，可求得结构的结点位移，并存于P中。 为求解单元应力，应从结构位移中取出单元的结点位移，进而求解单元应力。计算式为式（8-27），即 主应力和主平面角的计算公式可参照图9-1中的应力圆得出。 平面应力： 应力圆半径： 最大主应力：STMA=AST＋RST 最小主应力：STMI=AST－RST 主平面角： 弹性力学平面问题计算框图与程序 序标识符说明 平面应力问题和平面应变问题静力分析程序名PASSAP（Plane Stress/Strain Problem Analysis Program）。其主要标识符说明如下： TL（20）——算例标题。实型数组，输入参数。 NJ——结点总数。整形变量，输入参数。 NE——单元数组。整形变量，输入参数。 NZ——约束结点数。整形变量，输入参数。 NPJ——结点荷载数。整形变量，输入参数。 NPJ0——存放结点荷载数NPJ信息，以保证当 NPJ0=0时，亦可以按规定定义数组。 IPS——问题类型码。 ND——半带宽。整形变量。 E——弹性模量。实型变量，输入参数。 PR——|泊松比。实型变量，输入参数。 T——弹形体厚度。实型变量，输入参数。 V——容量。实型变量，输入参数。 LND（NE，3）——单元编号组数。LND（NE，1）、LND（NE，2）、LND（NE，3）分别为单元NE的三个结点编号。整型数组，输入参数。 X（NJ），Y（NJ）——结点坐标数组。X（I）、Y（I）分别为I号结点的x坐标、y坐标。实型数组，输入参数。 JN（NZ，3）——支承结点数组。JZ（I，1）为第I个支承的节点编号。JZ（I，2）、JZ（I，3）分别是第I个支承结点在u，v位移方向的约束信息。该信息为1时，表示有约束；为零时，则无约束。实型数组，输入参数。 PJ（NPJ，3）——结点载荷数组。PJ（I，1）表示第I个支承的结点编号。PJ（I，2）、PJ（I，3）分别表示该结点沿x，y方向作用的结点荷载数值。实型数组，输入参数。 N——自由度数 AKE（6，6）——单元刚度矩阵 B（3，6）——位移-应变转换矩阵（几何矩阵） D（3，3）——弹性矩阵 S（3，6）——应力矩阵，S=DB P（N）——结点荷载数组，存放结点荷载列向量；解方程后，存结点位移。 DE(6)——单元e的结点位移。实型数组。 ST(3)——单元赢利。实型数组。 READ——子程序名，求单元面积。 MKE——子程序名，形成单元刚度矩阵。 MA——子程序名，求单元面积。 MD——子程序名，计算弹性矩阵。 MB——子程序名，计算位移-应变转换矩阵。 MF——子程序名，求结点荷载。 RKR——子程序名，引入支承条件。 SLOV——子程序名，求结点位移并输出。 MADE——子程序名，计算应力并输出。 框图 弹性力学平面问题计算程序总框图如图9-2所示。 子程序MKE 子程序READ 开始 解方程并打