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

第２６卷第２期重庆交通学院学报ｖ０１．２６Ｎｏ．２２００７年４月ＪＯＵＲＮＡＬＯＦＣＨＯＮＧＱＩＮＧ：姒０７ｒＯＮＧＵＮⅣＥＲＳｒｒＹＡｐｒ．．２００７·＿＿＿＿Ｉ＿＿＿＿目一ｒｒ自目ｌ＿＝＝日日一ＭａｔＬａｂ在有限元刚度矩阵推导中的应用周水兴１’２，陈山林１（１．重庆大学，重庆４０００４５；２．重庆交通大学土木建筑学院，重庆４０００７４）摘要：ＭａｔＬａｂ具有强大的符号运算功能．以平面梁单元刚度矩阵为例，详细介绍ＭａｔＬａｂ在单元刚度矩阵推导中的具体应用，编写了平面梁单元刚度符号运算程序，运算结果与手工推导完全一致，该方法可进一步推广到其它单元甚至到更复杂非线性单元刚度矩阵推导中．关键词：ＭａｔＬａｂ；符号运算；刚度矩阵；推导中图分类号：０２４２．２１文献标识码：Ａ文章编号：１００１．７１６Ｘ（２００７）０２．００２９．０３ＭａｔＬａｂＡｐｐｌｉｃａｔｉｏｎｔｏｔｈｅＭａｔｒｉｘＤｅｄｕｃｔｉｏｎｏｆＦｉｎｉｔｅＥｌｅｍｅｎｔＳｔｉｆｆｎｅｓｓＺＨＯＵＳｈｕｉ—ｘｉｎ９１．．。

