MatLab在有限元刚度矩阵推导中的应用

第２６卷第２期重庆交通学院学报ｖ０１．２６Ｎｏ．２２００７年４月ＪＯＵＲＮＡＬＯＦＣＨＯＮＧＱＩＮＧ：姒０７ｒＯＮＧＵＮⅣＥＲＳｒｒＹＡｐｒ．．２００７·＿＿＿＿Ｉ＿＿＿＿目一ｒｒ自目ｌ＿＝＝日日一

ＭａｔＬａｂ在有限元刚度矩阵推导中的应用

周水兴１’２，陈山林１

（１．重庆大学，重庆４０００４５；２．重庆交通大学土木建筑学院，重庆４０００７４）摘要：ＭａｔＬａｂ具有强大的符号运算功能．以平面梁单元刚度矩阵为例，详细介绍ＭａｔＬａｂ在单元刚度矩阵推导中的具体应用，编写了平面梁单元刚度符号运算程序，运算结果与手工推导完全一致，该方法可进一步推广到其它单元甚至到更复杂非线性单元刚度矩阵推导中．

关键词：ＭａｔＬａｂ；符号运算；刚度矩阵；推导

中图分类号：０２４２．２１文献标识码：Ａ文章编号：１００１．７１６Ｘ（２００７）０２．００２９．０３

ＭａｔＬａｂＡｐｐｌｉｃａｔｉｏｎｔｏｔｈｅＭａｔｒｉｘＤｅｄｕｃｔｉｏｎｏｆＦｉｎｉｔｅＥｌｅｍｅｎｔＳｔｉｆｆｎｅｓｓ

ＺＨＯＵＳｈｕｉ—ｘｉｎ９１．．。ＣＨＥＮＳｈａｎ．１ｉｎｌ

（１．ＣｈｏｎｇｑｉｎｇＵｎｉｖｅｒｓｉｔｙ，Ｃｈｏｎｇｑｉｎｇ４０００４５，Ｃｈｉｎａ；

２．ＳｃｈｏｏｌｏｆＣｉｖｉｌＥｎｇｉｎｅｅｒｉｎｇ＆Ａｒｃｈｉｔｅｃｔｕｒｅ，ＣｈｏｎｇｑｉｎｇＪｉａｏｔｏｎｇＵｎｉｖｅｒｓｉｔｙ，Ｃｈｏｎｇｑｉｎｇ４０００７４，Ｃｈｉｎａ）Ａｂｓｔｒａｃｔ：ＭａｔＬａｂｈａｓｓｔｒｏｎｇｓｙｍｂｏｌｉｃｏｐｅｒａｔｉｏｎｆｕｎｃｔｉｏｎ．Ｗｉｔｈｔｈｅｅｘａｍｐｌｅｏｆｓｔｉｆｆｎｅｓｓｍａｔｒｉｘｏｆｐｌａｎｅ

ｂｅａｍｅｌｅｍｅｎｔ，ｔｈｅａｐｐｌｉｃａｔｉｏｎｏｆｄｅｄｕｃｉｎｇｅｌｅｍｅｎｔｓｔｉｆｆｎｅｓｓｍａｔｒｉｘｉｓｉｎｔｒｏｄｕｃｅｄｉｎｄｅｔａｉｌｆｏｒｔｈｅｆｉｒｓｔｔｉｍｅ，ａｎｄｔｈｅｓｙｍｂｏｌｉｃｏｐｅｒａｔｉｏｎｐｒｏｃｅｄｕｒｅ

ｏｆｐｌａｎｅｂｅａｍｅｌｅｍｅｎｔｉｓｄｅｖｄｏｐｅｄ，ａｎｄｔｈｅｏｐｅｒａｔｉｏｎｒｅｓｕｌｔｉｓ

ｃｏｍｐｌｅｔｅｌｙｉｎ

ａｃｃｏｒｄａｎｃｅｗｉｔｈｔｈｅｈａｎｄｉｃｒａｆｔｄｅｄｕｃｔｉｏｎ．Ｔｈｉｓｍｅｔｈｏｄｃａｎｆｕｒｔｈｅｒｅｘｐａｎｄｔｈｅａｐｐｌｉｃａｔｉｏｎｏｆｄｅｄｕｃｉｎｇｍａｔｒｉｘｔｏｏｔｈｅｒｅｌｅｍｅｎｔｓｍｏｒｅｃｏｍｐｌｉｃａｔｅｄｎｏｎｌｉｎｅａｒｅｌｅｍｅｎｔ．ＫｅｙｗｏｒｄｓｉＭａｔｌ．ａｂ，ｓｙｍｂｏｌｉｃｏｐｅｒａｔｉｏｎ，ｅｌｅｍｅｎｔｓｔｉｆｆｎｅｓｓｍａｔｒｉｘ，ｄｅｄｕｃｅ

