Matlab 有限元法计算分析程序编写

⎡ k 1,1 ⎢ ⎢k 2 ,1 ⎢ M ⎢ ⎢ k i ,1 ⎢ M ⎢ ⎢k n ,i ⎣
L L L L M O
k 1,i k 2 ,i M
L L 10 20 M O M L L k n ,i
L k 1, n ⎤ ⎧ a 1 ⎫ ⎧ ⎥ L k 2,n ⎥ ⎪ a 1 ⎪ ⎪ ⎪ ⎪ ⎪ O M ⎥⎪ M ⎪ ⎪ ⎪ ⎪ ⎪ ⎥⎨ ⎬ = ⎨ L k i , n ⎥ ⎪ a i ⎪ ⎪β i O M ⎥⎪ M ⎪ ⎪ ⎥⎪ ⎪ ⎪ L k n , n ⎥ ⎪a n ⎪ ⎪ ⎦⎩ ⎭ ⎩
3) 数值、变量 数值采用习惯的十进制表示，可以带小数点或负号，如 3 -99 0.001 9.456 1.3e-3 4.5e33 变量命名规则 变量名、函数名是对字母大小写敏感的 变量名、函数名的第一个字符必须是英文字母，最多可包 含31个字符 变量名中不能包含空格、标点 4) 表达式 MATLAB书写表达式的规则与手写算式几乎完全相同，具 体是 表达式由变量名、运算符合函数和组成 表达式将按与常规相同的优先级自左向右执行运算 优先级的规定：幂运算 > 乘除 > 加减 括号可以改变运算的次序
2) 简单矩阵的输入 (1) 在键盘上输入下列内容 A=[1,2,3;4,5,6;7,8,9] (2) 按【Enter】键，指令被执行 (3) 在指令被执行后，MATLAB指令窗中将显示以下结果 A = ⎡ 1 2 3⎤ 1 2 3 A = ⎢ 4 5 6⎥ ⎥ ⎢ 4 5 6 ⎢7 8 9 ⎥ ⎦ ⎣ 7 8 9 [说明]：在全部键入一个指令行内容后，必须按下【Eenter】键，该指令 才会被执行。 • 直接输入矩阵时，矩阵元素用空格或逗号‘，’分开；矩阵行用“；” 隔离，整个矩阵放在“[]”里。 • 在MATLAB里，不必事先对矩阵维数作任何说明，存储时将自动 配置 • 指令执行后，矩阵A被保存在MATLAB的工作空间中，以备后用。 如果用户不用Clear指令清除它或对它重新赋值，那么该矩阵会一直保 存在工作空间中，直到MATLAB指令窗被关闭。 • MATLAB对变量的大小写敏感。比如本例中的矩阵赋给了变量A， 而不是a。
计算应变和应力
• 对各单元计算 1）按自由度读出单元的节点位移 2）如果必要计算进行坐标转换（整体到局部） 3）计算应变和应力 5）输出结果
补充：刚度矩阵的基本性质
• 刚 阵 发 单பைடு நூலகம்• 刚 • 刚 • 刚 • 刚 为 阵 每 义 结 阵 阵 对 阵 运动 阵 异阵 时 绝对 异 确 阵 说 刚 线 刚 带 阵 结 编 规则 义 义 结构 在 态
计算荷载向量
对各节点集中荷载作如下的计算 a)分解成坐标方向的荷载分量 b)按自由度顺序叠加到荷载向量中 对非节点集中荷载作如下的计算 a)在单元局部坐标系下计算等效节点荷载 b)作坐标变换（如果需要） c)按自由度顺序叠加到荷载向量中
引入边界条件
• 最简单的方法是改变主元为十分大的值 （采用这种方法时，可以在刚度方程计算 过程中实施） • 修改刚度方程
a i = βi
⎡ k 1,1 ⎢ ⎢k 2 ,1 ⎢ M ⎢ ⎢ 0 ⎢ M ⎢ ⎢ k n ,i ⎣ L L 0 L k 1, n ⎤ ⎧ a 1 ⎫ ⎧ P1 ⎫ ⎥ L L 0 L k 2 , n ⎥ ⎪ a 1 ⎪ ⎪P2 ⎪ ⎪ ⎪ ⎪ ⎪ M O M O M ⎥⎪ M ⎪ ⎪ M ⎪ ⎪ ⎪ ⎪ ⎪ ⎥⎨ ⎬ = ⎨ ⎬ L L 1 L 0 ⎥ ⎪ a i ⎪ ⎪β i ⎪ M O M O M ⎥⎪ M ⎪ ⎪ M ⎪ ⎥⎪ ⎪ ⎪ ⎪ L L 0 L k n , n ⎥ ⎪a n ⎪ ⎪Pn ⎪ ⎦⎩ ⎭ ⎩ ⎭
MATLAB的使用方法
1) 最简单的计算器使用法 求[12＋2×(7－4)]÷32的算术运算结果 （1）用键盘在MATLAB指令窗中输入一下内容 (12+2*(7-4))/3^2 （2）在上述表达式输入完成后，按【Enter】键，该指令被执行 （3）在指令执行后，MATLAB指令窗中将显示一下内容 ans = 2 [说明] 加 + 减 乘 * 除 / 或 \ (这两个符号对于数组有不同的含义) 幂 ^ “ans”是answer的缩写，其含义是运算答案，它是MATLAB的一个默 认变量
6) M函数文件 与命令文件不同，函数文件从外界只能看到传给它的输入 量和送出来的计算结果，而内部运作是看不见的。它的特 点是 (1)从形式上看，与命令文件不同，函数文件的第一行总是以 “function”引导的“函数申明行”。 (2)从运行上看，与命令文件运行不同，每当函数文件运行， MATLAB就会专门为它开辟临时工作空间，所有中间变量 都存放在函数工作空间中，当执行完文件最后一条指令和 遇到return时，就结束该函数文件的执行，同时该临时函数 工作空间及其所有的中间变量立即被清除。 (3)对于函数文件中的变量，如果不作特别说明，默认为临时 局部变量，这些临时变量就存放在函数的临时工作空间中， 当函数结束时他们被立即清除。与之相对应的是全局变量， 他们是通过global指令进行特别申明，这些全局变量可被几 个不同的函数共享。 • 函数文件的编辑也可用MATLAB editor/debugger。
⎡ k 1,1 ⎢ ⎢k 2 ,1 ⎢ M ⎢ ⎢ k i ,1 ⎢ M ⎢ ⎢k n ,i ⎣
L L k 1,i L L k 2 ,i M O M L L k i ,i M M M L L k n ,i
L k 1, n ⎤ ⎧ a 1 ⎫ ⎧ P1 ⎫ ⎥ L k 2 , n ⎥ ⎪ a 1 ⎪ ⎪P2 ⎪ ⎪ ⎪ ⎪ ⎪ M M ⎥⎪ M ⎪ ⎪ M ⎪ ⎪ ⎪ ⎪ ⎪ ⎥⎨ ⎬ = ⎨ ⎬ L k i , n ⎥ ⎪ a i ⎪ ⎪ Pi ⎪ O M ⎥⎪ M ⎪ ⎪ M ⎪ ⎥⎪ ⎪ ⎪ ⎪ L k n , n ⎥ ⎪a n ⎪ ⎪Pn ⎪ ⎦⎩ ⎭ ⎩ ⎭
平面问题的有限元程序 function plane % This function analyze the problem of plane stress use triangle elements % the variables lists: % node ------- node definition % element ----- element definition % material ---- material definition % K ----------- global stiffness matrix -| % f ----------- global node force vector | [K]{delta}={f} % delta ------- global node displacement vector -| % bc ---------- boundary condition % nf ---------- node force % file_in ----- input file name % file_out ---- output file name % fid_in ------ input file ID % fid_out ----- output file ID % global node element material K file_in = input('input file name: ', 's' ) ; file_out = input( 'output file name: ', 's') ;
fid_in = fopen( file_in, 'r' ) ; % open input file node_number = fscanf( fid_in, '%d', 1 ) ; %read node number node = zeros( node_number, 2 ) ; for i=1:1:node_number nn = fscanf( fid_in, '%d', 1 ) ; node( i, : ) = fscanf( fid_in, '%f', [1,2] ) ; % read node definition end element_number = fscanf( fid_in, '%d', 1 ) ; % read element number element = zeros( element_number, 4 ) ; for i=1:1:element_number ne = fscanf( fid_in, '%d', 1 ) ; element( i, : ) = fscanf( fid_in, '%d', [1,4] ) ; % read element definition end
material_number = fscanf( fid_in, '%d', 1 ) ; % read material number material = zeros( material_number, 3 ) ; for i=1:1:material_number nm = fscanf( fid_in, '%d', 1 ) ; material( i, : ) = fscanf( fid_in, '%f', [1,3] ) ; % read materials definition end bc_number = fscanf( fid_in, '%d', 1 ) ; % read boundary conditions number bc = zeros( bc_number, 3 ) ; for i=1:1:bc_number % read boundary condition definition bc( i, 1 ) = fscanf( fid_in, '%d', 1 ) ; bc( i, 2 ) = fscanf( fid_in, '%d', 1 ) ; bc( i, 3 ) = fscanf( fid_in, '%f', 1 ) ; end
5) 命令文件 • 对于比较简单的问题和一次性的问题，可以通过 指令窗中直接输入一组指令去求解。但是当要解决 问题比较较复杂时，或当一组指令通过改变少量参 数就可以反复使用去解决不同问题时，直接在指令 窗中输入指令的方法就显得繁琐累赘。设计命令文 件就是解决这个矛盾的。 • 所谓命令文件是指：该文件中的指令形式和前后 位置，与解决同一个问题时在指令窗中输入的那组 指令没有任何区别；MATLAB在运行这个文件时， 只是简单地从文件中读取每一条指令，送到 MATLAB中去执行；与在指令窗中直接运行指令一 样，命令文件运行产生的变量也都是驻留在 MATLAB基本工作空间中。
阵对 阵 对 阵
MATLAB概要
• 对于有限元初学者来说，自己动手编写一个 有限元程序是学习、理解和掌握有限元法的 一条捷径。传统的有限元程序设计大多采用 FORTRAN语言，虽然FORTRAN语言在数值 计算方面享有盛誉，但是它的数据类型相对 单一，编写程序有些难度，而MATLAB则相 对容易。
有限元法计算分析程序编写
结构参数输入，包括
1）节点坐标值 2）单元类型以及连接信息 3）各单元的弹性模量、截面积（厚度）等 4）荷载形式以及作用位置、作用方向、荷载值 5）约束条件 6）输出信息
m j
对节点和单元分别编号 每个节点的自由度根据 节点号计算得到
i
y
o
x
计算结构的刚度矩阵
对各单元作如下的计算 a)计算单元刚度矩阵 b)计算坐标转换矩阵（如果需要） c)作坐标转换计算（如果需要） d)按自由度顺序叠加到总刚度矩阵中
⎫ ⎪ P2 ⎪ ⎪ M ⎪ ⎬ × 10 20 ⎪ ⎪ M ⎪ ⎪ Pn ⎭ P1
解方程（线性方程组）得到位移
⎧K 11a 1 + K 12 a 2 + L K 1, n −1a n −1 + K 1, n a n = P1 ⎪ ⎪K 21a 1 + K 22 a 2 + L K 2 , n −1a n −1 + K 2 , n a n = P2 ⎪ ⎨LL ⎪K a + K n −1,1 1 n −1, 2 a 2 + L K n −1, n −1a n −1 + K n −1, n a n = Pn −1 ⎪ ⎪K n ,1a 1 + K n , 2 a 2 + L K n , n −1a n −1 + K n , n a n = Pn ⎩