ＣＨＥＮＳｈａｎ．１ｉｎｌ（１．ＣｈｏｎｇｑｉｎｇＵｎｉｖｅｒｓｉｔｙ，Ｃｈｏｎｇｑｉｎｇ４０００４５，Ｃｈｉｎａ；２．ＳｃｈｏｏｌｏｆＣｉｖｉｌＥｎｇｉｎｅｅｒｉｎｇ＆Ａｒｃｈｉｔｅｃｔｕｒｅ，ＣｈｏｎｇｑｉｎｇＪｉａｏｔｏｎｇＵｎｉｖｅｒｓｉｔｙ，Ｃｈｏｎｇｑｉｎｇ４０００７４，Ｃｈｉｎａ）Ａｂｓｔｒａｃｔ：ＭａｔＬａｂｈａｓｓｔｒｏｎｇｓｙｍｂｏｌｉｃｏｐｅｒａｔｉｏｎｆｕｎｃｔｉｏｎ．Ｗｉｔｈｔｈｅｅｘａｍｐｌｅｏｆｓｔｉｆｆｎｅｓｓｍａｔｒｉｘｏｆｐｌａｎｅｂｅａｍｅｌｅｍｅｎｔ，ｔｈｅａｐｐｌｉｃａｔｉｏｎｏｆｄｅｄｕｃｉｎｇｅｌｅｍｅｎｔｓｔｉｆｆｎｅｓｓｍａｔｒｉｘｉｓｉｎｔｒｏｄｕｃｅｄｉｎｄｅｔａｉｌｆｏｒｔｈｅｆｉｒｓｔｔｉｍｅ，ａｎｄｔｈｅｓｙｍｂｏｌｉｃｏｐｅｒａｔｉｏｎｐｒｏｃｅｄｕｒｅｏｆｐｌａｎｅｂｅａｍｅｌｅｍｅｎｔｉｓｄｅｖｄｏｐｅｄ，ａｎｄｔｈｅｏｐｅｒａｔｉｏｎｒｅｓｕｌｔｉｓｃｏｍｐｌｅｔｅｌｙｉｎａｃｃｏｒｄａｎｃｅｗｉｔｈｔｈｅｈａｎｄｉｃｒａｆｔｄｅｄｕｃｔｉｏｎ．Ｔｈｉｓｍｅｔｈｏｄｃａｎｆｕｒｔｈｅｒｅｘｐａｎｄｔｈｅａｐｐｌｉｃａｔｉｏｎｏｆｄｅｄｕｃｉｎｇｍａｔｒｉｘｔｏｏｔｈｅｒｅｌｅｍｅｎｔｓｍｏｒｅｃｏｍｐｌｉｃａｔｅｄｎｏｎｌｉｎｅａｒｅｌｅｍｅｎｔ．ＫｅｙｗｏｒｄｓｉＭａｔｌ．ａｂ，ｓｙｍｂｏｌｉｃｏｐｅｒａｔｉｏｎ，ｅｌｅｍｅｎｔｓｔｉｆｆｎｅｓｓｍａｔｒｉｘ，ｄｅｄｕｃｅＭａｔＬａｂ是美国ＭａｔｈＷｏｒｋ公司开发的用于数值计算、概念设计、算法研究、建模仿真、实时实现理想的集成环境，因其完整的专业体系和强大的运算功能，已广泛应用于工业、电子、信号处理、控制、医疗、建筑、教学等各个领域…．ＭａｔＬａｂ本意是矩阵实验室（ＭａｔｒｉｘＬａｂｏｒａｔｏｒｙ），以向量和矩阵运算为基础，强大的数值计算功能是其最具代表性的特点．ＭａｔＬａｂ的数学计算包括数值计算和符号计算，前者主要用于矩阵、向量的各种数值计算；后者主要是计算式中带有符号变量、表达式的运算．该运算符由符号数学工具箱提供支持，有复合、简化、微分、积分及求解代数方程和微分方程的工具．ＭａｔＬａｂ符号数学工具箱中的工具建立在功能强大的Ｍａｐｌｅ软件基础上．有限元是近代数值计算最有效方法之一．有限元法的基础是单元刚度矩阵的推导，目前，刚度矩阵推导已有一个相对固定的模式，而烦琐、复杂的矩阵运算、微分、积分是刚度矩阵推导中的主要内容．通常，这种矩阵运算是由手工来完成的，工作量大，而且极易出错．利用ＭａｔＬａｂ丰富的符号运算功能，完成刚度矩阵推导过程中的符号运算，不但速度快，而且准确性高．Ｍａｆｌａｂ的命令和函数相当多，本文仅就涉及到ＭａｔＬａｂ符号运算作一些简要介绍，然后结合平面梁单元刚度矩阵的推导，用ＭａｔＬａｂ完成该单元刚度矩阵的符号运算，结果完全一致．１ＭａｔＬａｂ符号运算的一般规定１．１符号标变量和表达式的创建收稿日期：２００６－０２－０７；修订日期：２００６－０２－２３基金项目：重庆市教委项目（ＫＪ０５０４０７）；交通部西部交通建设科技项目（２００３３１８８１４２０２４）作者简介：周水兴（１９６７－），男，浙江嘉兴人，教授，主要从事桥梁设计理论研究．ｅ－ｍａｌｌ：ｓｈｕｉｘｉｎｇｚｈｏｕ＠ｃｑｕｃ．ｅｄｕ．ｃｎ３０重庆交通学院学报第２６卷ＭａｔＬａｂ中符号标量、变量和表达式是通过命令ｓｙｍ、ｓｙｍｓ来创建的．如：＞＞ｘ＝ｓｙｍ（’ｘ’）；＞＞ａｂｃ＝ｓｙｍ（’ａｌｐｈａ’）＞＞ｆｕｎ＝ｓｙｍ（’ａ宰ｘ＇２＋ｂ＇ｌｃＸ＋ｃ’）这里，＞＞是ＭａｔＬａｂ命令提示符，上述命令中前两个创建符号标量ｘ，ａｌｐｈａ，最后一个创建符号表达式ａ木ｘ２＋ｂ半ｘ＋ｃ．对符号表达式，还需指定表达式中的各个标量，如ａ，ｂ，Ｃ，Ｘ．ＭａｔＬａｂ中，除了上面介绍的方法外，还提供了一种简洁的方法，ｓｙｍｓ命令，如：＞＞ｓｙｍｓａｂＣＸ１．２符号矩阵在ＭａｔＬａｂ中，符号矩阵的元素为符号表达式，可用函数ｓｙｍ产生，如：＞＞Ｓ＝ｓｙｍ（’［ｘ，Ｙ，ｚ；ａ，ｂ，ｃ；ｕ，ｖ，Ｗ］’）该命令等价于：＞＞ｓｙｍｓｘＹｚａｂｃｕＶＷ＞＞Ｓ＝［ｘ，Ｙ，ｚ；ａ，ｂ，ｃ；ｕ，ｖ，ｗ］两者效果是一致的．对数值矩阵，也可以用ｓｙｍ把数值矩阵转换为符号形式．在转换过程中，函数ｓｙｍ尽可能地把数值矩阵中的非整元素指定为小的整数之比，以有理方式表示，如果元素是无理数，则在符号形式中ｓｙｍ将用符号浮点数表示元素．１．３符号矩阵的化简和格式化符号运算后，符号表达式有时非常复杂，为此，ＭａｔＬａｂ提供了能够化简复杂符号表达式的诸多命令，有ｐｒｅｔｔｙ（用于将符号表达式以常用方式显示）、ｃｏｌｌｅｃｔ（用于将表达式中同类项合并、合并后的多项式以变量的幂次方按大小依次排列）、ｅｘｐａｎｄ（用于展开符号表达式中的各项子式）、ｆａｃｔｏｒ（用于符号分解）、ｓｕｂｓ（较长子表达式的置换）、ｓｉｍｐｌｉｆｙ（用于符号表达式中进行等式的恒等替换）、ｓｉｍｐｌｅ等．１．４符号微积分Ｍａｔｌａｂ符号工具箱中提供了多种符号表达式的微分ｄｉｆｆ和积分ｉｎｔ函数形式，有：ｄｉｆｆ（ｆ），ｄｉｆｆ（ｆ，’ｔ’），ｄｉｌｌ（ｆ，ｎ），ｄｉｆｆ（ｆ，’ｔ’ｎ），ｉｎｔ（ｆ），ｉｎｔ（ｆ，’ｔ’），ｉｎｔ（ｆ，ａ，ｂ），ｉｎｔ（ｆ，’ｔ’，ａ，ｂ）和ｉｎｔ（ｆ，’ｍ’，’ｎ’）．有关微分ｄｉｆｆ和积分ｉｎｔ的详细内容可参考ＭａｔＬａｂ教程，这里不再赘述．２平面梁单元刚度矩阵为便于读者更好的理解ＭａｔＬａｂ在单元刚度矩阵推导中的应用，文中简要列出了平面梁单元刚度矩阵的推导过程．２．１平面梁单元位移模式‘２１对平面梁单元，单元的位移模式为Ⅱ＝ａｏ＋ａＩ石（１）秽＝ｂｏ＋ｂｌ省＋６２菇２＋ｂｓｘ３（２）单元节点位移列向量：…。

