有限元法的.ppt

《有限元法的.ppt》由会员分享，可在线阅读，更多相关《有限元法的.ppt（58页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。

1、第二章第二章 有限元法的预备知识有限元法的预备知识矩阵、线性代数和MATLAB在数学上，有限元法通过变分法或加权残数法把微分方程化为近似的代数方程组进行求解。如果采用矩阵形式进行中间过程的推导，就能使有限元公式简洁紧凑。为此，我们必须掌握矩阵和线性代数的基本知识，了解矩阵的详细操作，这非常有助于阅读和理解有限元法的基本理论并组织有效的算法 在叙述矩阵、线性代数基本知识的同时，将介绍使用MATLAB实现矩阵定义，运算等基本操作，这有助于增加感性认识，加深理解。矩阵简介A=1.4 0.0 -1.0 -0.2 0.0 1.4 -0.4 0.4 -1.0 -0.4 1.8 0.6 -0.2 0.4 0

2、.6 2.4 A=1.4000 0 -1.0000 -0.2000 0 1.4000 -0.4000 0.4000 -1.0000 -0.4000 1.8000 0.6000 -0.2000 0.4000 0.6000 2.4000在在MATLAB的命令窗口中输入：的命令窗口中输入：MATLAB在命令窗口显示如下：在命令窗口显示如下：表示刚才的输入正确，已经被表示刚才的输入正确，已经被MATLAB接受，下次可以引用矩阵接受，下次可以引用矩阵A元素之间可以用空格元素之间可以用空格或者逗号隔开或者逗号隔开行与行之间可行与行之间可以用回车或者以用回车或者分号隔开分号隔开列向量列向量b可以这样输入：可

3、以这样输入：b=1.0;0.0;0.0;0.0特殊矩阵1.零矩阵MATLAB中中定定义义零零矩矩阵阵十十分分方方便便，比比如如定定义义一一个个43的的零零矩阵的命令是矩阵的命令是 A=zeros(4,3)A=0 0 0 0 0 0 0 0 0 0 0 0MATLAB在命令窗口显示如下：在命令窗口显示如下：命令命令zeros接受可接受任接受可接受任意个参数，一般我们使意个参数，一般我们使用的都是用的都是2维的矩阵，前维的矩阵，前一个参数表示行数，后一个参数表示行数，后一个是列数。一个是列数。2.对称矩阵目前，在目前，在MATLAB中中没有定义对称矩阵的函没有定义对称矩阵的函数，或者说没有对称矩数

4、，或者说没有对称矩阵概念，这是阵概念，这是MATLAB中一个比较中一个比较大的缺憾！大的缺憾！MATLAB的函数的函数transpose是专门对矩阵进行转置操作的，是专门对矩阵进行转置操作的，例如例如 A=-1 2 3 4-3 1A=-1 2 3 4 -3 1 transpose(A)ans=-1 4 2 -3 3 1函数函数transpose可以用一个对应的操作符单引号可以用一个对应的操作符单引号 来代替，例如来代替，例如 Aans=-1 4 2 -3 3 1这里需要指出的是，函数这里需要指出的是，函数transpose与单引号对于实矩阵与单引号对于实矩阵的操作是完全一样的，而对于复的操作是

5、完全一样的，而对于复矩阵则不同。因为复矩阵在数学矩阵则不同。因为复矩阵在数学上的转置定义不仅包括行列互换，上的转置定义不仅包括行列互换，还要执行共轭操作。函数还要执行共轭操作。函数transpose只执行行列互换的操只执行行列互换的操作，而单引号则是数学意义上的作，而单引号则是数学意义上的转置。转置。3.单位矩阵MATLAB中用中用eye来定义单位矩阵，下面的命令定义来定义单位矩阵，下面的命令定义了一个三阶的单位矩阵了一个三阶的单位矩阵 eye(3)ans=1 0 0 0 1 0 0 0 1函数函数eye也能接受多个参数，当它们不相等时，给出的是对角线也能接受多个参数，当它们不相等时，给出的是

