新保存结构力学课设.docx

目录
结构力学课程设计任务书 2
第一章．绪论 7
（1）matlab简介 7
（2）matlab优势特点 7
（3）矩阵位移法简介 10
（4）矩阵位移法的特点 10
（5）矩阵位移法的基本原理 10
（6）矩阵位移法计算杆端力的步骤 11
第二章．Matlab解法 12
步骤一： 13
步骤二： 13
步骤三： 13
步骤四： 14
步骤五： 15
步骤六： 16
结论 19
附录程序 19
程序 19
子程序： 23
第三章．组内成员其他解法展示与对比 27
（1）Ansys有限元解法： 27
（2）矩阵位移法手算解法： 33
第四章．用 Matlab 完成所选矩阵方程题目计算 38
第五章．课程设计的心得体会。 40
参考文献 41

结构力学课程设计任务书
一、 设计题目
结构力学课程设计
二、 设计内容
1） 查阅资料并学习相关软件（Fortran或Matlab），对矩阵计算方程计算有所掌握，并完成所选矩阵方程题目计算；
2） 利用Fortran或Matlab软件，完成对所选结构力学题目的节点力及节点内力计算工作；
3） 查阅资料并学习相关软件（Ansys或Abaqus，Patran），对相关有限元软件的基本分析有所掌握，并完成所选结构力学题目的有限元分析计算；
4） 利用所学过的结构力学知识，完成对所选结构力学题目的手工计算工作，并完成相关内力图的绘制；
5） 完成对2）、3），4）三种计算结果进行对比分析。
三、 设计要求
1） 给出Fortran或Matlab矩阵计算程序的代码，并给出详细注释，并把计算结果抓图给出；
2） 给出Fortran或Matlab结构力学计算程序的代码，并给出详细注释，并把计算结果抓图给出；
3） 给出有限元软件计算结果（Mises应力、应变、位移、节点力）云图结果；
4） 完成题目的手工计算，并输入电脑，相关力学图的绘制可用Office的Visio组件或Autocad绘制；
5） 给出三种计算结果的对比（Excel曲线图或Origin）曲线图分析。
四、 注意事项
1） 时间：2016学年第2学期校历第17~18周（2016年6月20日至7月1日）；
2） 地点：所有任务均可自行选择完成地点；
3） 答疑地点为08A109；
4） 指导教师
王毅、代君、王晓璐、赵辉、莫振伟
第17周
周 一
周 二
周 三
周 四
周 五
教 师
王 毅
代 君
王晓璐
赵 辉
莫振伟
第18周
周 一
周 二
周 三
周 四
周 五
教 师
王 毅
代 君
王晓璐
赵 辉
莫振伟
5） 题目选择
每位学生需要在矩阵方程计算题目及结构力学计算题目中各选一道题目来完成；
6） 学生需要提交的课程设计报告应包括：前附《课程设计任务书》、设计要求所涉及的五项内容，参考文献；
7） 成绩评定及答辩
由于工作量较大，最终提交任务完成的时间节点为2016学年第2学期校历第19周周四（即2016年7月7日）。答辩时间暂定为2016年7月8日，成绩分为优秀、良好、中等、及格、不及格五等，具体时间另行通知。
五、 课程设计题目
一）矩阵方程
1. 利用全选主元的高斯约当（Gauss-Joadan）消去法求解如下方程组，并给出详细的程序注解和说明：

2. 利用追赶法求解如下方程组，并给出详细的程序注解和说明。

3. 利用全选主元的高斯约当（Gauss-Joadan）消去法如下求解大型稀疏矩阵的大型方程组，并给出详细注解及说明。

二） 结构力学
1. 试求解图示平面桁架各杆之轴力图，已知各材料性能及截面面积相同， 。（注：在有限元分析中，桁架杆的模拟只能选择Ansys的Link单元）。

2. 试求解图示平面刚架内力图（轴力图、剪力图和弯矩图），已知各材料性能及截面面积相同，，泊松比。

3. 试求解图示平面刚架内力图（轴力图、剪力图和弯矩图），已知各杆。

六、 工作计划