ＭａｔＬａｂ是美国ＭａｔｈＷｏｒｋ公司开发的用于数值计算、概念设计、算法研究、建模仿真、实时实现理想的集成环境，因其完整的专业体系和强大的运算功能，已广泛应用于工业、电子、信号处理、控制、医疗、建筑、教学等各个领域…．

ＭａｔＬａｂ本意是矩阵实验室（ＭａｔｒｉｘＬａｂｏｒａｔｏｒｙ），以向量和矩阵运算为基础，强大的数值计算功能是其最具代表性的特点．ＭａｔＬａｂ的数学计算包括数值计算和符号计算，前者主要用于矩阵、向量的各种数值计算；后者主要是计算式中带有符号变量、表达式的运算．该运算符由符号数学工具箱提供支持，有复合、简化、微分、积分及求解代数方程和微分方程的工具．ＭａｔＬａｂ符号数学工具箱中的工具建立在功能强大的Ｍａｐｌｅ软件基础上．

有限元是近代数值计算最有效方法之一．有限元法的基础是单元刚度矩阵的推导，目前，刚度矩阵推导已有一个相对固定的模式，而烦琐、复杂的矩阵运算、微分、积分是刚度矩阵推导中的主要内容．通常，这种矩阵运算是由手工来完成的，工作量大，而且极易出错．利用ＭａｔＬａｂ丰富的符号运算功能，完成刚度矩阵推导过程中的符号运算，不但速度快，而且准确性高．

Ｍａｆｌａｂ的命令和函数相当多，本文仅就涉及到ＭａｔＬａｂ符号运算作一些简要介绍，然后结合平面梁单元刚度矩阵的推导，用ＭａｔＬａｂ完成该单元刚度矩阵的符号运算，结果完全一致．

１ＭａｔＬａｂ符号运算的一般规定

１．１符号标变量和表达式的创建

收稿日期：２００６－０２－０７；修订日期：２００６－０２－２３

基金项目：重庆市教委项目（ＫＪ０５０４０７）；交通部西部交通建设科技项目（２００３３１８８１４２０２４）

作者简介：周水兴（１９６７－），男，浙江嘉兴人，教授，主要从事桥梁设计理论研究．ｅ－ｍａｌｌ：ｓｈｕｉｘｉｎｇｚｈｏｕ＠ｃｑｕｃ．ｅｄｕ．ｃｎ

３０重庆交通学院学报第２６卷

ＭａｔＬａｂ中符号标量、变量和表达式是通过命令ｓｙｍ、ｓｙｍｓ来创建的．如：

＞＞ｘ＝ｓｙｍ（’ｘ’）；

＞＞ａｂｃ＝ｓｙｍ（’ａｌｐｈａ’）

＞＞ｆｕｎ＝ｓｙｍ（’ａ宰ｘ＇２＋ｂ＇ｌｃＸ＋ｃ’）

这里，＞＞是ＭａｔＬａｂ命令提示符，上述命令中前两个创建符号标量ｘ，ａｌｐｈａ，最后一个创建符号表达式ａ木ｘ２＋ｂ半ｘ＋ｃ．对符号表达式，还需指定表达式中的各个标量，如ａ，ｂ，Ｃ，Ｘ．ＭａｔＬａｂ中，除了上面介绍的方法外，还提供了一种简洁的方法，ｓｙｍｓ命令，如：

＞＞ｓｙｍｓａｂＣＸ

１．２符号矩阵

在ＭａｔＬａｂ中，符号矩阵的元素为符号表达式，可用函数ｓｙｍ产生，如：

＞＞Ｓ＝ｓｙｍ（’［ｘ，Ｙ，ｚ；ａ，ｂ，ｃ；ｕ，ｖ，Ｗ］’）

该命令等价于：

＞＞ｓｙｍｓｘＹｚａｂｃｕＶＷ

＞＞Ｓ＝［ｘ，Ｙ，ｚ；ａ，ｂ，ｃ；ｕ，ｖ，ｗ］

两者效果是一致的．对数值矩阵，也可以用ｓｙｍ把数值矩阵转换为符号形式．在转换过程中，函数ｓｙｍ尽可能地把数值矩阵中的非整元素指定为小的整数之比，以有理方式表示，如果元素是无理数，则在符号形式中ｓｙｍ将用符号浮点数表示元素．１．３符号矩阵的化简和格式化