6、对角线元素为元素为1其余元素为零的非方矩阵，例如其余元素为零的非方矩阵，例如 eye(5,4)ans=1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1 0 0 0 0为什么为什么不用不用 I？4.带状稀疏矩阵对称带状矩阵的一维存储 在在MATLAB中中，我我们们可可以以用用命命令令sparse定定义义一一个个稀稀疏疏矩阵，比如矩阵，比如 A=sparse(4000,4000)命令窗口中显示命令窗口中显示A=All zero sparse:4000-by-4000用用来来表表示示A是是一一个个40004000的的元元素素全全部部为为零零的的稀稀疏疏矩矩阵。对矩阵元素的引用跟普通的矩

7、阵一样。阵。对矩阵元素的引用跟普通的矩阵一样。矩阵的运算1.相等2.加法与减法MATLAB中矩阵的运算跟单个数值的运算具有相同的形式，比如中矩阵的运算跟单个数值的运算具有相同的形式，比如A=-1 2 3;4-3 1;B=-1 1-1;-3 2-5;C=3-1 4;-2 6 3;A+Bans=-2 3 2 1 -1 -4B+Aans=-2 3 2 1 -1 -4(A+B)+Cans=1 2 6 -1 5 -1A+(B+C)ans=1 2 6 -1 5 -1操作操作方便！方便！数乘矩阵 A=-1 2 3;4-3 1;B=-1 1-1;-2 6 3;k=4;l=2;k*Aans=-4 8 12 16

8、 -12 4 A*kans=-4 8 12 16 -12 4 k*(l*A)ans=-8 16 24 32 -24 8(k*l)*Aans=-8 16 24 32 -24 8(k+l)*Aans=-6 12 18 24 -18 6 k*A+l*Aans=-6 12 18 24 -18 6 k*(A+B)ans=-8 12 8 8 12 16 k*A+k*Bans=-8 12 8 8 12 16操作操作方便！方便！矩阵相乘 矩阵相乘的性质 A=1;2;B=3 4;A*Bans=3 4 6 8 B*Aans=11我们用我们用MATLAB来验证上式，即来验证上式，即 矩阵相乘的性质下面用下面用MAT

9、LAB验证矩阵乘法的分配律。先定义三个随机矩阵验证矩阵乘法的分配律。先定义三个随机矩阵 A=rand(4,3);B=rand(4,3);C=rand(3,2);然后验证乘法的分配律，即然后验证乘法的分配律，即(A+B)*Cans=0.8051 1.8023 0.6336 1.2009 0.5684 0.9928 0.4681 0.9461 A*C+B*Cans=0.8051 1.8023 0.6336 1.2009 0.5684 0.9928 0.4681 0.9461操作操作方便！方便！矩阵不能约去！但是但是矩阵积的转置如果那么逆矩阵如果如果那么那么矩矩阵阵的的求求逆逆对对于于一一般般的的程

10、程序序语语言言必必须须编编写写一一段段复复杂杂的的程程序序，根根据据矩矩阵阵的的特特点点，选选择择合合适适的的算算法法，十十分分麻麻烦烦。而而对对于于MATLAB，则则是是轻轻而而易易举举，一一个个inv函函数数就就足足够够了了。例如，例如，A=rand(4,4)A=0.0153 0.4660 0.2026 0.6813 0.7468 0.4186 0.6721 0.3795 0.4451 0.8462 0.8381 0.8318 0.9318 0.5252 0.0196 0.5028 B=inv(A)B=0.2203 1.3026 -1.1110 0.5563 -5.0783 -5.4552

11、 5.5605 1.7991 0.1578 1.6202 -0.1151 -1.2461 4.8895 3.2202 -3.7443 -0.8724上面的矩阵上面的矩阵B就是就是A的逆矩阵，我们可以通过下面的命令可以验证它们的互逆关系的逆矩阵，我们可以通过下面的命令可以验证它们的互逆关系 A*Bans=1.0000 0 0 0 0 1.0000 0 0 0 0 1.0000 0 0 0 0 1.0000A和和B相乘的结果是单位阵，因此它们是互逆的。相乘的结果是单位阵，因此它们是互逆的。矩阵分块 MATLAB中分块矩阵是通过冒号来实现的。假设中分块矩阵是通过冒号来实现的。假设A是一个是一个55的