所在小组成员及分工情况
姓名
张世秋，
杨国，
朱春晖
顾凯强，
袁介周
刘永
孙子杨，
申正伟
李虎，
杨尚辉
任务
Matlab有限元分析
Ansys有限元分析
算法对比，误差分析
矩阵位移法手算
Matlab列主元高斯消元法解矩阵方程
我负责的是Matlab有限元分析，我们组分工的情况，如上所示。
我一开始先查相关的资料，了解matlab工作原理，查阅相关算例，仿照书上的算例，编制相关的程序，尝试计算。
但是程序编制并不是一蹴而就，在编制的过程中，出了很多的相关问题，找不到问题的所在，最终我听从老师的意见，然后终于查到问题所在，顺利解决问题。

课程负责人签名： 指导教师签名：

年 月 日 年 月 日

第一章．绪论
（1）matlab简介
MATLAB[1] 是美国MathWorks公司出品的商业数学软件，用于算法开发、数据可视化、数据分析以及数值计算的高级技术计算语言和交互式环境，主要包括MATLAB和Simulink两大部分。
MATLAB是matrix&laboratory两个词的组合，意为矩阵工厂（矩阵实验室）。是由美国mathworks公司发布的主要面对科学计算、可视化以及交互式程序设计的高科技计算环境。它将数值分析、矩阵计算、科学数据可视化以及非线性动态系统的建模和仿真等诸多强大功能集成在一个易于使用的视窗环境中，为科学研究、工程设计以及必须进行有效数值计算的众多科学领域提供了一种全面的解决方案，并在很大程度上摆脱了传统非交互式程序设计语言（如C、Fortran）的编辑模式，代表了当今国际科学计算软件的先进水平。
MATLAB和Mathematica、Maple并称为三大数学软件。它在数学类科技应用软件中在数值计算方面首屈一指。MATLAB可以进行矩阵运算、绘制函数和数据、实现算法、创建用户界面、连接其他编程语言的程序等，主要应用于工程计算、控制设计、信号处理与通讯、图像处理、信号检测、金融建模设计与分析等领域。
MATLAB的基本数据单位是矩阵，它的指令表达式与数学、工程中常用的形式十分相似，故用MATLAB来解算问题要比用C，FORTRAN等语言完成相同的事情简捷得多，并且MATLAB也吸收了像Maple等软件的优点，使MATLAB成为一个强大的数学软件。在新的版本中也加入了对C，FORTRAN，C++，JAVA的支持。
（2）matlab优势特点
1) 高效的数值计算及符号计算功能，能使用户从繁杂的数学运算分析中解脱出来；
2) 具有完备的图形处理功能，实现计算结果和编程的可视化；
3) 友好的用户界面及接近数学表达式的自然化语言，使学者易于学习和掌握；
4) 功能丰富的应用工具箱(如信号处理工具箱、通信工具箱等) ，为用户提供了大量方便实用的处理工具。
编程环境
MATLAB由一系列工具组成。这些工具方便用户使用MATLAB的函数和文件，其中许多工具采用的是图形用户界面。包括MATLAB桌面和命令窗口、历史命令窗口、编辑器和调试器、路径搜索和用于用户浏览帮助、工作空间、文件的浏览器。随着MATLAB的商业化以及软件本身的不断升级，MATLAB的用户界面也越来越精致，更加接近Windows的标准界面，人机交互性更强，操作更简单。而且新版本的MATLAB提供了完整的联机查询、帮助系统，极大的方便了用户的使用。简单的编程环境提供了比较完备的调试系统，程序不必经过编译就可以直接运行，而且能够及时地报告出现的错误及进行出错原因分析。
简单易用
Matlab是一个高级的矩阵/阵列语言，它包含控制语句、函数、数据结构、输入和输出和面向对象编程特点。用户可以在命令窗口中将输入语句与执行命令同步，也可以先编写好一个较大的复杂的应用程序（M文件）后再一起运行。新版本的MATLAB语言是基于最为流行的C++语言基础上的，因此语法特征与C++语言极为相似，而且更加简单，更加符合科技人员对数学表达式的书写格式。使之更利于非计算机专业的科技人员使用。而且这种语言可移植性好、可拓展性极强，这也是MATLAB能够深入到科学研究及工程计算各个领域的重要原因。
强处理能力
MATLAB是一个包含大量计算算法的集合。其拥有600多个工程中要用到的数学运算函数，可以方便的实现用户所需的各种计算功能。函数中所使用的算法都是科研和工程计算中的最新研究成果，而且经过了各种优化和容错处理。在通常情况下，可以用它来代替底层编程语言，如C和C++ 。在计算要求相同的情况下，使用MATLAB的编程工作量会大大减少。MATLAB的这些函数集包括从最简单最基本的函数到诸如矩阵，特征向量、快速傅立叶变换的复杂函数。函数所能解决的问题其大致包括矩阵运算和线性方程组的求解、微分方程及偏微分方程的组的求解、符号运算、傅立叶变换和数据的统计分析、工程中的优化问题、稀疏矩阵运算、复数的各种运算、三角函数和其他初等数学运算、多维数组操作以及建模动态仿真等。
图形处理
MATLAB自产生之日起就具有方便的数据可视化功能，以将向量和矩阵用图形表现出来，并且可以对图形进行标注和打印。高层次的作图包括二维和三维的可视化、图象处理、动画和表达式作图。可用于科学计算和工程绘图。新版本的MATLAB对整个图形处理功能作了很大的改进和完善，使它不仅在一般数据可视化软件都具有的功能（例如二维曲线和三维曲面的绘制和处理等）方面更加完善，而且对于一些其他软件所没有的功能（例如图形的光照处理、色度处理以及四维数据的表现等），MATLAB同样表现了出色的处理能力。同时对一些特殊的可视化要求，例如图形对话等，MATLAB也有相应的功能函数，保证了用户不同层次的要求。另外新版本的MATLAB还着重在图形用户界面（GUI）的制作上作了很大的改善，对这方面有特殊要求的用户也可以得到满足。
MATLAB对许多专门的领域都开发了功能强大的模块集和工具箱。一般来说，它们都是由特定领域的专家开发的，用户可以直接使用工具箱学习、应用和评估不同的方法而不需要自己编写代码。领域，诸如数据采集、数据库接口、概率统计、样条拟合、优化算法、偏微分方程求解、神经网络、小波分析、信号处理、图像处理、系统辨识、控制系统设计、LMI控制、鲁棒控制、模型预测、模糊逻辑、金融分析、地图工具、非线性控制设计、实时快速原型及半物理仿真、嵌入式系统开发、定点仿真、DSP与通讯、电力系统仿真等，都在工具箱（Toolbox）家族中有了自己的一席之地。
程序接口
新版本的MATLAB可以利用MATLAB编译器和C/C++数学库和图形库，将自己的MATLAB程序自动转换为独立于MATLAB运行的C和C++代码。允许用户编写可以和MATLAB进行交互的C或C++语言程序。另外，MATLAB网页服务程序还容许在Web应用中使用自己的MATLAB数学和图形程序。MATLAB的一个重要特色就是具有一套程序扩展系统和一组称之为工具箱的特殊应用子程序。工具箱是MATLAB函数的子程序库，每一个工具箱都是为某一类学科专业和应用而定制的，主要包括信号处理、控制系统、神经网络、模糊逻辑、小波分析和系统仿真等方面的应用。
应用软件开发
在开发环境中，使用户更方便地控制多个文件和图形窗口；在编程方面支持了函数嵌套，有条件中断等；在图形化方面，有了更强大的图形标注和处理功能，包括对性对起连接注释等；在输入输出方面，可以直接向Excel和HDF5进行连接。
（3）矩阵位移法简介
在结构力学的计算中，通过采用对结点位移作为基本未知量，进而通过矩阵的形式对各基本参数进行组织，编排，求出未知量的方法，称为矩阵位移法。
（4）矩阵位移法的特点
矩阵数学表达力强
矩阵数学表达力强，运算简洁方便并且适于计算机组织运算，是用计算机进行结构数值分析的最强有力的数学工具。
矩阵位移法与结构力学的力法和位移法相对应，也就是结构的矩阵分析方法。
矩阵位移法方便编制程序
矩阵位移法便于编制程序，因而在工程界得到广泛应用。
矩阵位移法并不因采用矩阵数学的描述手段，而改变位移法的基本原理。它与位移法的区别仅仅在于表达形式不同。
（5）矩阵位移法的基本原理
按位移法的基本原理运用矩阵计算内力和位移的方法。是结构矩阵分析方法中的一种，其基本未知数是结点位移，由于矩阵位移法较矩阵力法更适宜编制通用的计算程序，因而得到了更为广泛的应用。
结构矩阵分析方法首先把结构离散成有限数目的单元，然后再合成为原结构，因而也属于有限元法。矩阵位移法常用的单元形式为一直杆。对于曲杆，如拱结构，虽然也可取曲杆作为单元,但单元分析较烦,为简化起见，可将它化成折线来处理，每一直线段作为一单元。当单元承受非结点荷载时，可
用等效结点荷载代替。其方法是将单元间的分界结点作为固端求出固端反力，然后反其向作用在结点上。
（6）矩阵位移法计算杆端力的步骤
矩阵位移法计算杆端力的步骤为：
①划分单元，求出等效结点荷载；
②求单元刚度矩k(i),并转换为整体坐标的单元刚度矩阵;
③由刚度法求出整体刚度矩阵K；
④求出Krr；
⑤求出结点位ur，再求出杆端s，实际杆端力应再叠sf， 即确定。

