第三章_MATLAB有限元分析与应用.ppt

3、实例计算分析应用 E=210E6; A=0.003; L=2; k1= QuadraticBarElementStiffness(E,A,L); k2= QuadraticBarElementStiffness(E,A,L); n=5; K = zeros(n,n); K =QuadraticBarAssemble(K,k1,1,3,2); K =QuadraticBarAssemble(K,k2,3,5,4); U=zeros(4,1); F=[5;-10;-7;10]; KK=K; K=K(2:n,2:n); U=K\F; U=[0;U]; F=KK*U; u1=[U(1);U(3);U(2)]; f1= QuadraticBarElementForces(k1,u1); sigma1=QuadraticBarElementStresses(k1, u1, A); u2=[U(3);U(5);U(4)]; f2=QuadraticBarElementForces(k2,u2); sigma2=QuadraticBarElementStresses(k2, u2, A); §3-3 二次杆元 1、基本方程 §3-4 平面桁架元 平面桁架元是既有局部坐标又有总体坐标二维有限元，用线性函数描述 每个平面桁架元有二个节点（node） 单刚矩阵为： 总刚矩阵： 结构方程： 单元节点力： 2、MATLAB函数编写 %PlaneTrussElementLength This function returns the length of the % plane truss element whose first node has % coordinates (x1,y1) and second node has % coordinates (x2,y2). 2.1 计算单元长度 y = sqrt((x2-x1)*(x2-x1) + (y2-y1)*(y2-y1)); function y = PlaneTrussElementLength(x1,y1,x2,y2) §3-4 平面桁架元 2、MATLAB函数编写 %PlaneTrussElementStiffness This function returns the element % stiffness matrix for a plane truss % element with modulus of elasticity E, % cross-sectional area A, length L, and % angle theta (in degrees). % The size of the element stiffness % matrix is 4 x 4. 2.2 单元刚度矩阵的形成 x = theta*pi/180; C = cos(x); S = sin(x); y = E*A/L*[C*C C*S -C*C -C*S ; C*S S*S -C*S -S*S ; -C*C -C*S C*C C*S ; -C*S -S*S C*S S*S]; function y = PlaneTrussElementStiffness(E,A,L, theta) §3-4 平面桁架元 2、MATLAB函数编写 %PlaneTrussAssemble This function assembles the element stiffness % matrix k of the plane truss element with nodes % i and j into the global stiffness matrix K. % This function returns the global stiffness % matrix K aft