12、矩阵，那么的矩阵，那么 A(2,:)表示矩阵表示矩阵A的第二行的第二行 A(:,3)表示矩阵表示矩阵A的第三列的第三列 A(2:4,:)表示由矩阵表示由矩阵A的第的第2，3和和4行元素组成的一个行元素组成的一个35的矩阵，的矩阵，A(:,3:5)表示由矩阵表示由矩阵A的第的第3，4和和5列元素组成的一个列元素组成的一个53的矩阵的矩阵 A(2:4,3:5)表示由矩阵表示由矩阵A的第的第2行到第行到第4行，第行，第3列到第列到第5列元素组成的一个列元素组成的一个33的矩阵。的矩阵。A(1:2:5,:)表示由矩阵表示由矩阵A的第的第1，3和和5行元素组成的一个行元素组成的一个35的矩阵。的矩阵。A

13、=rand(5,5)A=0.7095 0.6822 0.3784 0.8998 0.3420 0.4289 0.3028 0.8600 0.8216 0.2897 0.3046 0.5417 0.8537 0.6449 0.3412 0.1897 0.1509 0.5936 0.8180 0.5341 0.1934 0.6979 0.4966 0.6602 0.7271 A(2,:)ans=0.4289 0.3028 0.8600 0.8216 0.2897 A(:,3)ans=0.3784 0.8600 0.8537 0.5936 0.4966 A(2:4,:)ans=0.4289 0.30

14、28 0.8600 0.8216 0.2897 0.3046 0.5417 0.8537 0.6449 0.3412 0.1897 0.1509 0.5936 0.8180 0.5341A=0.7095 0.6822 0.3784 0.8998 0.3420 0.4289 0.3028 0.8600 0.8216 0.2897 0.3046 0.5417 0.8537 0.6449 0.3412 0.1897 0.1509 0.5936 0.8180 0.5341 0.1934 0.6979 0.4966 0.6602 0.7271 A(:,3:5)ans=0.3784 0.8998 0.34

15、20 0.8600 0.8216 0.2897 0.8537 0.6449 0.3412 0.5936 0.8180 0.5341 0.4966 0.6602 0.7271 A(2:4,3:5)ans=0.8600 0.8216 0.2897 0.8537 0.6449 0.3412 0.5936 0.8180 0.5341 A(1:2:5,:)ans=0.7095 0.6822 0.3784 0.8998 0.3420 0.3046 0.5417 0.8537 0.6449 0.3412 0.1934 0.6979 0.4966 0.6602 0.7271矩阵的迹和行列式MATLAB中用两个

16、函数中用两个函数trace和和det分别求取矩阵的迹和行列式值。分别求取矩阵的迹和行列式值。A=3 4 6 3 0 2 -2 5 0 0 -2 3 0 0 0 4A=3 4 6 3 0 2 -2 5 0 0 -2 3 0 0 0 4 trace(A)ans=7 det(A)ans=-48太简单了！太简单了！MATLAB语言程序设计初步if I=J A(I,J)=2;elseif abs(I-J)=1 A(I,J)=-1;else A(I,J)=0;endIF语句MATLAB与与C/C+十分相十分相似！似！help ifIF statement condition.The general for

17、m of the IF statement is IF expression statements ELSEIF expression statements ELSE statements END The statements are executed if the real part of the expression has all non-zero elements.The ELSE and ELSEIF parts are optional.Zero or more ELSEIF parts can be used as well as nested IFs.The expressio

18、n is usually of the form expr rop expr where rop is=,=,or=.SWITCH语句method=Bilinear;switch lower(method)case linear,bilinear disp(Method is linear);case cubic disp(Method is cubic);case nearest disp(Method is nearest);otherwise disp(Unknown method.);end help switchSwitch among several cases based on

19、expression.The general form of the SWITCH statement is:SWITCH switch_expr CASE case_expr,statement,.,statement CASE case_expr1,case_expr2,case_expr3,.statement,.,statement .OTHERWISE,statement,.,statement END The statements following the first CASE where the switch_expr matches the case_expr are exe

20、cuted.When the case expression is a cell array(as in the second case above),the case_expr matches if any of the elements of the cell array match the switch expression.If none of the case expressions match the switch expression then the OTHERWISE case is executed(if it exists).Only one CASE is execut

21、ed and execution resumes with the statement after the END.The switch_expr can be a scalar or a string.A scalar switch_expr matches a case_expr if switch_expr=case_expr.A string switch_expr matches a case_expr if strcmp(switch_expr,case_expr)returns 1(true).Only the statements between the matching CA