符号运算后，符号表达式有时非常复杂，为此，ＭａｔＬａｂ提供了能够化简复杂符号表达式的诸多命令，有ｐｒｅｔｔｙ（用于将符号表达式以常用方式显示）、ｃｏｌｌｅｃｔ（用于将表达式中同类项合并、合并后的多项式以变量的幂次方按大小依次排列）、ｅｘｐａｎｄ（用于展开符号表达式中的各项子式）、ｆａｃｔｏｒ（用于符号分解）、ｓｕｂｓ（较长子表达式的置换）、ｓｉｍｐｌｉｆｙ（用于符号表达式中进行等式的恒等替换）、ｓｉｍｐｌｅ等．１．４符号微积分

Ｍａｔｌａｂ符号工具箱中提供了多种符号表达式的微分ｄｉｆｆ和积分ｉｎｔ函数形式，有：

ｄｉｆｆ（ｆ），ｄｉｆｆ（ｆ，’ｔ’），ｄｉｌｌ（ｆ，ｎ），ｄｉｆｆ（ｆ，’ｔ’ｎ），ｉｎｔ（ｆ），ｉｎｔ（ｆ，’ｔ’），ｉｎｔ（ｆ，ａ，ｂ），ｉｎｔ（ｆ，’ｔ’，ａ，ｂ）和ｉｎｔ（ｆ，’ｍ’，’ｎ’）．有关微分ｄｉｆｆ和积分ｉｎｔ的详细内容可参考ＭａｔＬａｂ教程，这里不再赘述．

２平面梁单元刚度矩阵

为便于读者更好的理解ＭａｔＬａｂ在单元刚度矩阵推导中的应用，文中简要列出了平面梁单元刚度矩阵的推导过程．

２．１平面梁单元位移模式‘２１

对平面梁单元，单元的位移模式为

Ⅱ＝ａｏ＋ａＩ石

（１）秽＝ｂｏ＋ｂｌ省＋６２菇２＋ｂｓｘ３（２）

单元节点位移列向量：

…。＝［ｕ；吼０ｉ叶吩岛］７

利用单元ｉｊ的边界条件，解出上述系数，有

ｕ＝［（１一ｘ／１）００ｘ／ｌ００］｛６｝。＝札（髫）｛艿｝。（３）ｔＪ＝［０（１—３ｘ２／１２＋２戈３／Ｚ３）（％一２ｘ２／ｌ＋菇３／２２）０（３ｘ２／１２—２ｘ３／１３）（一搿２／Ｚ＋

戈３／Ｚ２）］｛艿｝。＝［＾Ｌ（菇）］｛６｝。（４）２．２用结点位移表示应变和应力

忽略剪切影响，梁单元受到拉压和弯曲变形后的线应变为：

…＝ｈｄｕ／ｄ…ｘ：）＝一■…。＝

［Ｂ］｛艿｝。（５）式中，

［日］：∥引，１

Ｌ一∥：（髫）Ｊ

（６）［Ｂ］为应变矩阵．由胡克定律，单元中的应力表达式为：

｛矿｝＝［Ｄ］｛占｝＝［Ｄ］［Ｂ］｛艿｝‘（７）

２．３梁单元刚度矩阵

根据虚位移原理，梁单元刚度矩阵为

［ｋ］＝『ＪＪ［Ｂ］’［Ｄ］［Ｂ］ｄｙ（８）３符号矩阵的运算

３．１应变矩阵［Ｂ］的计算

应变矩阵中，有两个形函数表达式［札（菇）］、［见（石）］的微分，在Ｍａｔｌａｂ中定义应变矩阵和微分的程序段如下：

Ｎｕ＝［１一ｘ／Ｌ００ｘ／Ｌ００］％定义形函数矩阵Ｎｕ

Ｎｖ＝［０１—３，ｌｃｘ＇２／Ｌ＇２＋２，ｌｃｘ＇３／Ｌ３戈一

２拳ｘ＇２／Ｌ＋ｘ＾３／Ｌ＇２０３’ｌｃｘ＇２／Ｌ＂２—２木ｘ＾３／Ｌ＾３一ｘ＇２／Ｌ＋ｘ＾３／Ｌ＇。２］

Ｂ＝［ｄｉｆｆ（Ｎｕ，’菇’，１）；一Ｙ术ｄｉｆｆ（Ｎｖ，’戈’，２）］％两个形函数分别进行变量戈一阶微分，然后形成