有限元编程讲义.ppt

%定义转化矩阵 t=sqrt(l^2+n^2); if t==0 r=[0 m 0;-m 0 0;0 0 1]; else r=[l m n;-l*m/t t -m*n/t;-n/t 0 l/t]; end O=zeros(3,3); T=[r O;O r]; %整体坐标下的单刚矩阵 ke=E*A/L(i,1)*[1 0 0 -1 0 0;0 0 0 0 0 0 ;0 0 0 0 0 0 ;-1 0 0 1 0 0;0 0 0 0 0 0;0 0 0 0 0 0]; k=T*ke*T; %单刚矩阵的扩充，使之行数、列数与总刚对应 G=zeros(6,3*Nsum); I=[1 0 0;0 1 0;0 0 1]; G(1:3,3*p-2:3*p)=I; G(4:6,3*q-2:3*q)=I; K=G*k*G; Kz=Kz+K; %形成总刚矩阵 end %节点力矩阵的扩充 Rz=zeros(3*Nsum,1); for i=1:Psum j=P(i,1); Rz(3*j,1)=P(i,2); end %对角元素置1法处理刚度矩阵 与节点力矩阵 for i=1:3*Nsum if BCz(i)==0 Kz(i,:)=0; Kz(:,i)=0; Kz(i,i)=1; Rz(i)=0; end end u=Kz\Rz; %节点位移 L1=L; %单元长度矩阵 end 求解网架位移 具体算例 选择网架类型，0代表自由定义网架，1代表四角锥网架：1 输入网架上层节点行数：4 输入网架上层节点列数：6 输入网架上层节点列间距：2 输入网架上层节点行间距：2 输入网架上下层间距：2 请输入网架上层单元的截面面积：0.0003 请输入网架下层单元的截面面积：0.0003 请输入网架腹杆单元的截面面积：0.0004 定义节点荷载，若网架上层节点力与下层节点力均布，则输入0，否则输入1：0 请输入网架上层节点荷载：10000 请输入网架下层节点荷载：0 定义边界约束，若网架上层周边节点全约束，则输入0，若下层周边节点全约束，输入1，否则输入2：0 截面输入部分 * 有限元编程示例 题目描述: 如下图所示的平面桁架，杆件长度、弹性模量、截面积以及所受节点力P的大小可以自行定义。求节点位移及杆件轴力。 例一：桁架 解题思路： 建立模型 集成总刚 求解位移 求解杆件轴力 输出结果 建立模型： 定义节点坐标 Node = zeros(10,2) ; x=-1*L; %L为横杆长度 for i=1:2:10 x=x+L; Node(i,:)=[x 0]; end x=-1*L; for i=2:2:10 x=x+L; Node(i,:)=[x H];%H为竖杆长度 end 模型相关参数输入 H=input(竖杆长度(m):); L=input(水平杆长度(m):); E=input(杆件弹性模量(Gpa):); A=input(杆件截面积(m^2):); a=input(节点力P(kN)：); 节点编号方式 定义单元，即储存单元两端的节点号 Element=zeros(21,2); for i=1:2:7 Element(5/2*i-3/2,:)=[i,i+1]; Element(5/2*i-1/2,:)=[i,i+2]; Element(5/2*i+1/2,:)=[i,i+3]; end for i=2:2:8 Element(5*i/2-1,:)=[i,i+1]; Element(5*i/2,:)=[i,i+2]; end Element(21,:)=[9,10]; 加下划线的为单元编号 集成总刚： 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 ) ; 获取单元两端节点坐标 L = ( (xj-xi)^2 + (yj-yi)^2 )^(1/2) ; 计算杆件长度 形成等效荷载列阵 f=[0;0;0;a;0;0;0;a;0;0;0;a;0;0;0;a;0;0;0;a];%每个节点两个自由度，a为之前输入的节点力 计算从局部坐标到整体坐标的坐标转换矩阵T function T = TransformM