Matlab求解Timoshenko梁
第二帖求解Timoshenko梁
形函数：输入高斯积分点坐标，输出形函数和形函数对自然坐标的偏导数。function =shapetimoshenko(xi) %------------------------------------------------------------------------------if nargin ~= 1 error('Input error!');end %------------------------------------------------------------------------------syms XI;N1 = 1/2*(1-XI);N2 = 1/2*(1+XI); shapef = ;pdnf = diff(shapef,XI); %------------------------------------------------------------------------------shape =subs(shapef,XI,xi);pdn =subs(pdnf,XI,xi); end 雅可比矩阵：function = jacobimatrix(pdn, node) jm = pdn * node;invjm = inv(jm);pdc = invjm * pdn; end 应变矩阵：function B =strainmatrixtb(pdc, shape, type) if nargin == 2 B = pdc(1,1), shape(1,1), pdc(1,2),shape(1,2)];elseif nargin == 3 switch upper(type) case 'BENDING' B = 0, 0, 0, 0]; case 'SHEAR' B = pdc(1,1), shape(1,1),pdc(1,2), shape(1,2)]; otherwise error('Input error!') endelse error('Input error!');end end 单元刚度矩阵：function ke =elemstifftb(E, Nu, kapa, A, I, node) G = E/2/(1+Nu); D = ; =shapetimoshenko(0.577350269189626); =jacobimatrix(pdn1, node);Bb1 =strainmatrixtb(pdc1, shape1, 'bending');keb1 =Bb1'*D*Bb1*det(jm1); =shapetimoshenko(-0.577350269189626); =jacobimatrix(pdn2, node);Bb2 =strainmatrixtb(pdc2, shape2, 'bending');keb2 = Bb2'*D*Bb2*det(jm2); =shapetimoshenko(0); =jacobimatrix(pdn3, node);Bs3 =strainmatrixtb(pdc3, shape3, 'shear');kes =Bs3'*D*Bs3*det(jm3)*2; ke = keb1 + keb2 + kes; end 集成总刚度矩阵：function Fe =assembleforcetb(Fe, fe, i, j) Fe() =Fe() + fe; end
验证：E=10e7; Nu = 0.30;L= 1;thickness=0.001;I=thickness^3/12;A=thickness*1;EI=E*I;kapa=5/6;P=-1; nelem = 100;nnode = nelem+1;n=linspace(0,L,nnode);for i=1:nelem e(i,1)=i; e(i,2)=i+1;end K=zeros(nnode*2);F=zeros(nnode*2,1); for ei=1:nelem node =;ke = elemstifftb(E, Nu,kapa, A, I, node);fe =elemforcetb(P,node);K=assemblestifftb(K,ke,e(ei,1),e(ei,2));F=assembleforcetb(F,fe,e(ei,1),e(ei,2));end fixedNodes=find(abs(n(:,1))<1e-5);Dof=;activeDof=setdiff(',Dof);U=K(activeDof,activeDof)\F(activeDof);disp=zeros(nnode*2,1);disp(activeDof)=U;min(disp)
在求解Timoshenko梁的时候要注意剪力锁死问题。如验证中所采用的悬臂梁，若均采用两点积分，梁端位移仅为 -0.4539。若将剪力对刚度的贡献采用单点积分，则可以消除剪力锁死问题，梁端位移为-15.
fe =elemforcetb(P,node);缺少elemforcetb的function
页:
[1]