第二章．Matlab解法

为了便于计算和书写对图形做如下处理。其中p=20KN,L=2m。
因为结构是对称的，载荷也是对称的，为了简化计算，可以利用对称性只取结构的一半进行计算，如下图所示，因为节点2,4，没有水平位移，所以各加了一根水平铰链，又因为4,2，杆在对称轴上，所以它的横截面积取A/2，
各节点位移的编号如图所示，各基本数据的计算如下表所示。

步骤一：
离散化域，域分为8个单元5个节点，matlab中采用的单位是kN和m，下表给出了本题的单元连通性。

步骤二：
写出单元刚度矩阵，通过matlab函数PlaneTrussElementStiffness，得到8个单元的刚度矩阵k1,k2,k3,……,k8。每个矩阵都是4*4矩阵。
步骤三：
集成整体刚度矩阵。由于该结构一共有5个节点，所以整体刚度矩阵是10*10矩阵。因此为了得到整体刚度矩阵K,我们首先要生成一个10*10的零矩阵。由于该结构有8个单元，8次调用matlab的PlaneTrussAssemble函数就可以得到整体刚度矩阵K。每次对该函数调用都集成一个单元。
步骤四：
引入边界条件。用上一步得到的整体刚度矩阵，可以得到该结构的矩阵方程。

