Matlab 时程分析

0 动力平衡方程及相关参数取值

波浪、风载作用下的单桩动力反应计算（把结构简化为质点剪切型）

式中：u为结构水平位移；桩-结构集中质量矩阵；体系的阻尼矩阵；

体系的阻尼矩阵由结构和土体的阻尼矩阵集成，其中结构的阻尼按瑞雷阻尼理论，土体阻尼由材料阻尼和辐射阻尼组成。

体系的刚度矩阵；体系的刚度矩阵由结构和土体的刚度矩阵集成，土体刚度由动力P-Y曲线对Y求导得到。

地震、波浪、风载作用下的单桩动力反应计算

1 自由场地震分析（远离桩，取单位面积土柱）

土的刚度矩阵：； 

土的阻尼矩阵：；离桩较远，可采用刚度比例阻尼

2 桩土相互作用地震分析

也可写为

桩-结构集中质量矩阵；结构的阻尼矩阵；土体的阻尼矩阵；

结构的刚度矩阵；土体的刚度矩阵；基岩加速度。

自由场反应； 

1 结构动力响应数值方法

威尔逊法，隐式积分格式，一般取1.4时无条件稳定。=1即为常规线性加速度法。下面是威尔逊法的有限增量法，适用于非线性的动力求解。通过求解每一步积分的位移增量，计算速度、加速度的增量，最后求出每一步的位移、速度、加速度。

    非线性系统的威尔逊法  [摘自海洋工程结构动力学，唐友刚]

在一个延伸的计算步长内取：

          （>1.37）                      8-107 

在延伸的时间步长上，用标准的线性加速度法计算加速度增量，按此结果，用线性内插法求得在常规步长上的增量。当>1.37时，无条件稳定。

根据位移、速度、加速度的基本关系，可以得到时间步长的速度、位移：

                  8-108

其中，“”表明是与延伸的时间步长相应的增量。求解式8-108，得出表示的和，再代入运动方程：

式中：， 

，。

按延伸的时间步长写为：

                                           8-109a

                                   8-109b

        8-109c

求解拟静力关系式8-109a，得，并代入下式中，得到延伸时段内的加速度增量

                           8-110

由此，按线性内插法能得到常规步长的加速度增量。

                                           8-111

然后，按常规的步长由8-101，给出对应的速度和位移增量向量。

                        8-101

从而按照式8-105和8-106得到下一段的初始条件。

t+时刻的位移、速度向量、加速度。

                                     8-105

根据t+时刻的动力平衡条件求得，即

             8-106

式中：和分别代表从t+时刻的速度和位移（如材料特性与时程有关，则还包括过去的历史条件）求得的阻尼力和刚度力向量。

2 Matlab编程

***********************************************待参考

%威尔逊法主要程序

编程------------------------------------------------------------------------------------------------

符号定义

delt_Y---; delttao_Y---; disp---位移反应; velo---速度反应; acce---加速度反应;

time---计算时间;  delt_T---,计算时间间隔;  K---刚度矩阵;  C---阻尼矩阵; 

N---质点数;  theta---;  A(i)---加速度输入值

编程---------------------------------------------------------------------------------------------------

Function y=wilson(time,delt_T,K,C,N)

n=time/delt_T;

tao=theta*delt_T;

for t=delt_T: delt_T:time;

   i=t/ delt_T;

   for i=1:1:time/delt_T

F( : ,i)= M*A(i)*ones(N ,1);

A( :,i)=A(i)* ones(N ,1);

end

   if t= delt_T;

velo( : ,i)=zeros(N,1);        %结构速度初值为0

     disp( : ,i)= zeros(N,1);        %结构位移初值为0

acce( : ,i)=inv(M)*(F( : ,i)-K*disp( : ,i)-C*velo( : ,i));

   else

     tao_A( : ,i-1)=A(: ,i-1)+tao/delt_T*(A( : ,i)-A( : ,i-1));

     tao_K( : ,i-1)=K+M *6/tao^2 +C*3/tao;

     delttao_F( : ,i-1)=M*( tao_A( : ,i-1)- A(: ,i-1))+M*(6/tao*velo( : ,i-1)+3*acce( : ,i-1))+…

                 C*(3*velo( : ,i-1)+tao/2*acce( : ,i-1));

     delttao_disp( : ,i-1)=inv(tao_K( : ,i-1))* delttao_F( : ,i-1);

     delttao_acce( : ,i-1)= delttao_disp( : ,i-1)*6/tao^2 -velo( : ,i-1)*6/tao-acce( : ,i-1)*3;

     delt_acce( : ,i-1)=delttao_acce( : ,i-1)/theta; 

delt_velo( : ,i-1)=acce( : ,i-1)*delt_T+ delt_acce( : ,i-1)* delt_T/2;