＝［ｕ；吼０ｉ叶吩岛］７利用单元ｉｊ的边界条件，解出上述系数，有ｕ＝［（１一ｘ／１）００ｘ／ｌ００］｛６｝。

＝札（髫）｛艿｝。

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

＝［＾Ｌ（菇）］｛６｝。

（４）２．２用结点位移表示应变和应力忽略剪切影响，梁单元受到拉压和弯曲变形后的线应变为：…＝ｈｄｕ／ｄ…ｘ：）＝一■…。

＝［Ｂ］｛艿｝。

（５）式中，［日］：∥引，１Ｌ一∥：（髫）Ｊ（６）［Ｂ］为应变矩阵．由胡克定律，单元中的应力表达式为：｛矿｝＝［Ｄ］｛占｝＝［Ｄ］［Ｂ］｛艿｝‘（７）２．３梁单元刚度矩阵根据虚位移原理，梁单元刚度矩阵为［ｋ］＝『ＪＪ［Ｂ］’［Ｄ］［Ｂ］ｄｙ（８）３符号矩阵的运算３．１应变矩阵［Ｂ］的计算应变矩阵中，有两个形函数表达式［札（菇）］、［见（石）］的微分，在Ｍａｔｌａｂ中定义应变矩阵和微分的程序段如下：Ｎｕ＝［１一ｘ／Ｌ００ｘ／Ｌ００］％定义形函数矩阵ＮｕＮｖ＝［０１—３，ｌｃｘ＇２／Ｌ＇２＋２，ｌｃｘ＇３／Ｌ３戈一２拳ｘ＇２／Ｌ＋ｘ＾３／Ｌ＇２０３’ｌｃｘ＇２／Ｌ＂２—２木ｘ＾３／Ｌ＾３一ｘ＇２／Ｌ＋ｘ＾３／Ｌ＇。

