dynamic-Wilson-K.zip_wilson_混合动力

function [fequ,DisN,dD] = Fdis(Pos2,Dis1,Vel1,Acc1) %Fdis 利用刚度计算位移增量dis % 此处显示详细说明 global belta deltaT L Nnode d w EA ro ksoil theta Can=1.98;Cat=0.2;g=9.8;N=Nnode-1; Len0=L/N; for i=1:1:N Len(i)=sqrt((Pos2(1,i+1)-Pos2(1,i))^2+(Pos2(2,i+1)-Pos2(2,i))^2); Matm(i)=ro*Cat*pi*d^2/4*Len0; Manm(i)=ro*Can*pi*d^2/4*Len0; m(i)=w/g*Len0; phim(i)=atan((Pos2(2,i+1)-Pos2(2,i))/(Pos2(1,i+1)-Pos2(1,i))); end for i=2:N Mat(i)=0.5*(Matm(i)+Matm(i-1)); Man(i)=0.5*(Manm(i)+Manm(i-1)); phi(i)=0.5*(phim(i)+phim(i-1)); end %% 力的计算 Wt=[0.5*w*Len0,0.5*w*(Len0+Len0)*ones(1,N-1),0.5*w*Len0]; Fd=Fdrag_dyn(Len,phim,Vel1); z=Pos2(2,:); for i=2:N if z(i)>0 Fs(i)=0; else Fs(i)=0; % Fs(i)=0.5*ksoil*(Len(i-1)+Len(i))*z(i); end end %% 张力计算 for i=1:N Tx(i)=EA*(Pos2(1,i+1)-Pos2(1,i))/Len0-EA*(Pos2(1,i+1)-Pos2(1,i))/Len(i); Tz(i)=EA*(Pos2(2,i+1)-Pos2(2,i))/Len0-EA*(Pos2(2,i+1)-Pos2(2,i))/Len(i); end %% 惯性力计算 for i=2:N I1(i)=m(i)+Mat(i)*(sin(phi(i)))^2+Man(i)*(cos(phi(i)))^2; I3(i)=m(i)+Man(i)*(cos(phi(i)))^2+Mat(i)*(sin(phi(i)))^2; I2(i)=(Mat(i)-Man(i))*(sin(phi(i))*cos(phi(i))); Fix(i)=I1(i)*Acc1(1,i)+I2(i)*Acc1(2,i); Fiz(i)=I2(i)*Acc1(1,i)+I3(i)*Acc1(2,i); end %% 构建质量矩阵 Mnode=cell(1,N-1); for j=1:N-1 Mnode{1,j}=[I1(j+1),I2(j+1);I2(j+1),I3(j+1)]; end M=blkdiag(Mnode{:}); %% 构建函数 for i=2:N Fx(i)=Fd(1,i)+Tx(i)-Tx(i-1)-Fix(i); Fz(i)=Fd(2,i)+Tz(i)-Tz(i-1)-Wt(i)+Fs(i)-Fiz(i); end fequ=reshape([Fx(2:end);Fz(2:end)],2*(N-1),1); %% 构建刚度矩阵 Ksin1=cell(1,N);Ksin2=cell(1,N);Knode=cell(1,N-1); for i=1:N Ksin1{1,i}=[-EA/Len0+EA/Len(i)-EA*(Pos2(1,i+1)-Pos2(1,i))^2/(Len(i)^3),-EA*(Pos2(1,i+1)-Pos2(1,i))*(Pos2(2,i+1)-Pos2(2,i))/(Len(i)^3);... -EA*(Pos2(1,i+1)-Pos2(1,i))*(Pos2(2,i+1)-Pos2(2,i))/(Len(i)^3),-EA/Len0+EA/Len(i)-EA*(Pos2(2,i+1)-Pos2(2,i))^2/(Len(i)^3)]; Ksin2{1,i}=[EA/Len0-EA/Len(i)+EA*(Pos2(1,i+1)-Pos2(1,i))^2/(Len(i)^3),EA*(Pos2(1,i+1)-Pos2(1,i))*(Pos2(2,i+1)-Pos2(2,i))/(Len(i)^3);... EA*(Pos2(1,i+1)-Pos2(1,i))*(Pos2(2,i+1)-Pos2(2,i))/(Len(i)^3),EA/Len0-EA/Len(i)+EA*(Pos2(2,i+1)-Pos2(2,i))^2/(Len(i)^3)]; end for j=1:N-1 Knode{1,j}=Ksin1{1,j+1}+Ksin2{1,j}; end K=blkdiag(Knode{:}); %% 构建雅可比矩阵 % J=K-6*M/((theta*deltaT)^2); J=K-M/(belta*deltaT^2); dD=-inv(J)*(fequ); dDis=cat(2,[0;0],reshape(dD,2,[]),[0;0]); DisN=Dis1+dDis; end