东南大学 计算力学 中期作业
LSTR， 1， 2 LSTR， 2， 3 CM，_Y，LINE LSEL，，，， 1 CM，_Y1，LINE CMSEL，S，_Y CMSEL，S，_Y1 LATT，1，，1，，，，1 CMSEL，S，_Y CMDELE，_Y CMDELE，_Y1 FLST，5，1，4，ORDE，1 FITEM，5，1 CM，_Y，LINE LSEL，，，，P51X CM，_Y1，LINE CMSEL，，_Y LESIZE，_Y1，，，20，，，，，1 LMESH， 1 CM，_Y，LINE LSEL，，，， 2 CM，_Y1，LINE CMSEL，S，_Y CMSEL，S，_Y1 LATT，1，，1，，，，2 CMSEL，S，_Y CMDELE，_Y CMDELE，_Y1 FLST，5，1，4，ORDE，1 FITEM，5，2 CM，_Y，LINE LSEL，，，，P51X CM，_Y1，LINE CMSEL，，_Y LESIZE，_Y1，，，20，，，，，1 LMESH， 2 EPLOT FLST，2，1，1，ORDE，1 FITEM，2，1 D，P51X，，，，，，ALL，，，，， FLST，2，1，1，ORDE，1 FITEM，2，22 D，P51X，，，，，，UY，，，，， FLST，2，40，2，ORDE，2
6
东南大学 计算力学 中期作业
FITEM，2，1 FITEM，2，-40 SFBEAM，P51X，2，PRES，100000，，，，，，0 FINISH /SOL /STATUS，SOLU SOLVE FINISH /POST1 PRNSOL，U，Y !提取整体刚度矩阵和荷载向量 !FINISH !/SOL !/STATUS，SOLU !antype，7 !seopt，matname，1 !nsel，all !m，all，all !solve !selist，matname，3 FINISH /SOL /STATUS，SOLU SOLVE FINISH /POST1 PRNSOL，U，COMP !ETABLE，FXI，SMISC， 1 !ETABLE，FXJ，SMISC， 14 !ETABLE，SFYI，SMISC， 6 !ETABLE，SFYJ，SMISC， 19 !ETABLE，MZI，SMISC， 3 !ETABLE，MZJ，SMISC， 16 !PRETAB，FXI，FXJ，SFYI，SFYJ，MZI，MZJ Matlab程序： %总控 clear; Lk=matrix(); for i=1:length(Lk(:，1)) Disp=main( 1，Lk(i，:) );%1表示选用的为经典梁单元，0表示timoshenko梁 mid_deflec(i，1)=Disp(5);
function [ Lk ] = matrix( ) %%矩阵文件，单元的长度，以及剪切系数 Lk=[0.2 1.2 1.2 0.3 1.2 1.2 0.5 1.2 1.2 7
东南大学 计算力学 中期作业
Disp=main( 0，Lk(i，:) ); mid_deflec(i，2)=Disp(62); end fprintf('%s\\n\\t%s\\t\\t%s\\n'，'中点挠度计算结果:'，'经典梁元'，'Timoshenko梁元') for i=1:length(Lk(:，1)) 0.7 1.2 1.2 1.0 1.2 1.2 1.2 1.2 1.2 1.5 1.2 1.2 1.8 1.2 1.2 2.0 1.2 1.2 fprintf('%e\\t%e\\n'，mid_deflec(i，1)，mid_deflec(i，2)); 2.2 1.2 1.2 end 5.0 1.2 1.2 10 1.2 1.2 15 1.2 1.2 ]; end function [ Disp ] = main( flag，Lk ) %UNTITLED Summary of this function goes here % Detailed explanation goes here L=Lk(1); k1=Lk(2);k2=Lk(3); %经典梁里剪切系数取为1.2，能量法手算用的剪切系数取为1.2，Timoshenko梁这里也取1.2，而ansys是根据一定规则算出来个k %，不是手输，为了简便这里就认为取为1.2 h=0.4;b=0.2;E=2.1e11;u=0.3; if flag==1 %经典梁元，只需划分两个单元 n=1; else n=20; %Timoshenko梁元，必须把细化，否则结果不正确，总共化为2n段 end Joint=struct('NJoint'，[]，'NDOF'，[]，'Coord'，[]，'DOF'，[]，'Load'，[]); for i=1:(2*n+1) Joint.Coord(i，:)=[2*L/40*(i-1)，0]; end Joint.NJoint=2*n+1; Joint.DOF(:，1)=[1;1;1];Joint.DOF(:，2*n+1)=[0;1;0];%1表示受约束，另外每个结点位移都得输入，这里2结点正好能够默认为0. Elem=struct('NElem'，[]，'Def'，[]，'Load'，[]); for i=1:n Elem.Def(i，:)=[i i i+1 1.45*h*b， L/n， E， b*1.45^3*h^3/12，u，E/2/(1+u)，k1];%剪切矫正因子可能1.2617 Elem.Def(i+n，:)=[i+n n+i n+i+1 1.15*h*b， L/n， E， b*1.15^3*h^3/12，u，E/2/(1+u)，k2]; end Elem.NElem=2*n; for i=1:2*n
8
东南大学 计算力学 中期作业
Elem.Load(i，:)=[i 100000];%压力为正 end StructStiff = struct_Stiff( Joint，Elem，flag ); TotalLoad = PF_Load( Joint，Elem，flag ); Disp=Node_Disp(Joint，StructStiff，TotalLoad ); %求杆端力 timoshenko梁对不上啊 %for i=1:Elem.NElem; % InFL(:，i) = PF_EndInFl( i，Joint，Elem， Disp ，flag); %end end function [ StructStiff ] = struct_Stiff( Joint，Elem，flag ) %计算单元刚度矩阵，集成结构刚度矩阵 % Joint=结点信息，Elem=单元信息 % StructStiff=结构刚度矩阵StructStiff=zeros(3*Joint.NJoint); for i=1:Elem.NElem %计算单元刚度矩阵ElemStiff=Beam_Elem_Stiff(Elem.Def(i，:)，flag); %集成单元刚度矩阵 StructStiff=AssemElem(Elem.Def(i，:)，ElemStiff，StructStiff); end end function [ k ] = Beam_Elem_Stiff(ElemDef，flag ) %Struss_Elem_StiffGlobal函数是计算一个单元在结构坐标系（不是局部坐标系）中的单元刚度 %L=单元长度，Orient=位向A=ElemDef(4);L=ElemDef(5);E=ElemDef(6);I=ElemDef(7);u=ElemDef(8);G=ElemDef(9);k=ElemDef(10); %经典梁单元 if flag==1 As=A/k; b=12*E*I/G/As/L^2; k=[A*E/L 0 0 -A*E/L 0 0; 0 12*E*I/(1+b)/L^3 6*E*I/L^2/(1+b) 0 -12*E*I/L^3/(1+b) 6*E*I/L^2/(1+b); 0 6*E*I/L^2/(1+b) E*I*(4+b)/L/(1+b) 0 -6*E*I/L^2/(1+b) E*I*(2-b)/L/(1+b); -A*E/L 0 0 A*E/L 0 0; 0 -12*E*I/L^3/(1+b) -6*E*I/L^2/(1+b) 0 12*E*I/L^3/(1+b) -6*E*I/L^2/(1+b); 0 6*E*I/L^2/(1+b) E*I*(2-b)/L/(1+b) 0 -6*E*I/L^2/(1+b)
9
东南大学 计算力学 中期作业
E*I*(4+b)/L/(1+b)]; %timoshenko梁单元，自由度采用的是每个结点3个，而ansys里beam189每个自由度数为6个。 else dd=G*A/k/L;gg=E*I/L; %该单元刚度矩阵采用的是减缩积分后的，即修正后的，这也是ansys采用的单刚形式 k=[A*E/L 0 0 -A*E/L 0 0; 0 dd dd*L/2 0 -dd dd*L/2; 0 dd*L/2 dd*L^2/4+gg 0 dd*(-L/2) dd*L^2/4-gg; -A*E/L 0 0 A*E/L 0 0; 0 -dd dd*(-L/2) 0 dd dd*(-L/2); 0 dd*L/2 dd*L^2/4-gg 0 dd*(-L/2) dd*L^2/4+gg]; end end function [ B ] = AssemElem( iDef，A，B ) %AssemElem集成结构刚度矩阵和等效结点荷载列阵 %A单元刚度矩阵 %B结构刚度矩阵LocVec=[iDef(2)*3-2，iDef(2)*3-1，iDef(2)*3，iDef(3)*3-2，iDef(3)*3-1，iDef(3)*3]; ij=find(LocVec~=0); IJ=LocVec(LocVec~=0); if size(A，2)~=1 %集成单元刚度矩阵 B(IJ，IJ)=B(IJ，IJ)+A(ij，ij); else%集成单元等效结点荷载 B(IJ)=B(IJ)+A(ij); end end function [ Node_Disp ] = Node_Disp( Joint，StructStiff，TotalLoad ) %计算结点位移 v=find(Joint.DOF==1); v2=find(Joint.DOF==0); StructStiff(v，:)=[]; StructStiff(:，v)=[]; TotalLoad(v，:)=[]; Disp=StructStiff\\TotalLoad; Node_Disp(v2)=Disp; Node_Disp=Node_Disp'; end function [ InFL ] = PF_EndInFl( i，Joint，Elem， Disp，flag ) %计算杆端内力 % %计算第i个单元的杆端内力，Displ为此单元两个结点的位移Displ=Disp([3*Elem.Def(i，
10