四边形单元刚度矩阵是有限元分析中的一个重要概念，用于描述四边形单元在受力时的刚度特性。在有限元方法中，连续的求解域被离散为有限个单元的组合，每个单元都有其特定的刚度矩阵。
对于四边形单元，其刚度矩阵是一个方阵，用于将节点位移和节点力联系起来。在弹性力学中，刚度矩阵表示了材料在受力时的抵抗变形的能力。四边形单元的刚度矩阵通常通过对其形状函数和本构关系进行积分得到。
四边形单元的刚度矩阵具有对称性，这是由于材料的本构关系和平衡方程的性质决定的。刚度矩阵的元素反映了节点位移对节点力的影响程度，以及节点间相互作用的强弱。
在有限元分析中，四边形单元刚度矩阵的组装是求解问题的重要步骤之一。通过将所有单元的刚度矩阵按照特定的规则组装到一起，形成整体的刚度矩阵，可以进一步求解出整个结构的位移和应力分布。
四边形单元刚度矩阵的准确性和精度对于有限元分析的结果至关重要。因此，在实际应用中，需要对四边形单元的形状、大小、材料属性等进行合理的选择和描述，以确保分析结果的可靠性。
平面四节点单元的刚度矩阵的fortran编程
parameter(E0=2.0e06,u=0.3,h=2) !参数的选定
real,dimension(2,2):: kii,kij, kji,kjj, kim,kin, kjm,kjn,kmi,kmj,kmm,kmn,kni,knj,knm,knn
real,dimension(3,8)::B !单元总应变矩阵
real,dimension(2,4)::z ！节点的坐标
real,dimension(3,2)::Bi,Bj,Bm,Bn ！单元应变矩阵
real,dimension(8,8)::Kee !单元刚度矩阵
real,dimension(3,3)::D !矩阵 D为 3×3矩阵的定义
real,x,y
data z/xi,yi,xj,xm,ym,xn,yn/
a=z(1,2)-z(1,1)=xj-xi
b=z(2,3)-z(2,2)=ym-yj
Ni=0.25*(1+xi*x/a)*(1+yi*y/b) ！形函数的定义
Nj=0.25*(1+xj*x/a)*(1+yj*y/b)
Nm=0.25*(1+xm*x/a)*(1+ym*y/b)
Nn=0.25*(1+xn*x/a)*(1+yn*y/b)
Bi(1,1)=xi*(1+yi*y/b)/(4*a) ！ Bi矩阵的定义
Bi(1,2)=0
Bi(2,1)=0
Bi(2,2)=yi*(1+xi*x/a)/(4*b)
Bi(3,1)=yi*(1+xi*x/a)/(4*b)
Bi(3,2)=xi*(1+yi*y/b)/(4*b)
Bj(1,1)=xj*(1+yj*y/b)/(4*a)
Bj(1,2)=0
Bj(2,1)=0
Bj(2,2)=yj*(1+xj*x/a)/(4*b)
有限元程序设计
2 平面四边形4结点等参有限单元法
程序设计
1、程序功能及特点
a.该程序采用四边形4节点等参单元，能解决弹性力学的平面应力应变问题。
b.前处理采用网格自动划分技术，自动生成单元及结点信息。
b.能计算受集中力、自重体力、分布面力和静水压力的作用。
c.计算结点的位移和单元中心点的应力分量及其主应力。
d.后处理采取整体应力磨平求得各个结点的应力分量。
e.算例计算结果与ANSYS计算结果比较，并给出误差分析。
f.程序采用Visual Fortran 5.0编制而成。
2、程序流程及图框
启动输入原始数据自动划分网格形成MA，计算N，NH，MX形成整体刚度矩阵K形成荷载列向量RLU分解K=LU回代并求得结点位移输入结点位移计算单元应力及主应力等整体应力磨平结点应力停机
图2-１ 程序流程图 -3- MAINPROGRAMINPUTHUAFENCBANDSKOSTIFFDNXFUN8CONCRBODYRFACERDECOPFOBASTRESSGAUSSSTRESSSUMSSUMSTRSOUTDISTRE
图2-２子程序框图
其中，各子程序的主要功能为：
INPUT――输入原始数据
HUAFEN――自动网格划分，形成COOR(2,NP),X,Y的坐标值与单元信息
CBAND――形成主元素序号指示矩阵MA（*）
SKO――形成整体刚度矩阵[K]
CONCR――计算集中力引起的等效结点荷载{R}e
BODYR――计算自重体力引起的等效结点荷载{R}e
FACER――计算分布面力引起的等效结点荷载{R}e
DECOP――支配方程LU三角分解
FOBA――LU分解直接解法中的回代过程
OUTDISP――输出结点位移分量
STRESS――计算单元应力分量
OUTSTRE――输出单元应力分量
STIF――计算单元刚度矩阵
FDNX――计算形函数对整体坐标的导数TiiyNxN，i1,2,3,4。
/*函数功能：计算单元刚度矩阵
函数入口参数：实型数组nodeCoo[8][3]，存储单元的8个结点的总体自由度编号
实型数组elementK[24][24]，存储单元刚度矩阵
函数返回值：无
*/
void Getnumber(double nodeCoo[8][3],double elementK[24][24])
{
double Gauss[2]={0.57735026919,-0.57735026919}; //高斯点
double WeiG[2]={1.0,1.0}; //权重
int local[3][8]={{1,1,-1,-1,1,1,-1,-1},{-1,1,1,-1,-1,1,1,-1},{-1,-1,-1,-1,1,1,1,1}}; //结点局部坐标
double diffN[3][8]; //插值函数关于局部坐标偏导数
double diffNInverJ[3][8]={0}; //插值函数关于整体坐标偏导数
double DJ;
double J[3][3];
double *MatFirstadd;
double InverJ[3][3]; //存储J的逆矩阵
double B[6][24]={0}; //存储B矩阵
double BInver[24][6]; //存储B的转置
double *MatrixInverFirstadd;
double BInverD[24][6]={0}; //存储B的转置与D的乘积
int i,j,ksi,eta,zeta,m,n;
double Gaussksi,Gausseta,Gausszeta;
double G=E/(2*(1+v)),lam=E*v/((1+v)*(1-2*v)); //已知常数E，v
第２２卷第４期 ２０１１年１２月 苏州市职业大学学报 Ｊｏｕｒｎａｌ ｏｆ Ｓｕｚｈｏｕ Ｖｏｃａｔｉｏｎａｌ Ｕｎｉｖｅｒｓｉｔｙ Ｖｏ１．２２．ＮＯ．４ Ｄｅｃ．，２０１１
基于Ｍａｔｈｅｍａｔｉｃａ的单元刚度矩阵显式解析解推导
张永康 ，鲍四元 ，程宜康 ｂ
（１．苏州市职业大学ａ．机电工程系；ｂ．校长室，江苏苏，Ｊ＇ｌｉ ２１５１０４；２．苏州科技学院土木工程学院，江苏苏州 ２１５０１１）
摘 要：借助Ｍａｔｈｅｍａｔｉｃａ给出了平面四边形单元刚度矩阵的两种显式解析形式．第一种形式采用直 接积分方法，所得的刚度矩阵解析形式比较复杂；第二种结果基于高斯积分方法，所得的刚度矩阵解 析形式相对简单．解析解的推导过程中对复杂的常数项部分用简单的假设变量代替，如对Ｊａｃｏｂｉ行列 式所做的变量代换，刚度矩阵的积分形式得到简化．最后，通过具体算例验证了两种显式解析形式进
行单元刚度矩阵显式求解的有效性和实用性． 关键词：Ｍａｔｈｅｍａｔｉｃａ；刚度矩阵；显式；高斯积分 中图分类号：ＴＰ３９１．７５文献标志码：Ａ 文章编号：１００８—５４７５（２０１１）０４—００３２—０６
Ｅｘｐｌｉｃｉｔ Ａｎａｌｙｔｉｃａｌ Ｓｏｌｕｔｉｏｎｓ ｆｏｒ Ｅｌｅｍｅｎｔ Ｓｔｉｆｆｎｅｓｓ Ｍａｔｒｉｘ Ｂａｓｅｄ ｏｎ Ｍａｔｈｅｍａｔｉｃａ
ＺＨＡＮＧ Ｙｏｎｇ－ｋａｎｇｈ。ＢＡＯ Ｓｉ－ｙｕａｎ ，ＣＨＥＮＧ Ｍ—ｋａｎｇｍ
（１．ａ．Ｄｅｐａｒｔｍｅｎｔ ｏｆ Ｍｅｃｈａｎｏ—Ｅｌｅｃｔｒｏｎｉｃ Ｅｎｇｉｎｅｅｒｉｎｇ，ｂ．Ｐｒｅｓｉｄｅｎｔ Ｏｆｆｉｃｅ，Ｓｕｚｈｏｕ Ｖｏｃａｔｉｏｎａｌ Ｕｎｉｖｅｒｓｉｔｙ，Ｓｕｚｈｏｕ ２１５１０４，Ｃｈｉｎａｌ ２．Ｓｃｈｏｏｌ ｏｆ Ｃｉｖｉｌ Ｅｎｇｉｎｅｅｒｉｎｇ，Ｓｕｚｈｏｕ Ｕｎｉｖｅｒｓｉｔｙ ｏｆ Ｓｃｉｅｎｃｅ ａｎｄ Ｔｅｃｈｎｏｌｏｇｙ，Ｓｕｚｈｏｕ ２１５０１ １，Ｃｈｉｎａ）