F1x
F1y
F2x
F2y
F3x
F3y
F4x
F4y
F5x
F5y

U1x
U1y
U2x
U2y
U3x
U3y
U4x
U4y
U5x
U5y

[K]* =

0
-20
F2x
0
0
0
F4x
F4y
F5x
F5y

U1x
U1y
0
U2y
U3x
U3y
0
0
0
0

引入边界条件U2x=U4y= U4x= U5x=U5y=0，F1x=F2y=F3x=F3y=0，F1y=-20KN

[K]* =

步骤五：
用分解（手动），和高斯消元法（应用matlab）求解方程组。首先对方程组进行分解，提取整体刚度矩阵K的第1，2，4，5，6行及列的交集组成子矩阵[K1]。如下图所示：

从而得到：

0
-20
0
0
0

U1x
U1y
U2y
U3x
U3y

[K]* =

上述方程组的解用如matlab命令可以得到。
解得u =

-0.0000
-0.0024
-0.0008
0.0002
-0.0020

现在可以得知各节点位移。

步骤六：
在这一步中，首先我们建立结构节点位移矢量U
U=[0;-0.0024;0; -0.0008; 0.0002; -0.0020;0;0;0;0];
然后我们建立单元节点位移矢量u1,u2,u3,u4,u5,u6,u7,u8。
u1=[U(1);U(2);U(3);U(4)];
u2=[U(3);U(4);U(7);U(8)];
u3=[U(7);U(8);U(5);U(6)];
u4=[U(5);U(6);U(9);U(10)];
u5=[U(9);U(10);U(1);U(2)];
u6=[U(1);U(2);U(5);U(6)];
u7=[U(3);U(4);U(5);U(6)];
u8=[U(1);U(2);U(7);U(8)];

然后调用matlab的PlaneTrussElementStress函数计算出单元应力sigma1,sigma2, sigma3,sigma4, sigma5,sigma6, sigma7,sigma8。
sigma1 =

0

sigma2 =

84000

sigma3 =

21000

sigma4 =

-21000

sigma5 =

-1.2600e+05

sigma6 =

4.2000e+04

sigma7 =

-5.2500e+04

sigma8 =

1.2600e+05
结论
于是得到所有杆件的轴向应力：

附录程序
所有程序及其子程序如下：
程序
E=210*10^6;
A=90.7*10^-6;
A2=0.5*A;

