matlab有限元分析实例
MATLAB有限元分析与应用精选全文完整版
%SpringElementForces This function returns the element nodal force
%
vector given the element stiffness matrix k
%
and the element nodal displacement vector u.
2019/11/28
§2-1 弹簧元
u1=U(1:2); f1=SpringElementForces(k1,u1);
f1 = -15.0000 15.0000
u2=U(2:3); f2=SpringElementForces(k2,u2);
f2 = -15.0000 15.0000
12
§3-1 弹簧元
%
modulus of elasticity E, cross-sectional
%
area A, and length L. The size of the
%
element stiffness matrix is 2 x 2.
y = [E*A/L -E*A/L ; -E*A/L E*A/L];
2019/11/28
3.1 单元刚度矩阵的形成
function y = SpringElementStiffness(k)
%SpringElementStiffness This function returns the element stiffness %matrix for a spring with stiffness k. %The size of the element stiffness matrix is 2 x 2.
matlab有限元三角形单元编程
matlab有限元三角形单元编程
在MATLAB中进行有限元分析,可以使用其自带的有限元分析工具箱(FEATool)进行编程。
以下是一个简单的例子,演示如何使用三角形单元进行有限元分析:
1. 打开MATLAB,进入FEATool环境。
2. 创建新的有限元模型。
选择“Model”菜单下的“New Model”选项,进入“Model Builder”界面。
3. 在“Model Builder”界面中,选择“2D Triangle”单元类型,并在绘图区域中绘制出三角形网格。
4. 在“Model Builder”界面中,设置材料属性、边界条件和载荷等参数。
5. 运行有限元分析。
选择“Model”菜单下的“Solve”选项,进行有限元求解。
6. 查看结果。
选择“Model”菜单下的“Results”选项,可以查看位移、应力、应变等结果。
以上是一个简单的例子,演示了如何使用三角形单元进行有限元分析。
在实际应用中,还需要根据具体问题进行详细的建模和计算。
Matlab有限元分析操作基础
Matlab有限元分析操作基础Matlab 有限元分析20140226为了用Matlab 进行有限元分析,首先要学会Matlab 基本操作,还要学会使用Matlab 进行有限元分析的基本操作。
1. 复习:上节课分析了弹簧系统x 推导了系统刚度矩阵1122121200k k k k k k k k ----+??2. Matlab有限元分析的基本操作(1)单元划分(选择何种单元,分成多少个单元,标号)(2)构造单元刚度矩阵(列出…)(3)组装系统刚度矩阵(集成整体刚度矩阵)(4)引入边界条件(消除冗余方程)(5)解方程(6)后处理(扩展计算)3. Matlab有限元分析实战【实例1】分析:步骤一:单元划分步骤二:构造单元刚度矩阵>>k1=SpringElementStiffness(100) >>…?步骤三:构造系统刚度矩阵a) 分析SpringAssemble库函数function y = SpringAssemble(K,k,i,j)% This function assembles the element stiffness% matrix k of the spring with nodes i and j into the % global stiffness matrix K.% function returns the global stiffness matrix K% after the element stiffness matrix k is assembled. K(i,i) = K(i,i) + k(1,1);K(i,j) = K(i,j) + k(1,2);K(j,i) = K(j,i) + k(2,1);K(j,j) = K(j,j) + k(2,2);y = K;b) K是多大矩阵?今天的系统刚度矩阵是什么?因为11221212k kk kk k k k----+所以10001000200200 100200300----c) K=SpringAssemble(K,k1,1,2) function y = SpringAssemble(K,k,i,j) K(i,i) = K(i,i) + k(1,1);K(i,j) = K(i,j) + k(1,2);K(j,i) = K(j,i) + k(2,1);K(j,j) = K(j,j) + k(2,2); 1100100100100k -??=??-??10010001001000000K -??=-??K=SpringAssemble(K,k2,2,3) 1200200200200k -??10010001003002000200200K -??=---??1001000100010010030020002002000200200100200300----≠----步骤四:引入边界条件,消除冗余方程>>k=K(2:3,2:3)%构造不含冗余的方程>>f=[0;15]%构造外力列阵步骤五:解方程引例:已知1212u 31u u u +=??-=?,求 12u u 和解:类似求解KU=F ,输入下列Matlab 命令:>> K=[1 1;1,-1]>> F=[3;1]>> U=inv(K)*F>> U=K \F(继续弹簧系统求解)>>u=k \f %使用高斯消去法求解>>U=[0 ; u]%构造原方程组>>F=K*U %求出所有外力,含多余计算步骤六:后处理、扩展计算>>u1=[0;U(2)]%构造单元位移>>f1=SpringElementForces(k1,u1)%求单元1内力>>u2=[U(2) ; U(3)]%构造单元2位移>>f2=SpringElementForces(k2,u2)%求单元2内力4. 总结clck1=SpringElementStiffness(100)%创建单元刚度矩阵 1 k2=SpringElementStiffness(200)%创建单元刚度矩阵 2 K=zeros(3,3)%创建空白整体刚度矩阵K=SpringAssemble(K,k1,1,2)%按节点装入单元矩阵 1 K=SpringAssemble(K,k2,2,3)%按节点装入单元矩阵2 k=K(2:3,2:3)%构造不含冗余的方程f=[0;15]%构造外力列阵u=k\f%使用高斯消去法求解U=[0 ; u]%构造系统节点位移列阵F=K*U%求出所有外力,含多余计算u1=[0;U(2)]%构造单元位移f1=SpringElementForces(k1,u1)%求单元1内力u2=[U(2) ; U(3)]%构造单元2位移f2=SpringElementForces(k2,u2)%求单元2内力5. 练习1 Danyi 132 dan 34 3dan 35 4dan 35 dan5 54 dan6 42。
悬臂梁MATLAB有限元算例注释
用有限元法对悬臂梁分析的算例算例:如下图所示的悬臂梁,受均布载荷q =1N /mm 2作用。
E =2.1×105N /mm 2, μ=0.3厚度h =10mm 。
现用有限元法分析其位移及应力。
梁可视为平面应力状态,先按图示尺寸划分为均匀的三角形网格,共有8×10=80个单元,5×ll =55个节点,坐标轴以及单元与节点的编号如图。
将均布载荷分配到各相应节点上,把有约束的节点5l 、52、53、54、55视作固定铰链,建立如图所示的离散化计算模型。
程序计算框图:(续左)(接右)开 始 输入材料参数 计算具有代表性的单元刚阵 K<=0 将各单元刚阵按整体编号集成到整体刚阵 处理根部约束,修改【K 】【Q 】 求解[K][δ]=[Q] 整理[δ] 并画图计算单元应力,并输出结束程序中的函数功能介绍及源代码1. LinearTriangleElementStiffness(E,NU,t,xi,yi,xj,yj,xm,ym)――该函数用于计算平面应力情况下弹性模量为E、泊松比为NU、厚度为t、第一个节点坐标为(xi,yi)、第二个节点坐标为(xj,yj)、第三个节点坐标为(xm,ym)时的线性三角形元的单元刚度矩阵.该函数返回6×6的单位刚度矩阵k.2. LinearTriangleAssemble(K,k,i,j,m)――该函数将连接节点i,j,m的线性三角形元的单元刚度矩阵k集成到整体刚度矩阵K。
每集成一个单元,该函数都将返回2N×2N的整体刚度矩阵K.3. LinearTriangleElementStresses(E,NU,t,xi,yi,xj,yj,xm,ym,u)-- 该函数计算在平面应力情况下弹性模量为E、泊松比为NU、厚度为t、第一个节点坐标为(xi,yi)第二个节点坐标为(xj,yj)、第三个节点坐标为(xm,ym)以及单元位移矢量为u时的单元应力。
运用MATLAB对桁架单元进行有限元分析
np=3; ne=2; nload=1; nb=4; nu=0;% np为节点数,ne为单元数,nload为外力数,nb为约束数,nu 为泊松比np2=np*2; np3=np2-nb;% np2为不受约束时自由度是节点的两倍 ,np3为不受约束的节点自由度个数K=sym(zeros(np2,np2));% 定义受整体刚度空数组F=sym(zeros(np2,1));% 定义节点外力空数组KK=sym(zeros(np3,np3));% 预置自由度刚度空数组FF=sym(zeros(np3,1));% 预置自由度外力空数组syms h A P E L% 定义未知正常数为变量xy=[0,h;sqrt(3)*h,0;0,0];% 节点横纵坐标数组AA=[A;A];% 单元杆件面积load=[2,2,-P];% 荷载数组bound=[1,1,0;1,2,0;3,1,0;3,2,0];% 边界条件数组IJ=[1,2;3,2];% 各杆单元节点编号DW=zeros(1,4);for ie =1:neip=IJ(ie,1);jp=IJ(ie,2);DW(1)=ip*2-1; DW(2)=ip*2; DW(3)=jp*2-1; DW(4)=jp*2; % 给单元节点横纵方向编号x1=xy(ip,1); x2=xy(jp,1); y1=xy(ip,2); y2=xy(jp,2); % 杆单元节点x,y坐标L=sqrt((x2-x1)^2+(y2-y1)^2);% 杆单元长度c=(x2-x1)/L; s=(y2-y1)/L;% c为余弦, s为正弦T=[c,s,0,0;-s,c,0,0;0,0,c,s;0,0,-s,c];% 坐标转换矩阵(将局部坐标转换为整体坐标)A1=AA(ie);ke=[E*A1/L,0,-E*A1/L,0;0,0,0,0;-E*A1/L,0,E*A1/L,0;0,0,0,0];k=T.'*ke*T;% 将局部坐标下单元刚度转换为整体坐标下单元刚度 (转置后面加.可以去掉结果中的虚数)for i =1:4i1=DW(i);for j =1:4j1=DW(j);K(i1,j1)=K(i1,j1)+k(i,j);% 集成整体刚度矩阵endendendfor i =1:nloadi1=(load(i,1)-1)*2+load(i,2);F(i1,1)=F(i1,1)+load(i,3);% 将荷载按节点方向代码累加,计入外力荷载endFuu=sym(zeros(np2,1));% 节点位移NR=zeros(np2,1);for i =1:nbi1=(bound(i,1)-1)*2+bound(i,2);NR(i1)=-i1;% 将有约束的节点的横纵节点编号挑出来变为负数uu(i1)=bound(i,3);% 有约束的四个方向位移为0endj=0;for i =1:np2i1=NR(i);if i1==0j=j+1;NR(i)=j;endend% 此循环目的是为了将没有约束的节点以及方向挑出来for i =1:np2i1=NR(i);if i1>0FF(i1)=F(i);end% 将没有约束的节点方向的外力储存下来for j =1:np2j1=NR(j);if (i1>0 & j1>0)KK(i1,j1)=K(i,j);end% 将没有约束的节点的刚度储存下来endendKKFFfor i =1:np2i1=NR(i);for j =1:np2j1=NR(j);if (i1>0 & j1<0)jj1=abs(j1);FF(i1)=FF(i1)-K(i,jj1)*uu(jj1);end% 由于有约束方向的位移为0,所以无位移的地方FF[i1]不变endendu=sym(zeros(np2,1));u=KK\FF% 求解线性方程for i =1:np2i1=NR(i);if i1>0uu(i)=u(i1);end% 将求解的无约束方向的位移放进总的位移数组相应的位置中去enduudisp('位移')for ip =1:npstr=[ip;uu(ip*2-1);uu(ip*2)];disp(str)% 输出节点位移endfor ie =1:neip=IJ(ie,1);jp=IJ(ie,2);DW(1)=ip*2-1; DW(2)=ip*2; DW(3)=jp*2-1; DW(4)=jp*2; x1=xy(ip,1); x2=xy(jp,1); y1=xy(ip,2); y2=xy(jp,2);L=sqrt((x2-x1)^2+(y2-y1)^2);c=(x2-x1)/L; s=(y2-y1)/L;T=[c,s,0,0;-s,c,0,0;0,0,c,s;0,0,-s,c];A1=AA(ie);ke=[E*A1/L,0,-E*A1/L,0;0,0,0,0;-E*A1/L,0,E*A1/L,0; 0,0,0,0];d=sym(zeros(4,1));for i=1:4d(i)=uu(i);endde=T*d;% 将整体坐标位移转换为局部坐标位移Fe=ke*de;% 杆单元力for j=1:4N=Fe((ie-1)*2+1);% 杆单元轴力sigma=N/A1;enddisp('单元,轴力,应力')str=[ie;N;sigma];disp(str)endF3=K*uu;R=F3-F;disp('支反力');for i =1:nbi1=(bound(i,1)-1)*2+bound(i,2);str=[bound(i,1);bound(i,3);simplify(R(i1))]; disp(str)end。
matlab编程平面梁问题有限元分析
前处理程序clear;clcXY = [1, 0, 02, 2, 03, 2, -2]; %XY为N行3列,N为节点总数ELB = [1, 1, 2, 12, 2, 3, 1]; %ELB为MB行4列,MB为单元总数b = 4.5e-2;h = 2e-1; %截面宽b和高hA1 = b*h;IZ1 = 3e-5; %惯性矩EAIZ = [200e9, A1, IZ1]; %弹性模量与截面积和惯性矩ELPQ = [1, 1, -2e4, 22, 1, -1e4, 2]; %ELPQ第一列为受载荷单元号,第二列为长度C(见下表),三列为载荷大小,四列为载荷类型SU = [1,02, 03,0]; %SU第一列为已知节点位移对应的方程号,第2列为节点位移数值,约束即等于0子程序1function [PO, ii, jj] = dxjd(ELB, XY, ELPQ1)PO = zeros(6,1);k = ELPQ1(1);jj = ELB(k, 3);dxy = [XY(ii, 2), XY(ii, 3)XY(jj, 2), XY(jj, 3)];DY = dxy(2,2) - dxy(1,2);DX = dxy(2,1) - dxy(1,1);L = sqrt(DX^2 + DY^2);C = ELPQ1(2);Q = ELPQ1(3);type = ELPQ1(4);switch typecase 1PO(2) = PO(2) + 0.5*Q*C*(2-2*C^2/L^2 + C^3/L^3);PO(3) = PO(3) + Q*C^2*(6-8*C/L+3*C^2/L^2)/12;PO(5) = PO(5)+Q*C-PO(2);PO(6) = PO(6) - Q*C^3*(4-3*C/L)/12/L;case 2D = L - C;PO(2) = PO(2) + Q*(L+2*C)*D^2/L^3;PO(3) = PO(3) + Q*C*D^2/L^2;PO(5) = PO(5) + Q-PO(2);PO(6) = PO(6) + Q*D*C^2/L^2;case 3D = L - C;PO(1) = PO(1) + Q*D/L;PO(4) = PO(4) + Q*C/L;case 4PO(2) = PO(2) + 7*Q*L/20;PO(3) = PO(3) + Q*L^2/20;PO(5) = PO(5) + 3*Q*L/20;PO(6) = PO(6) + Q*L^2/30;end子程序2function [KE, T] = ketb(dxy, E, A, IZ)DY = dxy(2,2) - dxy(1,2);DX = dxy(2,1) - dxy(1,1);L = sqrt(DX^2 + DY^2);S = DY/L;C = DX/L;a1 = IZ/L;a2 = a1/L;a3 = E/L;KE = a3*[A, 0, 0, -A, 0, 0 0, 12*a2, 6*a1, 0, -12*a2, 6*a10, 6*a1, 4*IZ, 0, -6*a1, 2*IZ-A, 0, 0, A, 0, 00, -12*a2, -6*a1, 0, 12*a2, -6*a10, 6*a1, 2*IZ, 0, -6*a1, 4*IZ]; t = [C, S, 0; -S, C, 0; 0, 0, 1];t1 = zeros(3,3);T = [t, t1; t1, t];end子程序function KZ = kztb(XY, ELB, EAIZ)[N,M] = size(XY);KZ =zeros(9, 9);[MB, m] = size(ELB);for k = 1 :2ii = ELB(k, 2);jj = ELB(k, 3);LTB = ELB(k, 4);dxy = [XY(ii,2), XY(ii,3)XY(jj,2), XY(jj,3)];E = EAIZ(LTB, 1);A = EAIZ(LTB, 2);IZ = EAIZ(LTB, 3);[KE, T] = ketb(dxy, E, A, IZ);CN = [3*ii-2, 3*ii-1, 3*ii, 3*jj-2, 3*jj-1, 3*jj];KE = (T')*KE*T;for i = 1:6for j = 1:6KZ(CN(i), CN(j)) = KZ(CN(i), CN(j)) + KE(i, j);endendend子程序function U = weiyi(KZ, P, SU)[LR, m] = size(SU);for k =1:LRi = SU(k, 1);P(i) = 1e10*SU(k, 2);KZ(i,i) = 1e10*KZ(i, i);endU = KZ\P;end子程序function P = ydx(XY, ELB, EAIZ, ELPQ)[N, m] = size(XY)P = zeros(3*N, 1)[LPQ, m] = size(ELPQ)for j = 1:LPQELPQ1 = ELPQ(j, :)k = ELPQ(j, 1)[PO, ii, jj] = dxjd(ELB, XY, ELPQ1)LTB = ELB(k, 4);dxy = [XY(ii,2), XY(ii,3)XY(jj,2), XY(jj, 3)];E = EAIZ(LTB, 1);A = EAIZ(LTB, 2);IZ = EAIZ(LTB, 3);[KE, T] = ketb(dxy, E, A, IZ );PO = (T')*PO;CN = [3*ii-2, 3*ii-1, 3*ii, 3*jj-2, 3*jj-1, 3*jj]for i = 1:6P(CN(i)) = P(CN(i)) + PO(i);endend后处理程序clc;[N, M] = size(XY);DU = zeros(N,4);for i = 1:NDU(i,1) = i;DU(i,2) = U(3*i-2);DU(i,3) = U(3*i-1);DU(i,4) = U(3*i);end[MB, M] = size(ELB);for i = 1:MBP1(i, 1) = i;end'节点号水平位移竖直位移转角'DU主程序>>GJINPUT;>>KZ = kztb(XY, ELB,EAIZ);>>P = ydx(XY, ELB, EAIZ, ELPQ);>>U = weiyi(KZ, P, SU);>>GJOUTPUT;运行结果:ans ='节点号水平位移竖直位移转角' DU =1.0000 -0.0000 -0.0000 -0.00002.0000 -0.0000 -0.0111 -0.01003.0000 -0.0231 -0.0111 -0.0125。
matlab桁架结构有限元计算
matlab桁架结构有限元计算摘要:一、引言二、MATLAB 在桁架结构有限元计算中的应用1.桁架结构的离散化与编号2.组装刚度矩阵3.求解器求解4.后处理三、MATLAB 桁架有限元分析实例1.空间桁架有限元分析2.基于MATLAB 的三维桁架有限元分析3.五杆桁架有限元分析四、用MATLAB 编写程序求解桁架结构内力问题1.静定结构绘制简图2.计算结构力学3.机动分析五、结论正文:一、引言在工程领域中,桁架结构由于其优越的力学性能和简单的结构形式,被广泛应用于桥梁、塔架等大型建筑结构中。
为了确保桁架结构的安全性和稳定性,对其进行有限元分析是非常必要的。
MATLAB 作为一种强大的数学软件,可以方便地进行桁架结构的有限元计算。
本文将介绍如何使用MATLAB 进行桁架结构有限元计算。
二、MATLAB 在桁架结构有限元计算中的应用1.桁架结构的离散化与编号在进行有限元分析之前,首先需要对桁架结构进行离散化处理,将连续的桁架杆件划分为有限个小杆件。
同时,为了方便后续计算,需要对各个杆件和节点进行编号。
2.组装刚度矩阵根据桁架结构的几何参数和材料性能,可以计算出各个杆件的刚度矩阵。
将这些刚度矩阵组装成一个总的刚度矩阵,用于描述整个桁架结构的刚度特性。
3.求解器求解利用MATLAB 的求解器,可以对桁架结构进行有限元求解。
求解器会根据刚度矩阵和施加的边界条件,计算出节点的位移和单元的应力。
4.后处理在求解完成后,需要对计算结果进行后处理。
这包括对计算结果的保存、可视化以及结果的验证等。
三、MATLAB 桁架有限元分析实例1.空间桁架有限元分析例如,可以针对一个空间桁架结构,使用MATLAB 进行有限元分析。
假设该桁架结构由铝制成的垂直和水平部分和钢制成的对角桁架构件组成。
结构承受荷载P,需要计算节点位移、单元应力以及支反力。
2.基于MATLAB 的三维桁架有限元分析还可以利用MATLAB 进行三维桁架有限元分析。
matlab有限元分析实例
matlab有限元分析实例问题描述:考虑一平面有界区域,设其边界为[。
我们求解泊松方程之狄利克雷边值问题。
问题的强形式为一椭圆型偏微分方程当之几何形状稍复杂时,一般无法求得其解析解。
我们可应用有限元法来求其数值解。
通常我们先写出该问题的抽象弱形式:求使得其中为检验函数为一适当索伯列夫空间(在本例中也是一希尔伯特空间)为一双线性型为一线性型。
其具体表达式为有限元空间离散我们采用最简单的二维单元离散单元即三节点线性三角形单元,其插值基函数(即形函数为一次多项式。
有限元之核心思想为:使用离散的函数空间来分片逼近连续的函数空间。
于是所求之近似解可写成基函数之线性组合其中为待求系数,常称为自由度。
将此近似解之表达式代入前述问题之弱形式,并取检验函数,可得写成矩阵形式,便成为常见之有限元方程由于历史原因,通常采用固体力学中的习惯命名:[公式] 为刚度矩阵,[公式] 为自由度向量,[公式] 为载荷向量。
由于空间已分片离散,上面的有限元方程只在各单元内部成立。
为了求解方便,通常我们将所有单元的有限元方程连立起来求解,于是需要将各单元之刚度矩阵组装成总体刚度矩阵。
求解器MATLAB 代码及简释本求解器之十行代码如下:function u=fem(nds,els,bcs)nnd=size(nds,1); u=zeros(nnd,1); K=zeros(nnd,nnd); f=zeros(nnd,1);for j=1:size(els,1)K(els(j,:),els(j,:))=K(els(j,:),els(j,:))+stima(nds(els(j,:),:));f(els(j,:))=f(els(j,:))+ones(3,1)*det([1,1,1;nds(els(j,:),:)'])/ 6;endfreends=setdiff(1:nnd,bcs);u(freends)=K(freends,freends)\f(freends);function stima=stima(vertices)ndim=size(vertices,2);J=[ones(1,ndim+1);vertices'];B=J\[zeros(1,ndim);eye(ndim)];stima=det(J).*B*B'/prod(1:ndim);输入数据格式:唯一的一个函数stima 用来计算各单元刚度矩阵。
一维梁的MATLAB有限元法分析
一维梁的M A T L A B有限元法分析问题如下:梁A B在A和B两端固定,中间点表示为C,在中间区域承受均匀分布的载荷q,如图所示。
梁A B的抗弯刚度为E I。
1,使用R i t z法确定点C处的位移和弯矩,并讨论随着包含更多基本函数的准确性。
提示:根据偏转曲线的形状,可以选择基函数的形式,其中系数,并且应该由点A或B处的边界条件确定。
2,采用一维有限元法解决问题,并讨论网格越细时的准确性。
提示:使用1-D梁单元。
1R i t z法一维欧拉-伯努利梁的势能如下:设选择基函数,容易看出基函数满足边界条件设,代入势能表达式得到由于三角函数是正交函数系,所以得到令q=10N/m m,E=200000M P a,I=10000,L=200m m在M A T L A B中计算A k前十项得到A=-0.008400754770396-0.000320811945459-0.000049922673823 -0.000020050746591-0.000009258470170-0.000003960641302 -0.000001943426801-0.000001253171662-0.000000837688763 -0.000000513299113计算C点位移,使用1-10个试函数结果如下:0.0084007547703960.0084007547703960.0085006001180410.0085006001180410.0085191170583800.0085191170583800.0085230039119830.0085230039119830.0085246792895100.008524679289510计算C点弯矩,,使用1-10个试函数结果如下:0.8291212625436840.7024697829907620.746814316705690 0.7151514468174600.7379958063005830.723923419683592 0.7333220380018130.7254063205297550.732103122461494 0.727037063279377可以看到,位移收敛是很快的,弯矩收敛速度慢于位移。
matlab有限元分析实例
求解器之使用
可以使用文本文件存储输入数据,用MATLAB 提供的 load 命令载入。仅需下面三行代码即可完成读取输入数据,调用求解器求解,绘制结果云图诸个步骤。
ndim=size(vertices,2); J=[ones(1,ndim+1);vertices'];
B=J\[zeros(1,ndim);eye(ndim)]; stima=det(J).*B*B'/prod(1:ndim);
输入数据格式:
nds:节点列表矩阵。对于 [公式] 个节点,为一 [公式] 矩阵。第一列为 [公式] 坐标值,第二列为 [公式] 坐标值。
其中[公式] 为检验函数,[公式] 为一适当索伯列夫空间(在本例中也是一希尔伯特空间)。[公式] 为一双线性型,[公式] 为一线性型。其具体表达式为
[公式]
有限元空间离散
我们采用最简单的二维[公式] 单元离散 [公式]。[公式] 单元即三节点线性三角形单元,其插值基函数(即形函数) [公式] 为一次多项式。有限元之核心思想为:使用离散的函数空间 [公式] 来分片逼近连续的函数空间 [公式]。于是所求之近似解 [公式] 可写成基函数之线性组合
[公式]
其中[公式] 为待求系数,常称为自由度。将此近似解之表达式代入前述问题之弱形式,并取检验函数 [公式],可得
[公式]
写成矩阵形式,便成为常见之有限元方程
[公式]
由于历史原因,通常采用固体力学中的习惯命名:[公式] 为刚度矩阵,[公式] 为自由度向量,[公式] 为载荷向量。
matlab有限元分析实例
有限元:
在数学中,有限元法是一种为求解偏微分方程边值问题近似解的数值技术。
求解时对整个问题区域进行分解,每个子区域都成为简单的部分,这种简单部分就称作有限元。
它通过变分方法,使得误差函数达到最小值并产生稳定解。
类比于连接多段微小直线逼近圆的思想,有限元法包含了一切可能的方法,这些方法将许多被称为有限元的小区域上的简单方程联系起来,并用其去估计更大区域上的复杂方程。
它将求解域看成是由许多称为有限元的小的互连子域组成,对每一单元假定一个合适的(较简单的)近似解,然后推导求解这个域总的满足条件(如结构的平衡条件),从而得到问题的解。
这个解不是准确解,而是近似解,因为实际问题被较简单的问题所代替。
MATLAB有限元分析与应用:
《MATLAB有限元分析与应用》是2004年4月清华大学出版社出版的图书,作者是卡坦,译者是韩来彬。
内容简介:
《MATLAB有限元分析与应用》特别强调对MATLAB的交互应用,书中的每个示例都以交互的方式求解,使读者很容易就能把MATLAB用于有限分析和应用。
另外,《MATLAB有限元分析与应用》还提供了大量免费资源。
《MATLAB有限元分析与应用》采用当今在工程和工程教育方面非常流行的数学软件MATLAB来进行有限元的分析和应用。
《MATLAB有限元分析与应用》由简单到复杂,循序渐进地介绍了各种有限元及其分析与应用方法。
书中提供了大量取自机械工程、土木工程、航空航天工程和材料科学的示例和习题,具有很高的工程应用价值。
《有限元基础教程》_【MATLAB算例】3.2.5(2)__四杆桁架结构的有限元分析(Bar2D2Node)
【MATLAB 算例】3.2.5(2) 四杆桁架结构的有限元分析(Bar2D2Node)如图3-8所示的结构,各个杆的弹性模量和横截面积都为4229.510/E N mm =⨯, 2100A mm =。
试基于MATLAB 平台求解该结构的节点位移、单元应力以及支反力。
图3-8 四杆桁架结构解答:对该问题进行有限元分析的过程如下。
(1) 结构的离散化与编号对该结构进行自然离散,节点编号和单元编号如图3-8所示,有关节点和单元的信息见表3-1~表3-3。
(2)计算各单元的刚度矩阵(基于国际标准单位)建立一个工作目录,将所编制的用于平面桁架单元分析的4个MA TLAB 函数放置于该工作目录中,分别以各自函数的名称给出文件名,即:Bar2D2Node_Stiffness ,Bar2D2Node_Assembly ,Bar2D2Node_Stress ,Bar2D2Node_Forces 。
然后启动MATLAB ,将工作目录设置到已建立的目录中,在MATLAB 环境中,输入弹性模量E 、横截面积A ,各点坐标x1,y1,x2,y2,x3,y3,x4,y4,角度alpha 1, alpha 2和alpha 3,然后分别针对单元1,2,3和4,调用4次Bar2D2Node_Stiffness ,就可以得到单元的刚度矩阵。
相关的计算流程如下。
>>E=2.95e11;>>A=0.0001;>>x1=0;>>y1=0;>>x2=0.4;>>y2=0;>>x3=0.4;>>y3=0.3;>>x4=0;>>y4=0.3;>> alpha1=0;>> alpha2=90;>> alpha3=atan(0.75)*180/pi;>> k1=Bar2D2Node_Stiffness (E,A,x1,y1,x2,y2,alpha1)k1 = 73750000 0 -73750000 00 0 0 0-73750000 0 73750000 00 0 0 0>> k2=Bar2D2Node_Stiffness (E,A,x2,y2,x3,y3,alpha2)k2 = 1.0e+007 *0.0000 0.0000 -0.0000 -0.00000.0000 9.8333 -0.0000 -9.8333-0.0000 -0.0000 0.0000 0.0000-0.0000 -9.8333 0.0000 9.8333>> k3=Bar2D2Node_Stiffness (E,A,x1,y1,x3,y3,alpha3)k3 = 1.0e+007 *3.7760 2.8320 -3.7760 -2.83202.8320 2.1240 -2.8320 -2.1240-3.7760 -2.8320 3.7760 2.8320-2.8320 -2.1240 2.8320 2.1240>> k4=Bar2D2Node_Stiffness (E,A,x4,y4,x3,y3,alpha1)k4 = 73750000 0 -73750000 00 0 0 0-73750000 0 73750000 00 0 0 0(3) 建立整体刚度方程由于该结构共有4个节点,因此,设置结构总的刚度矩阵为KK (8×8),先对KK 清零,然后四次调用函数Bar2D2Node _Assembly 进行刚度矩阵的组装。
有限元的matlab编程
本程序的网架位移求解函数附在主程序后面,主程序运行时调用该函数。
完整ppt
14
几何建模 定义荷载
用自定义输入
加下划线的为单元编号
完整ppt
5
形成等效荷载列阵
f=[0;0;0;a;0;0;0;a;0;0;0;a;0;0;0;a;0;0;0;a];%每个节点两个自由度,a 为之前输入的节点力
集成总刚:
获取单元两端节点坐标
xi = Node( Element( ie, 1 ), 1 ) ;%ie为单元号,以下相同 yi = Node( Element( ie, 1 ), 2 ) ; xj = Node( Element( ie, 2 ), 1 ) ; yj = Node( Element( ie, 2 ), 2 ) ;
节点号:',num2str(Element(ie,2)),' 轴力:',num2str(nodef(1))] ) ;
完整ppt
12
end
例二:网架
完整ppt
13
思路分析
几点说明
网架是由多根杆件按照一定的网格形式通过节点连结而成的空间结构。构成 网架的基本单元有三角锥,三棱体,正方体,截头四角锥等。鉴于网架的形 式较多,本程序提供一种通用的网架输入方法,但录入较为繁琐,同时提供 一种正放四角锥网架的简易输入方法作为典型。
有限元编程示例
完整ppt
1
例一:桁架
题目描述:
如下图所示的平面桁架,杆件长度、弹性模量、截 面积以及所受节点力P的大小可以自行定义。求节 点位移及杆件轴力。
matlab有限元分析实例
matlab有限元分析实例有限元是一种数值求解偏微分方程的方法。
基本过程大致是设置形函数,离散,形成求解矩阵,数值解矩阵,后处理之类的。
MATLAB要把这些过程均自己实现,不过在数值求解矩阵时可以调用已有函数。
可以理解为MATLAB是一个通用的计算器,当然它的功能远不止如此。
而ANSYS之类的叫做通用有限元软件,针对不同行业已经将上述过程封装,前后处理也比较漂亮,甚至不太了解有限元理论的人也能算些简单的东西,当然结果可靠性又另说了。
比较两者,ANSYS之类的用起来容易得多,但灵活性不如MATLAB。
MATLAB用起来很困难,也有人做了一些模块,但大多数只能解决一些相对简单的问题。
对于大多数工程问题,以及某些领域的物理问题,一般都用通用有限元软件,这些软件还能添加一些函数块,用以解决一些需要额外设置的东西。
但是对于非常特殊的问题,以及一般性方程的有限元解,那只能用MATLAB或C,Fortran之类的了。
matlab的优势是提供一个矩阵化的编程语言,matlab采用COO格式存储稀疏矩阵,在组装全局刚度阵的时候,可以由下标直接组装访问元素,处理单个单元的时候,对于二维模型,需要至少二重循环,加上对单元的遍历,至少循环层数较多。
众所周知,matlab在处理循环时速度较慢,特别是循环嵌套,速度更慢。
采用parfor,也不成熟。
优势就是在矩阵求解,雅可比矩阵与行列式求解的时候,有函数,又快又香。
对于C++,一般采用CSR格式存储稀疏阵,矩阵的行列,非零元,都要注意其存储位置。
编程求解比较费劲。
在处理单个单元的的时候行列式,小矩阵运算,都要自己编程,这个时候逻辑一定要清晰,代码相关的注释一定要写好。
否则编着编者,就把自己绕进去。
采用C++的福音是做并行时更容易控制,采用openmpi并行处理灵活。
在线性方程求解时,intel MKL是另外一个福音,intel的处理器还是intel理解,真的超级快,有机会你可以试一下,体会一下CPU 所有线程飙到100%的快感。
matlab有限元分析实例
为了使用matlab进行有限元分析,我们必须首先学习MATLAB 的基本操作以及使用MATLAB进行有限元分析的基本操作。
复习:最后一课分析弹簧系统,得出系统刚度矩阵matlab有限元分析20140226 matlab有限元分析20140226为了使用matlab进行有限元分析,我们首先必须学习MATLAB的基本操作,还要学习使用matlab进行有限元分析的基本操作。
1.回顾:最后一堂课分析了弹簧系统并推导出系统刚度矩阵。
2.为了使用matlab进行有限元分析,我们应该首先学习MATLAB的基本操作,然后再学习使用MATLAB进行有限元分析的基本操作。
Matlab有限元分析20140226为了使用Matlab进行有限元分析,首先必须学习MATLAB的基本操作以及使用MATLAB进行有限元分析的基本操作。
1.回顾:在上一课中对弹簧系统进行了分析,并得出了系统刚度矩阵。
2. MATLAB有限元分析的基本运算单位划分(选择哪个单位和划分多少个元素)分为两部分:单位选择和单位数回顾:最后一课分析了运算的基础MATLAB弹簧系统有限元分析matlab有限元分析20140226为了使用matlab进行有限元分析,我们首先应该学习MATLAB的基本操作,也要学会使用matlab进行有限元分析的基本操作。
1.回顾:最后一堂课分析了弹簧系统并推导出系统刚度矩阵。
2. Matlab导出系统刚度矩阵。
Matlab有限元分析的操作基础Matlab有限元分析20140226为了使用MATLAB进行有限元分析,首先必须学习MATLAB的基本操作,还必须学会使用Matlab 进行有限元分析的基本操作。
matlab 有限元 编程
MATLAB 是一种高级编程语言和交互式环境,常用于数值计算、数据分析、算法开发、数据可视化以及应用开发等。
在有限元分析(FEA)中,MATLAB 可以用来编写和运行有限元程序。
以下是一个简单的 MATLAB 有限元分析编程示例:这个例子假设你正在解决一个一维拉普拉斯方程,其形式为 -d^2u/dx^2 = f,在区间 [0, 1] 上。
```matlab% 参数定义L = 1; % 长度h = 0.01; % 步长N = L/h; % 元素数量x = linspace(0, L, N+1); % x 向量% 创建有限元网格nodes = x(1:N+1);elements = [nodes(2:end-1) nodes(1:end-1) nodes(2:end)]; % 定义载荷向量 ff = sin(pi*x);% 定义刚度矩阵 A 和载荷向量 FA = zeros(N, N);F = zeros(N, 1);for i = 1:Nfor j = 1:i-1A(i, j) = A(j, i) = h^2/(6*(x(i)-x(j))^2);endF(i) = -f(i)*h/2;end% 使用 MATLAB 的线性方程求解器u = A\F;% 绘制结果plot(x, u);xlabel('x');ylabel('u');title('有限元解');```这个代码首先定义了问题的参数,然后创建了一个有限元网格。
然后,它定义了刚度矩阵 A 和载荷向量 F。
最后,它使用 MATLAB 的线性方程求解器来求解方程,并绘制结果。
请注意,这只是一个非常简单的例子。
在实际的有限元分析中,问题可能会更复杂,并且需要更多的编程工作。
有限元的MATLAB解法
有限元的MATLAB解法1.打开MATLAB。
2.输入“pdetool”再回车,会跳出PDE Toolbox的窗口(PDE意为偏微分方程,是partial differential equations的缩写),需要的话可点击Options菜单下Grid命令,打开栅格。
3.完成平面几何模型:在PDE Toolbox的窗口中,点击工具栏下的矩形几何模型进行制作模型,可画矩形R,椭圆E,圆C,然后在Set formula栏进行编辑并(如双脊波导R1+R2+R3改为RI-R2-R3,设定a、b、s/a、d/b的值从而方便下步设定坐标)用算术运算符将图形对象名称连接起来,若还需要,可进行储存,形成M文件。
4.用左键双击矩形进行坐标设置:将大的矩形left和bottom都设为0,width是矩形波导的X轴的长度,height是矩形波导的y轴的长度,以大的矩形左下角点为原点坐标为参考设置其他矩形坐标。
5.进行边界设置:点击“Boundary”中的“Boundary Mode”,再点击“Boundary”中的“Specify Boundary Conditions”,选择符合的边界条件,Neumann为诺曼条件,Dirichlet为狄利克雷条件,边界颜色显示为红色。
6.进入PDE模式:点击"PDE"菜单下“PDE Mode”命令,进入PDE 模式,单击“PDE Specification”,设置方程类型,“Elliptic”为椭圆型,“Parabolic”为抛物型,“Hyperbolic”为双曲型,“Eigenmodes”为特征值问题。
7.对模型进行剖分:点击“Mesh”中“Initialize Mesh”进行初次剖分,若要剖的更细,再点击“Refine Mesh”进行网格加密。
8.进行计算:点击“Solve”中“Solve PDE”,解偏微分方程并显示图形解,u值即为Hz或者Ez。
9.单击“Plot”菜单下“Parameters”选项,打开“Plot Selection”对话框。
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
1.物理现象:这个对工程师来说是直观的物理现象和物理量,温
度多少度,载荷是多大等等。
通常来说,用户界面中呈现的、用户对工程问题进行设置时输入的都是此类信息。
2.数学方程:将物理现象翻译成相应的数学方程,例如流体对应
的是NS方程,传热对应的是传热方程等等;大部分描述这些现象的方程在空间上都是偏微分方程,偶尔也有ODE(如粒子轨迹、化学反应等)。
在这个层面,软件把物理现象“翻译”
为以解析式表示的数学模型。
3.数值模型:在定义了数学模型,并执行了网格剖分后,商业软
件会将数学模型离散化,利用有限元方法、边界元法、有限差分法、不连续伽辽金法等方法生成数值模型。
软件会组装并计算方程组雅可比矩阵,并利用求解器求解方程组。
这个层面的计算通常是隐藏在后台的,用户只能通过一些求解器的参数来干预求解。
有限元是一种数值求解偏微分方程的方法。
基本过程大致是设置形函数,离散,形成求解矩阵,数值解矩阵,后处理之类的。
MATLAB要把这些过程均自己实现,不过在数值求解矩阵时可以调用已有函数。
可以理解为MATLAB是一个通用的计算器,当然它的功能远不止如此。
而ANSYS之类的叫做通用有限元软件,针对不同行业已经将上述过程封装,前后处理也比较漂亮,甚至不太了解有限元理论的人也能算些简单的东西,当然结果可靠性又另说了。
比较两者,ANSYS之类的用起来容易得多,但灵活性不如MATLAB。
MATLAB用起来很困难,也有人做了一些模块,但大多数只能解决一些相对简单的问题。
对于大多数工程问题,以及某些领域的物理问题,一般都用通用有限元软件,这些软件还能添加一些函数块,用以解决一些需要额外设置的东西。
但是对于非常特殊的问题,以及一般性方程的有限元解,那只能用MATLAB或C,Fortran之类的了。