下载文档原格式
第三章连续梁程序的编制与使用入结构力学领域中而产生的一种方法，而Matlab语言正是进行矩阵运算的强大工具，因此，用Matlab语言编写结构力学程序有更大的优越性。
本章将详细介绍如何利用Matlab语言编制连续梁结构的计算程序。
矩阵位移法的解题思路是将结构离散为单元（杆件），建立单元杆端力与杆端位移之间的关系－单元刚度方程；再将各单元集成为原结构，在满足变形连续条件和平衡条件时，建立整体刚度方程；在边界条件处理完毕后，由整体刚度方程解出节点位移，进而求出结构内力。
用矩阵位移法计算连续梁的步骤如下：1）整理原始数据，如材料性质、荷载条件、约束条件等，并进行编码：单元编码、结点编码、结点位移编码、选取坐标系。
2）建立局部坐标系下的单元刚度矩阵。
3）建立整体坐标系下的单元刚度矩阵。
4）集成总刚。
5）建立整体结构的等效节点荷载和总荷载矩阵6）边界条件处理。
7）解方程，求出节点位移。
8）求出各单元的杆端内力。
实际上，上述步骤也是编制Matlab程序的基本步骤，在求出计算结果后，还可以利用Matlab的绘图功能绘制结构图、内力图、变形图等等。
图3-1程序流程图3.1 程序说明%******************************************************************* % 矩阵位移法解连续梁主程序%******************************************************************* ●功能：运用矩阵位移法解连续梁的基本原理编制的计算主程序。
●基本思想：结点（结点位移）编码默认为从左至右，从1开始顺序进行；杆端弯矩的方向默认为逆时针。
●荷载类型：可计算结点荷载，每单元作用的跨中集中力和均布荷载。
●说明：主程序的作用是通过赋值语句、读取和写入文件、函数调用等完成算法的全过程，即实现程序流程图的程序表达。
%-----------------------------------------------------------------------------------------------------1 程序准备format short e %设定输出类型clear all %清除所有已定义变量clc %清屏●说明：format short e －设定计算过程中显示在屏幕上的数字类型为短格式、科学计数法；clear all －清除所有已定义变量，目的是在本程序的运行过程中，不会发生变量名相同等可能使计算出错的情况；clc －清屏，使屏幕在本程序运行开始时%-----------------------------------------------------------------------------------------------------2 打开文件FP1=fopen('input.txt','rt'); %打开输入数据文件存放初始数据FP2=fopen('output.txt','wt'); %打开输出数据文件存放计算结果●说明：FP1=fopen('input.txt','rt'); －打开已存在的输入数据文件input.txt，且设置其为只读格式，使程序在执行过程中不能改变输入文件中的数值，并用文件句柄FP1来FP2=fopen('output.txt','wt'); －打开输出数据文件，该文件不存在时，通过此命令创建新文件，该文件存在时则将原有内容全部删除。
1、 问题的描述算例1 用MATLAB 写无限长梁在受集中力作用下弯矩，剪力，挠度，转角的表达式，要求绘出相应的分布图明之处。
图3-8在硬性粘土的天然地基上，有一无限长梁，基床系数k =4x104KN/m 4，无限长梁的宽度b 为0.5m ，s k =k b=2 x104KN/m 4，E C =3x104MPa,梁的高度为0.8m,I =2.1X10-2跑m 4,λ= 0.298/m ，则L 为特征长度取3.36m ，考虑到无限长梁的特征，取梁长为32m ，集中力P 0=30KN ，作用于梁的跨中位置，满足16m ﹥10.55m 。
2、 无限长梁的理论分析44d 0dxs w k w += （3-21）四阶常系数线性常微分方程（3-21）的通解为：1234(C cos sin )(C cos sin )x x w e x C x e x C x λλλλλλ-=+++（3-22）式（3-22）中含有四个积分常数C 1、C 2、C 3、C 4 ，可由荷载情况和边界条件确定。
由材料力学可知，对梁的挠度w 求一次导师可求出梁的转角ϕ,求两次导数可求出梁的弯矩M,求三次导数可得梁的剪力Q ,即2233dw dxc cd w E I M dx d wE I Q dx ϕ⎫=⎪⎪⎪-=⎬⎪⎪-=⎪⎭（3-23）1）无限长梁受集中荷载P 作用下的解答无限长梁见图3-8),对于这样的无限长梁可以得到三个边界条件，并由此求出积分常数C 。
（1）当x →∞时，0w =。
把这个边界条件代入式（3-22），即1234(C cos sin )(C cos sin )0x x x w e x C x e x C x λλλλλλ-→∞=+++=在上式中，当x →∞时e 0,e-0,x x λλ≠=由此得到待定常数C1=C2=0，式（3-22）成为: 34(C cos sin )x w e x C x λλλ-=+ （3-24） （2）当x=0,时，dw==0dxϕ。
/prep7L1=30 !设置变量L2=30h=-25K, 1, 0, 0, 0,K, 2, L1, 0, 0,K, 3, L1, L2, 0,K, 4, 0, L2, 0,KWPA VE, 1 !将工作平面原点定义在1号点RECTNG, 0, L1, 0, L2,wpro, , -90, !将工作平面绕X轴Z到Y方向90度RECTNG, 0, L1, 0, -h,KWPA VE, 4 !将工作平面原点定义在4号点RECTNG, 0, L1, 0, -h,wpro, , ,90 !将工作平面绕y轴x到z方向90度RECTNG, 0, L2, 0, -h,KWPA VE, 3 !将工作平面原点定义在3号点RECTNG, 0, L2, 0, -h,AGLUE, all !粘结所有面ET, 1, SHELL43 !ET，ITYPE，Ename，KOPT1，~，KOPT6，INOPR（定义单元）!KOPT1～KOPT6为元素特性编码!shell43 4 节点塑性大应变单元ET, 2, COMBIN14 !COMBIN14弹簧-阻尼器Spring-Damper MPTEMP,,,,,,,, !删除系统中已存在的温度表MPTEMP, 1, 0 !定义一个温度表MPDATA, EX, 1, , 2.4E10 !指定与温度相应的材料性能数据弹性模量MPDATA, PRXY, 1, , 0.15 !主泊松比ESIZE, 1, 0 !指定单元边长AMESH, ALL !划分面生成面单元NSEL, S, LOC, Z, 0 !选择一组节点子集创建新集ESLN, S !选择已选节点上的单元NSEL, S, LOC, Z, -1 !选择z坐标值为-1的---ESLN, U !从已选集中删除此时剩下只支撑板CM, STRUT, ELEM !将选择集命名STRUT生成元件alls !all sel 全选CMSEL, U, STRUT !去除STRUT元件CM, W ALL, ELEM !将选择集命名wall生成元件NSEL, S, LOC, X, 0.1, L1-0.1 !选择一组节点子集创建新集NPLOT !显示节点NSEL, R, LOC, Y, 0 !从当前集选择一组节点子集ESLN, S !从已选集中选择NSEL, S, LOC, Y, 1 !从当前集选择一组节点子集ESLN, U !从已选集中删除ENSYM, , , , ALL !反转壳单元法线方向NSEL, S, LOC, Y, 0.1, L2-0.1 !选择一组节点子集创建新集NPLOT !显示节点NSEL, R, LOC, X, 0 !从当前集选择一组节点子集ESLN, S !从已选集中选择NSEL, S, LOC, X, 1 !从当前集选择一组节点子集ESLN, U !从已选集中删除ENSYM , , , , ALL !反转壳单元法线方向ALLSNUMCMP, ALL !所有实体进行重新编号!直接生成节点*DO, i, 1, L1-1 ! 从1到29进行循环CSYS, 0 !激活默认笛卡尔坐标系N, 100000+2*i-1, i, 0, -1 ! 节点编号后面为坐标N, 100000+2*i, i, 2.5, -1 ! 节点编号后面为坐标*enddo*DO, i, 1, L1-1 ! 从1到29进行循环!Modeling>Creat>Elements>Elem AttributesTYPE, 2 !设置单元类型属性指示器MAT , 1 !MP命令中的MA T即材料性能REAL, 0 !材料实常数ESYS, 0 !材料坐标系统属性指示器EN, 100000+i, 100000+2*i-1, 100000+2*i !根据给定的单元号和节点号生成单元*enddoNUMCMP , ALL !所有实体进行重新编号*DO, i, 1, L1-1CSYS, 0 !激活默认笛卡尔坐标系N, 100000+2*i-1, i, L2, -1N, 100000+2*i, i, L2-2.5, -1*enddo*DO, i, 1, L1-1TYPE, 2MAT , 1REAL, 0ESYS, 0EN, 100000+i, 100000+2*i-1, 100000+2*i*enddoNUMCMP , ALL !所有实体进行重新编号*DO, i, 1, L2-1CSYS, 0N, 100000+2*i-1, 0, i, -1N, 100000+2*i, 2.5, i, -1*enddo*DO, i, 1, L2-1TYPE, 2MAT , 1REAL, 0ESYS, 0EN, 100000+i, 100000+2*i-1, 100000+2*i*enddoNUMCMP , ALL !所有实体进行重新编号*DO, i, 1, L2-1CSYS, 0N,100000+2*i-1, L1, i, -1N,100000+2*i, L1-2.5, i, -1*enddo*DO, i, 1, L2-1TYPE, 2 $ MAT , 1 $ REAL , 0ESYS, 0EN, 100000+i, 100000+2*i-1, 100000+2*i*enddoNUMCMP , ALL !所有实体进行重新编号ESEL, S, TYPE, , 2 !选择单元类型号为2的单元EGEN, 25, 100000, ALL, , , , , , , , , , -1,ESEL, S, TYPE, , 2 !选择单元类型号为2的单元CM, SPRING, ELEM !生成一个叫SPRING由单元组成的元件ALLSALLSEL, ALL !选择所有实体NUMMRG, NODE, 0.01, 0.01, LOW !节点合并距离小于0.01则同保留编号底的点NUMCMP, ALL !所有实体进行重新编号*DO, i, 1, 25NSEL, S, LOC, Z, -1*i !从当前集选择一组节点子集深度-1以下ESLN, S !从已选集中选择NSEL, R, TYPE, , 2 !选择单元类型号为2的单元R, i, m1*i*b*h, , , !单位面积内受的力随深度增加而增加EMODIF, ALL, REAL, i, !对已存在单元进行修改*enddoR, 101, 0.6, 0.6, 0.6, 0.6, , !识别号+实常数R, 102, 0.1, 0.1, 0.1, 0.1, ,CMSEL, S, WALLEMODIF, ALL, REAL, 101, !将实常数101组赋给墙CMSEL, S, STRUTEMODIF, ALL, REAL, 102, !赋值给支撑NSEL, S, LOC, Z-25 !约束墙底竖向位移D, ALL, , , , , , UZ, , , , ,CMSEL, S, SPRING !选择土弹簧单元NSLE, S !以下命令从已选弹簧集合中选CMSEL, S, WALL !在已选集中选墙单元NSLE, U !从集合中删除刚选择的单元，即与墙有关的单元D, ALL, , , , , , ALL, , , , , !约束土弹簧单元端点的所有位移ALLSSA VE/SOLUALLSANTYPE, STATIC, NEW !分析模式静力NROPT, FULL !指定计算模式*AFUN, DEG !指定角度单位为度Q=2.0E4 !Q为超载GAMA1=0.9E4 !浮重度FAI1=20 !内摩擦角C1=10e3 !粘聚力m1=1500e3 !比例系数b=1 !单元宽度h=1 !墙体单元高度!无支撑开挖1m 坑内水位-1.0 坑外0NSEL, S, LOC, Z, 0, -0.9 !0到0.9的位置即开挖面以上ESLN, S !以下命令在当前集里选择CMSEL, R, wall !在当前集选wll 单元NSLE, S !以下命令在当前集里选择*GET, ZMIN, NODE, , MNLOC, Z, , , , !墙最浅节点处的位置*GET, ZMAX, NODE, , MXLOC, Z, , , , !墙最深节点处的位置LOCZ1=abs (ZMAX) !取绝对值LOCZ2=abs (ZMIN)KA= (TAN(45.0 - FAI1/2))**2 !主动土压力系数!采用水土分算PA1=(Q + GAMA1*LOCZ1)*KA-2.0*C1*SQRT(KA) !最深主动土压力公式* IF , PA1, LT, 0, THEN !去除小于零值的可能PA1=0*ENDIFPA2=10*1E3*LOCZ1 !水压力计算PA=PA1+PA2 !总应力SA1=(Q+GAMA1*LOCZ2)*KA-2*C1*SQRT(KA)SA2=10*1E3*LOCZ2WA=SA2SA=SA1+SA2pressure=SASLZER=ZMAXSLOPE=(SA-PA)/(LOCZ2-LOCZ1) !增长率SFGRAD, PRES, 0, z, SLZER, SLOPE !沿z方向从-1m开始的面力减少SFE, all, 1, PRES, , -PA, , , !alls!开挖面以下加载NSEL, S, LOC, Z, -1.1, -24.9ESLN, SGMSEL, R, wallNSLE, S*GET, ZMIN, NODE, , MNLOC, Z, , , ,*GET, ZMAX, NODE, , MXLOC, Z, , , ,SLZER=ZMAXSLOPE=0SFGRAD, PRES, O, z, SLZER, SLOPESFE, all, 1, PRES, , -pressure, , , !alls !将集合扩大到全集cmsel, S, strut !选择内支撑全部杀死EKILL, ALLallsNSEL, S, LOC, Z, 0.1, -1.1ESLN, SESEL, R, TYPE, , 2 !选择单元类型号为2的单元EKILL, ALL*do, i, 2, 25z=i-1 !开挖了1m 要在被动区减去相应的土压力K=m1*z*b*hR, i, K, , ,*enddoTIME, 1 !载荷步1allssolveNSEL, S, LOC, Z, 0, -4.9ESLN, SCMSEL, R, wallNSLE, S*GET, ZMIN, NODE, , MNLOC, Z, , , ,*GET, ZMAX, NODE, , MXLOC, Z, , , ,LOCZ1=abs (ZMAX)LOCZ2=abs (ZMIN)KA= (TAN(45.0 - FAI1/2))**2PA1=(Q + GAMA1*LOCZ1)*KA-2.0*C1*SQRT(KA)* IF , PA1, LT, 0, THENPA1=0*ENDIFPA2=10*1E3*LOCZ1PA=PA1+PA2SA1=(Q+GAMA1*LOCZ2)*KA-2*C1*SQRT(KA)SA2=10*1E3*LOCZ2WA=SA2SA=SA1+SA2pressure=SASLZER=ZMAXSLOPE=(SA-PA)/(LOCZ2-LOCZ1)SFGRAD, PRES, 0, z, SLZER, SLOPESFE, all, 1, PRES, , -PA, , ,allsNSEL, S, LOC, Z, -5.1, -24.9ESLN, SGMSEL, R, wallNSLE, S*GET, ZMIN, NODE, , MNLOC, Z, , , ,*GET, ZMAX, NODE, , MXLOC, Z, , , ,SLZER=ZMAXSLOPE=0SFGRAD, PRES, O, z, SLZER, SLOPESFE, all, 1, PRES, , -pressure, , ,allscmsel, s, strutealive, ALLallsNSEL, S, LOC, Z, 0.1, -5.1ESLN, SESEL, R, TYPE, , 2 !选择单元类型号为2的单元EKILL, ALL*do, i, 6, 25 z=i-5K=m1*z*b*h R, i, K, , ,*enddo TIME, 2allssolve。
MATLAB有限元编程03梁单元引言平面纯弯梁单元描述有限元基本格式描述节点位移、节点力位移场：其中，为梁单元内任一点位移，为单元形函数，注意到这里已经用到了“归一化”思想，为以后讲解等参单元做基础。
应变场：其中，为单元的几何矩阵，为所测点以中性层为起点的y方向的坐标。
应力场：其中，为弹性模量，为单元的应力矩阵。
单元刚度矩阵：单元刚度方程：模型描述Matlab编程：% 公众号：易木木响叮当，回复梁单元，即可自动获取。
function k =Beam1D2Node_Stiffness(E,I,L)% 直接组装梁单元刚度k = E*I/(L*L*L)*[12 6*L -12 6*L6*L 4*L*L -6*L 2*L*L-12 -6*L 12 -6*L6*L 2*L*L -6*L 4*L*L];function z = Beam1D2Node_Assembly(KK,k,i,j)% 该函数进行整体刚度矩阵的组装% 输入单元刚度矩阵k，单元的节点编号i、j% 输出整体刚度矩阵KKDOF(1)=2*i-1;DOF(2)=2*i;DOF(3)=2*j-1;DOF(4)=2*j;for n1=1:4for n2=1:4KK(DOF(n1),DOF(n2))= KK(DOF(n1),DOF(n2))+k(n1,n2); endendz=KK;function v = Beam1D2Node_Deflection(x,L,u)% 该函数计算单元内某点的挠度% 输入所测点距梁单元左节点的水平距离x% 输入梁单元的长度L，节点位移列阵ue=x/L;N1=1-3*e*e+2*e*e*e;N2=L(e-2*e*e+e*e*e);N3=3*e*e-2*e*e*e;N4=L(e*e*e-e*e);N=[N1,N2,N3,N4];v=N*u% -----------主程序----------------format compactE = 200*10^9;I = 118.6*10^-6;L1 = 5;L2 = 2.5;% cal element stiffnessk1 = Beam1D2Node_Stiffness(E,I,L1);k2 = Beam1D2Node_Stiffness(E,I,L2);% stiffness assemblyKK = zeros(6,6);KK = Beam1D2Node_Assembly(KK,k1,1,2);KK = Beam1D2Node_Assembly(KK,k2,2,3)% cal displacementk = KK(4:6,4:6);p = [39062;-31250;13021];u = k\p% cal node forceU = [0;0;0;u];P = KK*U% cal RFF = [-62500;-52083;-93750;39062;-31250;13021];RF = P-F最终计算得到的支反力、节点位移值与书中解析解一致，可自行验证。
基于VB和MATLAB混合编程下的弹性地基梁分析
同霄;张笑
【期刊名称】《山西建筑》
【年(卷),期】2010(036)012
【摘要】利用Active X技术,基于VB和MATLAB混合编程,对弹性地基梁受集中力和集中力偶两种情况进行分析,并在MATLAB下绘制弯矩图,解决了VB在数值计算方面的欠缺,弥补MATLAB在人机互动方面的缺陷,为今后进一步分析提供新的
平台.
【总页数】2页(P364-365)
【作者】同霄;张笑
【作者单位】兰州交通大学土木工程学院硕士研究生,甘肃,兰州,730000;长江大学地球科学学院硕士研究生,湖北,荆州,434500
【正文语种】中文
【中图分类】TP393
【相关文献】
1.基于VB与Matlab混合编程的电连接器接触件优化设计 [J], 叶珍珍;徐雷;王强
2.基于COM组件的VB与MATLAB混合编程技术在地形变数据分析中的应用 [J], 胡静;吴云;张燕
3.基于COM组件的VB与MATLAB混合编程技术在地形变数据分析中的应用 [J], 胡静;吴云;张燕
4.基于MATLAB与VB混合编程的麻花钻仿真 [J], 郭天骏;赵鸿生
5.基于VB和Matlab混合编程的离心泵测试数据分析 [J], 王宏伟;王蕊;商德勇;于治福
因版权原因，仅展示原文概要，查看原文内容请购买。