L1=2;theta1=0;
L2=2;theta2=90;
L3=2;theta3=0;
L4=2; theta4=0;
L5=PlaneTrussElementLength(0,0,2,2);theta5=45;
L6=2;theta6=90;
L7=L5;theta7=45;
L8=L5;theta8=135;

k1=PlaneTrussElementStiffness(E,A,L1,theta1);

k2=PlaneTrussElementStiffness(E,A2,L2,theta2);

k3=PlaneTrussElementStiffness(E,A,L3,theta3);

k4=PlaneTrussElementStiffness(E,A,L4,theta4);

k5=PlaneTrussElementStiffness(E,A,L5,theta5);

k6=PlaneTrussElementStiffness(E,A,L6,theta6);

k7=PlaneTrussElementStiffness(E,A,L7,theta7);

k8=PlaneTrussElementStiffness(E,A,L8,theta8);

K=zeros(10,10);
K=PlaneTrussAssemble(K,k1,1,2);
K=PlaneTrussAssemble(K,k2,4,2);
K=PlaneTrussAssemble(K,k3,3,4);
K=PlaneTrussAssemble(K,k4,5,3);
K=PlaneTrussAssemble(K,k5,5,1);
K=PlaneTrussAssemble(K,k6,3,1);
K=PlaneTrussAssemble(K,k7,3,2);
K=PlaneTrussAssemble(K,k8,4,1);

k=K([1,2,4,5,6],[1,2,4,5,6])
f=[0;-20;0;0;0];
u=k\f

u=[0; -0.0024; -0.0008; 0.0002; -0.0020];

U=[0;-0.0024;0; -0.0008; 0.0002; -0.0020;0;0;0;0 ];

u1=[U(1);U(2);U(3);U(4)];
u2=[U(3);U(4);U(7);U(8)];
u3=[U(7);U(8);U(5);U(6)];
u4=[U(5);U(6);U(9);U(10)];
u5=[U(9);U(10);U(1);U(2)];
u6=[U(1);U(2);U(5);U(6)];
u7=[U(3);U(4);U(5);U(6)];
u8=[U(1);U(2);U(7);U(8)];
sigma1=PlaneTrussElementStress(E,L1,theta1,u1)
sigma2=PlaneTrussElementStress(E,L2,theta2,u2)
sigma3=PlaneTrussElementStress(E,L3,theta3,u3)
sigma4=PlaneTrussElementStress(E,L4,theta4,u4)
sigma5=PlaneTrussElementStress(E,L5,theta5,u5)
sigma6=PlaneTrussElementStress(E,L6,theta6,u6)
sigma7=PlaneTrussElementStress(E,L7,theta7,u7)
sigma8=PlaneTrussElementStress(E,L8,theta8,u8)

子程序：
function y=PlaneTrussElementStress(E,L,theta,u)
%该函数根据弹性模量E,长度L,角度theta(单位是度)以及单位节点位移矢量u计算单元应力。它返回单位应力大小，返回值是一个标量而不是矢量。
x=theta*pi/180;
C=cos(x);
S=sin(x);
y=E/L*[-C -S C S]*u;

function y=PlaneTrussInclinedSupport(T,i,alpha)%该函数根据倾斜支柱的节点号i,以及倾斜角alpha计算得到倾斜支柱的变换矩阵。它返回2n*2n的变换矩阵。
x=alpha*pi/180;
T(2i-1,2i-1)=cos(x);
T(2i-1,2i)=sin(x);
T(2i-1,2i-1)=-sin(x);
T(2i,2i)=cos(x);
y=T;

function y=PlaneTrussElementStiffness(E,A,L,theta)
%该函数根据每个平面桁架元的弹性模量E，横截面积A,长度L,
%以及角度theta(单位是度)计算得到单元刚度矩阵，它返回4*4的单元刚度矩阵。
x=theta*pi/180;
C=cos(x);
S=sin(x);
y=E*A/L*[C*C C*S -C*C -C*S; C*S S*S -C*S -S*S
-C*C -C*S C*C C*S; -C*S -S*S C*S S*S];

function y=PlaneTrussElementLength(x1,y1,x2,y2) %该函数根据给出的第一个坐标（x1,y1）和第二个坐标（x2,y2）计算返回单元长度。
y=sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1));

function y=PlaneTrussElementForce(E,A,L,theta,u)%该函数根据弹性模量E,长度L,角度theta(单位是度)以及单元节点位移矢量u计算单元节点力。它以标量形式返回单元节点力。
x=theat*pi/180;
C=cos(x);
S=sin(x);
y=E*A/L*[-C -S C S]*u;