22、SE and the next CASE,OTHERWISE,or END are executed.Unlike C,the SWITCH statement does not fall through(so BREAKs are unnecessary).FOR 语句N=10;for I=1:N for J=1:N A(I,J)=1/(I+J-1);endend help forFOR Repeat statements a specific number of times.The general form of a FOR statement is:FOR variable=expr,s

23、tatement,.,statement END The columns of the expression are stored one at a time in the variable and then the following statements,up to the END,are executed.The expression is often of the form X:Y,in which case its columns are simply scalars.Some examples (assume N has already been assigned a value)

24、.WHILE语句E=0*A;F=E+eye(size(E);N=1;while norm(E+F-E,1)0 E=E+F;F=A*F/N;N=N+1;end help whileWHILE Repeat statements an indefinite number of times.The general form of a WHILE statement is:WHILE expression statements END The statements are executed while the real part of the expression has all non-zero e

25、lements.The expression is usually the result of expr rop expr where rop is=,=,or=.The BREAK statement can be used to terminate the loop prematurely.sinpn1=1.0;for k=2:15 n=2k;sinpn=sqrt(2)/2*sqrt(1-sqrt(1-sinpn12);sinpn1=sinpn;pn=n*sinpn;fprintf(n=%7d pi=%.10fn,n,pn);end一个简单的例子计算圆周率的近似方法如下：计算圆周率的近似方

26、法如下：用用MATLAB实现上述公式的代码为：实现上述公式的代码为：FUNCTION语句function mean,stdev=stat(x)%main function n=length(x);mean=avg(x,n);stdev=sqrt(sum(x-avg(x,n).2)/n);function mean=avg(x,n)%AVG sub-function mean=sum(x)/n;假设硬盘上有文件假设硬盘上有文件stat.m，其中有这样两个函数，其中有这样两个函数函数函数stat可以在其可以在其他他m文件中被调文件中被调用，而用，而avg只能在只能在文件文件stat.m中有中有效。

27、效。弹性力学基础（补充内容）弹性力学基础（补充内容）基本假设：基本假设：（1）理想弹性：物体在撤除外力后完全恢复到初始状态）理想弹性：物体在撤除外力后完全恢复到初始状态（2）连续性：材料是连续的，原子结构不被考虑）连续性：材料是连续的，原子结构不被考虑（3）小变形：变形很小，方程是线性的）小变形：变形很小，方程是线性的补充假设：补充假设：（1）均匀性：物体的各个部分的材料特性是一致的）均匀性：物体的各个部分的材料特性是一致的（2）各向同性：各个方向的材料特性是一致的）各向同性：各个方向的材料特性是一致的在有限元中后面的补充假设其在有限元中后面的补充假设其实是不必要的，因为对于有限实是不必要的，

28、因为对于有限元来说，是否均匀和各向同性元来说，是否均匀和各向同性对有限元实现的难度几乎是一对有限元实现的难度几乎是一样的！样的！弹性力学研究内容弹性力学研究内容位位 移移（3）应应 变变（6）应应 力力（6）几何关系几何关系（6）本构关系本构关系（6）平衡方程平衡方程（3）边界条件边界条件用矩阵表示的基本物理量用矩阵表示的基本物理量位移位移应力应力应变应变用矩阵表示的几何方程用矩阵表示的几何方程用矩阵表示的本构方程用矩阵表示的本构方程用矩阵表示的平衡方程用矩阵表示的平衡方程用矩阵表示的边界条件用矩阵表示的边界条件位移的边界条件位移的边界条件力的边界条件力的边界条件二维：二维：三维：三维：弹性力

29、学控制方程平衡方程几何方程本构关系边界条件最小势能原理可能位移满足边界条件的位移真实位移满足边界条件和控制方程的位移最小势能原理在所有的可能位移中，真实位移使系统的总势能取最小值用位移表示的势能有限元分析过程概述建立力学模型结构离散选择位移模式用结点位移表示的势能应用最小势能原理集合所有单元的平衡方程，建立整个结构的平衡方程1.将各单元的刚度矩阵集合成整个结构的整体刚度矩阵2.将各单元的等效结点力集合成整个结构的等效结点力求解未知结点位移，计算单元内力或应力这里要注意的是，最后求解前要修改前面得到的平衡方程，以考虑结构的边界条件。否则刚度矩阵是奇异的。矩阵导数的概念（补充）1.定义：函数对矩阵的导数例子如果对称如果对称