基于Matlab App Designer功能实现结构力学矩阵位移法

叶洪河,胡琼静

摘要：基于matlab的app designer功能，利用矩阵位移法思想，进行了平面杆系结构的矩阵位移法程序设计，实现矩阵位移法的GUI操作，并利用一个常见的例子进行演示。

关键词：matlab；矩阵位移法；有限元

引言

基于矩阵位移法的基本思想，可以采用不同的语言来进行程序实现，而matlab作为“矩阵实验室”，在矩阵位移法的编程实现上有着其得天独厚的优势。app designer作为matlab中可以实现图形用户界面（GUI）可视化与App行为编程功能的集成式环境，为广大用户提供了设计专属于自己的软件平台。综上所述，本文将介绍用matlab app designer模块实现矩阵位移法的程序，并举一个例子来加以说明。

1、程序描述

1.1步骤

以位移法为基础的矩阵位移法的基本步骤为：

基本信息的输入计算单元刚度矩阵整合整体刚度矩阵K计算等效节点力整合整体等效节点力求单元节点反力支座约束处理求节点位移求单元节点反力画弯矩图。

1.2基本信息的输入

在进行矩阵位移法求解时，我们需要利用有限元的思想，将结构离散化，划分单元，将单元、节点进行编号，然后向程序提供节点坐标，材料等信息。总结下来，我们需要向程序提供的信息如下：

节点：node（x坐标 y坐标 节点号）

单元：element（起始节点 终止节点 节点号）

材料：material（单元号 弹性模量E 截面惯性矩I 面积A ）

节间荷载：jjhz（单元号 节点1处值 节点2处值）

节点荷载：jdhz（节点号 自由度号 大小）

支座约束：zzys（节点号 自由度号 约束值）

耦合oh（节点1 节点2 耦合方向）

在本程序中，基本信息通过app界面中各子窗口来进行输入，这一功能的实现是通过matlab中的回调函数来实现的。

以节点输入为例，点击主界面中的“结点”按钮后，在弹出的窗口中输入节点信息即可。右边的空白部分便会画出相应位置处的节点和节点编号。相应的，在输入单元、材料、节点荷载等信息后，程序便会画出相应的模型。这时，基本信息的输入过程完成。

1.3计算单元刚度矩阵

单元刚度的建立通过一个单元刚度矩阵（function dg= dygdjz(dyh)）来完成，限于篇幅原因，代码会放在附件中。在读取单元的信息后，dygdjz子程序会根据节点坐标和材料信息来生成单元刚度矩阵。在这里，如果材料信息中弹性模量E为无穷大（即为桁架），可以在材料信息中输入“-1”，这时单元刚度矩阵为一个的矩阵。

此外，由于单元刚度矩阵的建立是在局部坐标中建立的，在整合之前，需要调用一个单元转化矩阵子程序（function T= dyzhjz(dyh)）将局部坐标转换为整体坐标。

1.4整合整体刚度矩阵

整体刚度矩阵的建立通过以下循环即可获得。不过需要注意的问题是仍需要判断是否为杠杆，这时循环需要做出一定的改变。

for i=1:1:2

        for j=1:1:2

            for p=1:1:3

                for q =1:1:3

                    m = (i-1)*3+p ;

                    n = (j-1)*3+q ;

                    M = (element(dyh,i)-1)*3+p ;

                    N = (element(dyh,j)-1)*3+q ;

                    zg(M,N) = zg(M,N) + dg(m,n);

                end

            end

        end

    end

1.5计算等效节点力

外力分为节点荷载和节间荷载两种，对于节间荷载，需要通过单元等效节点力子程序（function enf=dydxjdl(dyh,p1,p2)）将节间荷载转换成单元等效节点。对于节点荷载，直接在整合整体等效节点力时加上即可。

节间荷载的等效节点力可通过以下式子计算：

当节间荷载垂直于杆件时：

  enf(2) = (7*p1+3*p2)*cd/20;

  enf(3) = (3*p1+2*p2)*cd^2/60;

  enf(5) = (3*p1+7*p2)*cd/20;

  enf(6) = -(2*p1+3*p2)*cd^2/60;

1.6添加支座和耦合约束

约束的添加反映在整体刚度矩阵中就是划行划列，这一过程用代码表示出来比较复杂。我们可以利用乘大数法，我们以下面的例子来说明：

假设我们在自由度2处引入一个位移为2的约束：在刚度矩阵（2，2）位置乘一个大数，在荷载列阵处乘以2大数。如下式：

可得：.忽略和，即,

耦合约束的添加需要对荷载列阵和刚度矩阵进行以下的处理：

对于荷载列阵：

fp(m2) = fp(m2)+fp(m1);

    fp(m1) = 0;

对于刚度矩阵：

    zg(m2,:) = zg(m2,:)+zg(m1,:);

    zg(:,m2) = zg(:,m2)+zg(:,m1);

    zg(m1,:) = zeros(1, size(gNode,1)*3);

    zg(:,m1) = zeros(size(gNode,1)*3, 1);

    zg(m1,m1) = 1.0;

    zg(m1,m2) = -1;

    zg(m2,m1) = -1;

    zg(m2,m2) = zg(m2,m2)+1;

1.7 求解

由所得的整体刚度矩阵[K]、荷载列阵[]，通过式[K][X]=[],可求得节点处的位移。由所得位移列阵利用位移法可求得节点反力，从而可画出相应的内力图。

2、实例

一刚架受一节间荷载与节点荷载，画出其内力图，将该钢架划分成三个单元，将其节点进行编号。

通过各子窗口输入基本信息，完成建模。

这里需要输入的信息有;

节点：Node1:（0 0 1）Node2:(10 0 2) Node3:(0 10 3) Node4:(10 10 4)

单元：Element1:（1 3 1） Element1:（3 4 2） Element1:（4 2 3）

材料：material（1 10）

节间荷载：jjhz（2 10 20）

节点荷载：jdhz（2 1 10）

支座约束：zzys（1 1 0）（1 2 0）（1 3 0）（2 1 0）（2 2 0）（2 3 0）

用结构力学求解器验证可得结果正确。

3、总结

Matlab作为一款功能强大的编译器，能够在学术研究方面发挥很大的作用。此次利用其app designer实现了结构力学矩阵位移法，能够解决大多数问题。但此程序仍存在改进之处，比如：没有文件读取功能，数据输入比较麻烦等。通过日后的学习，希望可以将该程序完善。