function y=PlaneTrussAssemble(K,k,i,j)%该函数将连接i和j的平面桁架元的刚度矩阵k集成到整体刚度矩阵K。
%每集成一个单元，该函数都会返回2n*2n的整体刚度矩阵K.
K(2*i-1,2*i-1)=K(2*i-1,2*i-1)+k(1,1);
K(2*i-1,2*i)=K(2*i-1,2*i)+k(1,2);
K(2*i-1,2*j-1)=K(2*i-1,2*j-1)+k(1,3);
K(2*i-1,2*j)=K(2*i-1,2*j)+k(1,4);

K(2*i,2*i-1)=K(2*i,2*i-1)+k(2,1);
K(2*i,2*i)=K(2*i,2*i)+k(2,2);
K(2*i,2*j-1)=K(2*i,2*j-1)+k(2,3);
K(2*i,2*j)=K(2*i,2*j)+k(2,4);

K(2*j-1,2*i-1)=K(2*j-1,2*i-1)+k(3,1);
K(2*j-1,2*i)=K(2*j-1,2*i)+k(3,2);
K(2*j-1,2*j-1)=K(2*j-1,2*j-1)+k(3,3);
K(2*j-1,2*j)=K(2*j-1,2*j)+k(3,4);

K(2*j,2*i-1)=K(2*j,2*i-1)+k(4,1);
K(2*j,2*i)=K(2*j,2*i)+k(4,2);
K(2*j,2*j-1)=K(2*j,2*j-1)+k(4,3);
K(2*j,2*j)=K(2*j,2*j)+k(4,4);
y=K;

第三章．组内成员其他解法展示与对比
（1）Ansys有限元解法：
第一步：进入 ANSYS(设定工作目录和工作文件)
进入ANSYS菜单路径 “程序 >ANSYS >ANSYS12.0”
设置工作文件名
菜单路径 “file > Change Jobname”，弹出“Change Jobname”对话框，输入“hengjia”，单击OK确定并关闭对话框。
设置计算类型
菜单路径 “ANSYS Main Menu: Preferences”，在弹出的对话框中选择“Structural” ，单击OK确定并关闭对话框。
选择单元类型
菜单路径 “ANSYS Main Menu: Preprocessor >Element Type>Add/Edit/Delete…”，在弹出“Library of Element Types”对话框中按照参数选择，单击OK确定并关闭对话框。
第二步：定义实常数
菜单路径 “ANSYS Main Menu: Preprocessor >Real Constants… >Add/Edit/Delete ”，在弹出的对话框中单击“Add > OK”， 弹出“Real Constant …”对话框，参数设置“AREA 0.00325”，单击OK确定并关闭对话框。
第三步：定义材料参数
菜单路径 “ANSYS Main Menu: Preprocessor > Material Props > Material Models”，在弹出的菜单中打开“Structural > Linear > Elastic > Isotropic”，弹出“Linear Isotropic Material…”对话框，并设置参数。
第四步：定义梁的截面
菜单路径 “ANSYS Main Menu: Preprocessor > Sections > Beam > Common Sections”，弹出“Beam Tool”对话框，并设置，单击OK确定关闭对话框。
生成几何模型，结果如图所示。

第五步：划分网格，生成有限元模型
菜单路径 “ANSYS Main Menu: Preprocessor > Meshing > Mesh Tool”， “Mesh Tool”对话框。单击其中Lines旁边的Set按钮，打开 “Element Size on Picked Lines”选择对话框，然后输入参数：
第六步：在mesh框中单击Pick All按钮，生成有限元模型如下
第七步:模型施加约束和外载
在1点处施加自由度为零的约束
菜单路径 “ANSYS Main Menu: Solution > Define Loads > Apply >Structural > Displayment > On Keypoints”，弹出的对话框，选择关键点1，单击Apply弹出“Apply U，ROT on KPs”对话框，选择“ALL DOF”， 单击OK确定关闭对话框。
第八步：在7点处施加Y方向的约束方法同上。设置重力加速度g
第九步：分析计算
菜单路径 “ANSYS Main Menu: Solution > Solve > Current LS”，在弹出的对话框中单击OK，并关闭文字窗口 。
第十步：显示变形图
菜单路径 “ANSYS Main Menu: General Postproc > Plot Results >Deformed Shape … ”，弹出“Plot Deformed Shape”对话框，设置“Def +Undeformed” 单击OK确定关闭对话框。
第十一步：显示位移云图
菜单路径 “ANSYS Main Menu: General Postproc > Plot Results > Nodal Solu”，弹出 “Contour Nodal Solution Data”对话框，并设置，单击OK查看位移云图。结果如下图所示。