delt_disp( : ,i-1)=velo( : ,i-1)* delt_T+ acce( : ,i-1)*(delt_T^2)/2+ delt_acce( : ,i-1)* (delt_T^2)/6;

disp( : ,i)=disp( : ,i-1)+delt_diso( : ,i-1);

     velo( : ,i)=velo( : ,i-1)+ delt_velo( : ,i-1);

     acce( : ,i)=inv(M)*(M*A( : ,i)-K*disp( : ,i)-C*velo( : ,i)

SD( : ,i)=disp( : ,i);

SV( : ,i)=velo( : ,i);

SA( : ,i)=acce( : ,i);

end

***********************************************

以下未编好

% 地震自由场分析---远离桩结构，土柱弹性

编程------------------------------------------------------------------------------------------------

符号定义

delt_Y_f ---; delttao_Y_f ---; disp_f ---土位移反应; velo_f ---土速度反应; acce_f ---土加速度反应;

time---计算时间;  delt_T---,计算时间间隔; Mtu—土质量;  Ktu---土刚度矩阵;  Ctu---土阻尼矩阵; 

N---质点数;  theta---;  A(i)---加速度输入值

h---土层厚度; rou---土层密度; G---土层剪切模量;w1,w2---自由场基频;Df---自由场阻尼比

编程---------------------------------------------------------------------------------------------------

Function y=EQ_FREE_wilson(time,delt_T ,N,rou,h,G,w1,w2,Df)

n=time/delt_T;

tao=theta*delt_T;

for t=delt_T: delt_T:time;

   i=t/ delt_T;

   for i=1:1:time/delt_T

F_f( : ,i)= M*A(i)*ones(N ,1);

A( :,i)=A(i)* ones(N ,1);

end

   if t= delt_T;

velo_f( : ,i)=zeros(N,1);        %结构速度初值为0

     disp_f( : ,i)= zeros(N,1);        %结构位移初值为0

acce_f( : ,i)=inv(Mtu)*(F( : ,i)-Ktu *disp( : ,i)-Ctu *velo( : ,i));

   else

     tao_A( : ,i-1)=A(: ,i-1)+tao/delt_T*(A( : ,i)-A( : ,i-1));

     tao_K_f( : ,i-1)=Ktu +Mtu *6/tao^2 +Ctu *3/tao;

     delttao_F_f( : ,i-1)=Mtu *( tao_A( : ,i-1)- A(: ,i-1))+Mtu *(6/tao*velo( : ,i-1)+3*acce( : ,i-1))+…

                 Ctu *(3*velo( : ,i-1)+tao/2*acce( : ,i-1));

     delttao_disp_f( : ,i-1)=inv(tao_K_f( : ,i-1))* delttao_F_f( : ,i-1);

     delttao_acce_f( : ,i-1)= delttao_disp_f( : ,i-1)*6/tao^2 -velo_f( : ,i-1)*6/tao-acce_f( : ,i-1)*3;

     delt_acce_f( : ,i-1)=delttao_acce_f( : ,i-1)/theta; 

delt_velo_f( : ,i-1)=acce_f( : ,i-1)*delt_T+ delt_acce_f( : ,i-1)* delt_T/2;

delt_disp_f( : ,i-1)=velo_f( : ,i-1)* delt_T+ acce_f( : ,i-1)*(delt_T^2)/2+ delt_acce_f( : ,i-1)* (delt_T^2)/6;

disp_f( : ,i)=disp_f( : ,i-1)+delt_diso_f( : ,i-1);

     velo_f( : ,i)=velo_f( : ,i-1)+ delt_velo_f( : ,i-1);

     acce_f( : ,i)=inv(Mtu)*(Mtu*A( : ,i)-Ktu*disp( : ,i)-Ctu*velo( : ,i)

D_f( : ,i)=disp_f( : ,i);

V_f( : ,i)=velo_f( : ,i);

A_f( : ,i)=acce_f( : ,i);

end

function [Mtu Ktu Ctu]=tu(rou(i),h(i),G(i),w1,w2,Df)

for i=1:1:N

   h(0)=0;h(N+1)=0;

   rou(0)=0;rou(N+1)=0;

mtu(i)=0.5*h(i-1)*rou(i-1)+ 0.5*h(i)*rou(i);

   ktu(i)=G(i)/h(i);

end

Mtu=diag(mtu);

ktu0=[0 ktu(1:n-1)];

Ktu= diag(ktu)+ diag(ktu0)- diag(ktu0,1)- diag(ktu0,-1);

Ctu=2/(w1+w2)*Df*Ktu;

% 桩土地震分析---结构弹性，土反力-P-Y，质点取结构质量

编程------------------------------------------------------------------------------------------------

符号定义

delt_Y ---; delttao_Y ---; disp_f ---土位移反应; velo_f ---土速度反应; acce_f ---土加速度反应;

time---计算时间;  delt_T---,计算时间间隔; Mtu—土质量;  Ktu---土刚度矩阵;  Ctu---土阻尼矩阵; 

N---质点数;  theta---;  A(i)---加速度输入值

h---土层厚度(下部桩段长度);  rou---土层密度;  G---土层剪切模量; w1,w2---自由场基频; Df---自由场阻尼比

disp---位移反应;  velo---速度反应;  acce---加速度反应;

Mp---桩柱质量矩阵;  Kp---桩柱刚度矩阵;  Cp---桩柱阻尼矩阵; NN---上部桩质点数; hh---上部桩柱每段长度, R---桩柱线密度; EI---桩柱刚度

Function y=EQ _wilson(time,delt_T ,NN,R,hh,G,w1,w2,Df,L,R，EI)

n=time/delt_T;

tao=theta*delt_T;

for t=delt_T: delt_T:time;

   i=t/ delt_T;

   for i=1:1:time/delt_T

F ( : ,i)= M*A(i)*ones(N ,1);

A( :,i)=A(i)* ones(N ,1);

end

% 下面是计算土的刚度，阻尼

for i=1:1:N

Cmat( : ,i)=2*Dtumax* (Kpy*rou(i))^0.5*(1-Kpy( : ,i)/Kh)

h(0)=0;h(N+1)=0;

rou(0)=0;rou(N+1)=0;

mtu(i)=0.5*h(i-1)*rou(i-1)+ 0.5*h(i)*rou(i);

   ktu(i)=G(i)/h(i);

end

Mtu=diag(mtu);

ktu0=[0 ktu(1:n-1)];

Ktu= diag(ktu)+ diag(ktu0)- diag(ktu0,1)- diag(ktu0,-1);

Ctu=2/(w1+w2)*Df*Ktu;

% 上面是计算土的刚度阻尼

% 下部桩柱的质量，刚度

for i=1:1:N

   h(0)=0;h(N+1)=0;

   R(0)=0;R(N+1)=0;

m2 (i)=0.5*h(i-1)*R(i-1)+ 0.5*h(i)*R(i);

   k2(i)=12*EI/h(i)^3;

end

%  上部桩柱的质量，刚度

for i=1:1:NN

   hh(0)=0;hh(NN+1)=0;

   R(0)=0;R(NN+1)=0;

m1 (i)=0.5*hh(i-1)*R(i-1)+ 0.5*hh(i)*R(i);

   k1(i)=12*EI/hh(i)^3;

end

Mp=[m1 m2]

M=diag(Mp);

ktu0=[0 ktu(1:n-1)];

Ktu= diag(ktu)+ diag(ktu0)- diag(ktu0,1)- diag(ktu0,-1);

Ctu=2/(w1+w2)*Df*Ktu;

   if t= delt_T;

velo_f( : ,i)=zeros(N,1);        %结构速度初值为0

     disp_f( : ,i)= zeros(N,1);        %结构位移初值为0

acce_f( : ,i)=inv(M)*(F( : ,i)-K *disp( : ,i)-C *velo( : ,i));

   else

     tao_A( : ,i-1)=A(: ,i-1)+tao/delt_T*(A( : ,i)-A( : ,i-1));

     tao_K ( : ,i-1)=K +M *6/tao^2 +C *3/tao;

     delttao_F ( : ,i-1)=M *( tao_A( : ,i-1)- A(: ,i-1))+M *(6/tao*velo( : ,i-1)+3*acce( : ,i-1))+…

                 C *(3*velo( : ,i-1)+tao/2*acce( : ,i-1));

     delttao_disp ( : ,i-1)=inv(tao_K ( : ,i-1))* delttao_F ( : ,i-1);

     delttao_acce ( : ,i-1)= delttao_disp ( : ,i-1)*6/tao^2 -velo ( : ,i-1)*6/tao-acce ( : ,i-1)*3;

     delt_acce ( : ,i-1)=delttao_acce ( : ,i-1)/theta; 

delt_velo ( : ,i-1)=acce ( : ,i-1)*delt_T+ delt_acce ( : ,i-1)* delt_T/2;

delt_disp ( : ,i-1)=velo ( : ,i-1)* delt_T+ acce ( : ,i-1)*(delt_T^2)/2+ delt_acce ( : ,i-1)* (delt_T^2)/6;

disp( : ,i)=disp( : ,i-1)+delt_diso( : ,i-1);

     velo( : ,i)=velo( : ,i-1)+ delt_velo( : ,i-1);

     acce( : ,i)=inv(M)*(M*A( : ,i)-K*disp( : ,i)-C*velo( : ,i)

D_P( : ,i)=disp( : ,i);

V_P( : ,i)=velo( : ,i);

A_P( : ,i)=acce( : ,i);

end