２］Ｂ＝［ｄｉｆｆ（Ｎｕ，’菇’，１）；一Ｙ术ｄｉｆｆ（Ｎｖ，’戈’，２）］％两个形函数分别进行变量戈一阶微分，然后形成第２期周水兴等：ＭａｔＬａｂ在有限元刚度矩阵推导中的应用３１应变矩阵［Ｂ］上述３个语句已经形成应变矩阵［Ｂ］，接下来就可以对刚度矩阵［Ｋ］进行积分．从语句Ｂ＝［ｄｉｆｆ（Ｎｕ，’并’，１）；一Ｙ水ｄｉｆｆ（Ⅳ匆，’茁’，２）］可以看出，ＭａｔＬａｂ由低维矩阵组成高维矩阵的方法非常简单，而且在形式上没有变化，容易被初学者理解，可见其卓越的功能．３．２刚度矩阵【Ｋ］积分和化简Ｄ＝［Ｅ０；ＯＥ］％定义弹性矩阵Ｋ＝Ｂ’枣Ｄ木Ｂ％矩阵相乘［Ｂ］＋［Ｄ］［Ｂ］Ｋ＝ｉｎｔ（Ｋ，’ｘ’，’Ｏ’，’Ｌ’）％对矩阵Ｋ中的每个元素沿长度方向积分Ｋ＝ｉｎｔ（Ｋ，’Ｙ’，’一ｈ／２’，’ｈ／２’）％矩阵Ｋ沿截面高度ｙ方向积分Ｋ＝ｉｎｔ（Ｋ，’Ｋ’，’一Ｗ／２’，’Ｗ／２’）％矩阵Ｋ沿截面宽度积分Ｋ＝ｓｕｂｓ（Ｋ，’ｈ‘３，ｌｃＷ’，’１２乖Ｉ’）％将矩阵中的ｈ３木Ｗ用１２１替换Ｋ＝ｓｕｂｓ（Ｋ，’ｈ水Ｗ’，’Ａ’）％将矩阵中的ｈ宰Ｗ用Ａ替换Ｋ＝ｓｉｍｐｌｉｆｙ（Ｋ）％积分后作恒等变换结果如下：［１／Ｌ枣Ｅ球Ａ，０，０，一１／Ｌ水Ｅ木Ａ，０，０Ｊ［０。

１２／Ｌ＾３术Ｅ书Ｉ，６／Ｌ‘２书ＥｊｌｃＩ，Ｏ，一１２／Ｌ＾３，ＩｃＥ术Ｉ，６／Ｌ＇２＇ｌｃＥ，ｌｃＩ］～［０，６／Ｌ＇２水Ｅ冰Ｉ，４／Ｌ木Ｅ木Ｉ，０，一６／Ｌ‘２水Ｅ，ＩｃＩ，２／Ｌ，ＩｃＥ木Ｉ］Ｋ＝．［一１／Ｌ木Ｅ，Ｉ‘Ａ，０，０，１／Ｌ水Ｅ，ＩｃＡ，０，０］［０，一１２／Ｌ＾３水ＥｊＩｃＩ，一６／Ｌ。

２木Ｅ木Ｉ，０，１２／Ｌ＾３术Ｅ宰Ｉ，一６／Ｌ‘２＝ｌｃＥ水Ｉ］［０，６／Ｌ＂２术Ｅ木Ｉ，２／Ｌ宰Ｅ木Ｉ，０，一６／Ｌ。

‘２水ＥＩ，４／Ｌｊｌ：Ｅ木Ｉ］运算结果与平面梁单元刚度矩阵完全一致心］．用该方法，已经成功地导出了空间梁单元几何非线同理，只需再增加几个语句，就可以快速地推导出等性大位移矩阵、初应力矩阵，为编写有限元程序奠定效荷载矩阵．了基础．特别指出的是，空间梁单元几何非线性中的上述程序可以用Ｍ文件保存，稍加修改可应用大位移矩阵，由于符号运算工作量巨大，在目前的书到其它单元刚度矩阵的推导中．籍中只介绍其理论，没有给出其显式刚度矩阵．对４结语此，我们将另文作专门介绍．本文以平面梁单元为例，详细介绍了ＭａｔＬａｂ在有限元刚度矩阵推导中的应用，不但为工程技术人员节省了大量的时间和精力，而且极大提高工作效率，对初学者深刻理解和运用有限元非常具有帮助．虽然文中以平面梁单元为例介绍其方法，但其中的原理和方法可十分方便地推广到其它单元刚度矩阵，甚至是非线性刚度矩阵的推导中，对此，我们利参考文献：［１］飞思科技产品研发中心．ＭａｔＬａｂ７基础与提高［Ｍ］．北京：电子工业出版社，２００５．［２］谢贻权，何福保．弹性和塑性力学中的有限单元法［Ｍ］．北京：机械工业出版社，１９８１．●●●－●－……－＿●－…－，－●－一●－－●●…一…－……－－●●－＿●－●●－●●一●●＿－●－●●－●－…●●●－●。