第十二步：显示轴力云图

节点力图

最终输出轴力的大小如下所示：
PRINT ELEMENT TABLE ITEMS PER ELEMENT

***** POST1 ELEMENT TABLE LISTING *****

STAT CURRENT
ELEM SMIS1
1 1826.5
2 -11559.
3 -3653.1
4 5166.2
5 -1826.5
6 -11559.
7 -7306.1
8 -1826.5
9 -11559.
10 1826.5
11 5166.2
12 -3653.1
13 -11559.
14 0.0000
15 0.0000

MINIMUM VALUES
ELEM 2
VALUE -11559.

MAXIMUM VALUES
ELEM 4
VALUE 5166.2
将其转化为如下表所示：

（2）矩阵位移法手算解法：
1.试求解图示平面桁架各杆之轴力图，已知各材料性能及截面面积相同，。
为了便于计算和书写对图形做如下处理。其中p=20KN,L=2m。
因为结构是对称的，载荷也是对称的，为了简化计算，可以利用对称性只取结构的一半进行计算，如下图所示，因为节点2,4，没有水平位移，所以各加了一根水平铰链，又因为4,2，杆在对称轴上，所以它的横截面积取A/2，
各节点位移的编号如图所示，各基本数据的计算如下表所示。

应用式：

Ui vi uj vj

Xij
Yij
Xji
Yji

K=EF/L(cosθ)^2cosθsinθcosθsinθ(sinθ)^2-(cosθ)^2-cosθsinθ-cosθsinθ-(sinθ)^2-(cosθ)^2-cosθsinθ-cosθsinθ-(sinθ)^2(cosθ)^2cosθsinθcosθsinθ(sinθ)^2
计算各单元劲度矩阵如下：

K12=K34=K53=EF/L1000-1000-10001000
K31=EF/L0001/2000-1/2000-1/20001/2
K42=EF/L0001/2000-1/2000-1/20001/2
K32=K51=EF/L2/42/42/42/4-2/4-2/4-2/4-2/4-2/4-2/4-2/4-2/42/42/42/42/4
K41=EF/L2/4-2/4-2/42/4-2/42/42/4-2/4-2/42/42/4-2/42/4-2/4-2/42/4
用对号入座的方法组成结构劲度矩阵。

K=EF/L1.7071001.707100000-100000-10.8536-0.3536-0.3536-0.35362.35360.3536-0.35360.35361.3536
列出结构劲度方程如下：

EF/L1.7071001.707100000-100000-10.8536-0.3536-0.3536-0.35362.35360.3536-0.35360.35361.3536*u1v1v2u3v3=0-p000
由上列方程解得未知结点位移：
u1v1v2u3v3=PL/EF0-1.1561-0.36540.0914-0.9735

根据式

Nij=EF/L[COSθ(uj-ui)+sinθ(vj-vi)]

列表计算出如下表所示的桁架各杆的内力。

由于上图所示的桁架是由原来的桁架截取的一半，所以杆42的真实内力应将表中的数值乘以2.这样如下表所示桁架的各杆的内力。

（3）三种解法的对比与分析：

三种解法的数据对比显示于下列表格和折线图中：

从图表中很容易得出结论，整体来看，三种计算方法得出的结果是一致的，但使用 ansys
有限元方法进行分析，得出的计算结果更加接近于手算得出的结果。而且在实际的练习中发 现使用 ansys 有限元方法比较生动形象，更加易于操作和求解，省去了很多编程计算所不必 要的麻烦。
而使用 matlab 求解此次出现的误差较大主要是在编程的过程中，为了减少不必要的麻 烦和出错的可能性，选择的计算方法略为粗糙，所以得到的计算结果不甚精确。

第四章．用 Matlab 完成所选矩阵方程题目计算

3.利用全选主元的高斯约当（Gauss-Joadan）消去法如下求解大型稀疏矩阵的大型方程组， 并给出详细注解及说明。

Matlab 计算程序如下所示：

clear
A=input('输入系数矩阵A：');
b=input('输入b向量（按