基于双折线本构模型的能量计算

1双折线本构模型

2耗能计算的方法

3单自由度体系的耗能计算

1双折线本构模型

双折线本构模型是在弹性本构模型的基础上引入屈服强度和屈服后的第二刚度对弹性本构模型进行修正。本构模型如下：

单元状态与位移的历程时刻相关，本构关系是确定单元或结构所处状态的关键，需要以本构关系为基础来编写确定当前关系状态的函数。

以后文使用的随动强化非线性单轴二折线本构为例，简单说明非弹性单元状态确定的思路：

（1）获取本构基本信息，如：初始刚度k0，屈服强度Fy，材料强化准则，强化系数等，构建初始本构曲线。

（2）根据上一时刻的已知条件F1，D1,以及当前时刻的位移D2，计算出抗力下限值Fa，抗力上限值Fb（屈服面的范围），和按照初始刚度不变的原则计算出当前位移D2的抗力值Fc，最终计算出正确抗力值F2。计算公式如下：

Fa=-Fy1-b*αn+bk0D2

Fb=Fy1-b*αp+bk0D2

Fc=F1+k0(D2-D1)

F2=max[Fa,min(Fb,Fc)]

未考虑受滞回影响的正负向强度调整系数αn，αp，其值取1。

（3）由于强化后的屈服强度与初始屈服强度有所区别，实际上应根据记录的历史抗力峰值，确定当前时刻对应的正负向屈服强度，每个时刻重新计算正负向强度调整系数αn，αp，重组本构曲线。

2耗能计算的方法

弹性能量=动能+弹性应变能；

塑性能量=阻尼耗能+非弹性应变耗能

总能量=弹性能量+塑性能量

对阻尼进行考虑（常阻尼与变阻尼的对比）：

实际情况中，结构在动力荷载作用下进入弹塑性状态，刚度与自振频率会发生变化，在此基础上，阻尼系数c也会随其变化；常阻尼以初始刚度，质量，阻尼比计算，不随时刻变化而变化；变阻尼计算理论与常阻尼一致，但刚度随时刻的变化而变化，所以阻尼亦随之变化。通过研究表明：分别采用常阻尼与变阻尼的结构体系计算得到的结果相差不大，为了简化计算，提高效率，均采用常阻尼。

3单自由度体系的耗能计算

选用NGA-68地震动记录，其加速度时程如下：

选定的单自由度体系模型，本构关系，阻尼比默认为0.05%，阻尼采用常阻尼，采用结果如下：

单自由度体系模型（质量m=1000kg）：

常阻尼：（c）

%计算结构的自振频率和阻尼系数

wn=sqrt(k0/m);

c=2*kesi*m*wn;

本构关系(屈服力为fy=1e6N,屈服强度k0=2e5，屈服第二刚度系数为0.2)；

计算结果如下：

弹性能量的时程分布：

动能Ek时程	弹性应变能Es时程

塑性能量的时程分布：

阻尼耗能Ed时程	非弹性应变耗能Ep时程

能量输入封闭体系中，总能量等于弹性时程的能量与塑性耗能的和，由上述时程图得出：动能以及弹性应变能在时程末为0，则表明输入的总能量都在时程末转化为塑性耗能。

下图为总塑性耗能：

本构关系(屈服力为fy=1e6N,屈服强度k0=2e5，屈服第二刚度系数为0)；

计算结果如下：

弹性能量的时程分布

动能Ek时程	弹性应变能Es时程

塑性能量的时程分布：

阻尼耗能Ed时程	非弹性应变耗能Ep时程

下图为总塑性耗能：

对比这两个不同本构关系的各个应变耗能，不同主要在弹性耗能时程变化方面，但总体的塑性耗能差别不大,约为5%左右。

参考内容：

1采用matlab简单编程计算，计算程序如下：

%%本构模型采用强化系数为0.2的双折线本构，模型为单自由度体系。
%%弹塑性单自由度体系的能量分析；
clear,clc
format shortEng
%输入结构参数（阻尼比，质量）
kesi=0.05;
m=1000;
%导入地震动记录，设定时间步数和步长
data=load('NGA-68.txt');
ug=9800*data';
n1=length(ug);
dt=0.01
%定义非线性单元本构；
global fy k0 b
fy=1000000;%屈服力
k0=2e5;%屈服强度
b=0.2;%折线段第二刚度系数
%计算结构的自振频率和阻尼系数
wn=sqrt(k0/m);
c=2*kesi*m*wn;
%指定控制参数γ和β的值，定义积分常数
gama=0.5;beta=1/4;
AA1=m/(beta*dt)+(gama/beta)*c;
AA2=m/(2*beta)+dt*(gama/2/beta-1)*c;
%求解结构反应
%定义初始条件
n2=6000;  %分析步数
Time(1,1)=0;%定义时间向量Time,初始值为0；
u(1,1)=0;%定义位移向量u，初始值为0；
v(1,1)=0;%定义速度向量v，初始值为0；
Rs(1,1)=0;%定义恢复力向量Rs，初始值为0；
k(1,1)=k0;%定义刚度向量k，初始值为k0；
P(1,1)=-m*ug(1,1);%定义外荷载向量P
aa(1,1)=(P(1,1)-c*v(1,1)-Rs(1,1))/m;%定义加速度向量aa
Ek(1,1)=0.5*m*v(1,1)^2;%动能
Es(1,1)=0.5*Rs(1,1)^2/k0;%势能
Ed(1,1)=0;%阻尼耗能
Er(1,1)=0;%总应变能
Ep(1,1)=0;%塑性耗能
%逐步积分
for i=2:n2
Time(1,i)=(i-1)*dt;
%求等效刚度kk
kk(1,i-1)=k(1,i-1)+(1/beta/dt/dt)*m+(gama/beta/dt)*c;
%求等效荷载增量dltPP
if i<=n1
P(1,i)=-m*ug(1,i);%激励中
else
ug(1,i)=0;%激励停止后自由振动，外荷载为0；
P(1,i)=0;
end
dltP(1,i-1)=P(1,i)-P(1,i-1);
dltPP(1,i-1)=dltP(1,i-1)+AA1*v(1,i-1)+AA2*aa(1,i-1);
%开始积分步内迭代
%初始化积分数据
utrial0=u(1,i-1);%初始迭代位移
Rtrial0=Rs(1,i-1);%初始迭代恢复力
dltR0=dltPP(1,i-1);%初始迭代残余力
KKt0=kk(1,i-1);%初始迭代等效切线刚度
dltU=1e-8;%变量
dltu=1e-8;%变量
while(dltu/dltU)>1e-4
dltu=dltR0/KKt0;
utrial=utrial0+dltu;
[Rtrial,Ktrial]=stateDetermineNSDOF(Rtrial0,utrial0,utrial);
KKt0=Ktrial+(gama/beta/dt)*c+(1/beta/dt/dt)*m;
dltF=Rtrial-Rtrial0+((gama/beta/dt)*c+(1/beta/dt/dt)*m)*dltu;

dltR0=dltR0-dltF;
dltU=dltU+dltu;
utrial0=utrial;
Rtrial0=Rtrial;
end
u(1,i)=utrial;
Rs(1,i)=Rtrial;
k(1,i)=Ktrial;
dltv=(gama/beta/dt)*(u(1,i)-u(1,i-1))-(gama/beta)*v(1,i-1)+(1-gama/2/beta)*dt*aa(1,i-1);
dltaa=(1/beta/dt/dt)*(u(1,i)-u(1,i-1))-(1/beta/dt)*v(1,i-1)-(1/2/beta)*aa(1,i-1);
v(1,i)=v(1,i-1)+dltv;
aa(1,i)=aa(1,i-1)+dltaa;

Ek(1,i)=0.5*m*v(1,i)^2;
Es(1,i)=0.5*Rs(1,i)^2/k0;
Ed(1,i)=Ed(1,i-1)+c*0.5*(v(1,i)+v(1,i-1))*(u(1,i)-u(1,i-1));
Er(1,i)=Er(1,i-1)+0.5*(Rs(1,i)+Rs(1,i-1))*(u(1,i)-u(1,i-1));
Ep(1,i)=Er(1,i)-Es(1,i);
end
%绘制动能Ek时程
figure(1)
plot(Time,Ek)
xlabel('Time(s)')
ylabel('Ek')
title('Ek vs Time')
saveas(gcf,'Ek vs Time,png')
%绘制弹性应变能Es时程
figure(2)
plot(Time,Es)
xlabel('Time(s)')
ylabel('Es')
title('Es vs Time')
saveas(gcf,'Es vs Time,png')
%绘制阻尼耗能Ed时程
figure(3)
plot(Time,Ed)
xlabel('Time(s)')
ylabel('Ed')
title('Ed vs Time')
saveas(gcf,'Ed vs Time,png')
%绘制非弹性应变耗能Ep时程
figure(4)
plot(Time,Ep)
xlabel('Time(s)')
ylabel('Ep')
title('Ep vs Time')
saveas(gcf,'Ep vs Time,png')
%绘制耗能时程
E=Ep+Ed;
figure(5)
plot(Time,E)
xlabel('Time(s)')
ylabel('E')
title('E vs Time')
saveas(gcf,'E vs Time,png')
%保存数据
dlmwrite('Ek.csv',Ek,'delimiter','/n')
dlmwrite('Es.csv',Es,'delimiter','/n')
dlmwrite('Ed.csv',Ed,'delimiter','/n')
dlmwrite('Ep.csv',Ep,'delimiter','/n')

文件:单自由度体系的耗能计算.txt。计算过程参考崔济东老师的书籍《结构地震反应分析》。

2公茂盛，翟长海，谢礼立等，地震动滞回能量谱衰减规律研究