第1页 / 共11页
第2页 / 共11页
第3页 / 共11页
第4页 / 共11页
第5页 / 共11页
第6页 / 共11页
第7页 / 共11页
第8页 / 共11页
ABAQUS粘弹性边界及地震荷载施加的简单实现.pdf
ABAQUS 粘弹性边界及地震荷载施加的简单实现
粘弹性边界因为能够考虑地基辐射阻尼而使得结构抗震的计算结果更趋于合理，所以在需要考虑结构地基
相互作用的结构抗震计算时，是较为常用的地基边界处理和地震荷载施加方法。而 ABAQUS 软件是经常用
来进行结构响应分析的有限元软件。下面介绍一种在 ABAQUS 中实现粘弹性边界及地震荷载施加的方法。
粘弹性边界是通过在有限元模型的地基边界节点上施加弹簧阻尼器实现的，在 ABAQUS 中的实现有以下几
种方法：第一种，通过 ABAQUS 自有的弹簧单元 spring 单元和阻尼单元 dashpot 实现，具体的单元参数可
以参考文献[1]，这种较为精确；第二种是通过 ABAQUS 的 UEL 子程序实现，可以看下文献[2]；还有一种是
等效单元替代的方法，就是在地基周围加一层单元，然后设置近似的材料参数，参考文献[3]，这一种精度
较差，但实现起来较为简单。我采用的是第一种方法，但操作起来较为繁琐，具体程序及过程后面介绍。
采用粘弹性边界，其配套的地震荷载输入方法就是在已知输入地震位移和速度的情况下，计算各个时刻地
基边界各个结点上应当施加的集中力荷载，然后施加荷载，一步一步的进行计算。地震荷载的施加在 ABAQUS
中也有两种不同的思路，文献[2]中的方法是通过 ABAQUS 的 DLOAD 和 UTRACLOAD 两个子程序实现。DLOAD
子程序用于施加边界面的法向荷载，UTRACLOAD 用于施加边界面的切向荷载。而文献[1]中则是将边界上每
一个节点每个时刻的力都计算出来，然后导入到 ABAQUS 中作为幅值数据，对每个对应节点施加。
我最初的想法是两篇文章的思路各取一半，用文献[1]的方法实现粘弹性边界，用文献[2]的方法施加地震荷
载。然而尝试了很久，发现这样做的效果并不是太好，可能我编的程序哪儿还是有问题吧。最后放弃了，
统一采用文献[1]的方法实现，具体实现采用 MATLAB 语言生成 ABAQUS 的 input 文件，然后将生成的 input
文件在模型文件的指定位置插入，用 ABAQUS 运行即可。
MATLAB 程序与运行前的准备工作
首先需要准备一些必要的数据文件（上图中红色框内的文件），其余黄色框内为模型文件，蓝色框内的文件
为程序运行后的生成文件，
boundary1~5.rpt 是从 ABAQUS 反力文件中提取的反力文件，其值代表某一节点的控制面积，可以通过在地
基边界（5 个面）施加值为 1 的压力荷载，即可提取得到这些反力文件；
coord_point.rpt 为 5 个边界面上节点的坐标文件，提取方法可以在很容易的百度到；
DIS.txt 是三个方向地震波的位移文件；
VEL.txt 是三个方向地震波的速度文件；
job-996.inp 是模型的 inp 文件；
Amplitude.inp 里面是计算过程中边界节点上随时间要施加的所有集中力荷载，文件较大；
load.inp 是将 Amplitude.inp 里面的幅值施加到对应节点的荷载命令；
springs&dashpot.inp 是模型地基边界施加的弹簧阻尼器文件；
三个 input 文件在模型 inp 文件中的插入位置：
记住 springs&dashpot.inp 在 Assembly 部分，所以搜索到关键字*End Assembly，把 springs&dashpot.inp 放在
*End Assembly 之前，Amplitude.inp 放在*End Assembly 之后，load.inp 放在 load 部分即可，如下所示
·································
·································
*Include, Input= springs&dashpot.inp
*End Assembly
*Include, Input=Amplitude.inp
·································
·································
** LOADS
**
*Include, Input=load.inp
**
** OUTPUT REQUESTS
以下为 MATLAB 程序，记得依据模型修改标红的参数准备并必要的文件
%%%%
%%%%-----------------------------说明------------------------------------
%%%%
% 1.本程序用于 ABAQUS 隐式计算的粘弹性边界的 inp 文件及荷载输入文件的生成
% 2.需要准备 5 个边界的节点反力文件（用于节点提取控制面积）、地震动位移和速度文件以及边
界节点坐标文件
% 3.输出为三个文件，springsanddashpot.inp 用于施加粘弹性边界;Amplitude.inp 为各节点幅
值;load.inp 为荷载文件
% 4.首先生成 Boundray_point，保存节点编号、节点坐标、控制面积
% 5.再生成弹簧阻尼器的 input 文件
% 6.再生成荷载文件等
% by w_tao
% 2018/10/01
%%%%---------------------------------------------------------------------
%%%%
%%%%-------------------------地基及相关计算参数---------------
%%%%
d_time=0.01;
%地震波的时间间隔，也将是计算步长
Z0=-600
%地基底面坐标，Z 方向为垂直方向
Z2=0
%地基地面坐标
H=Z2-Z0
%地基高度
X0=0
%地基水平 X 向坐标
X1=400
%地基水平 X 向坐标
LLLX=X1-X0
%地基 X 向宽度
Y0=0
%地基水平 Y 向坐标
Y1=400
%地基水平 Y 向坐标
LLLY=Y1-Y0
%地基 Y 向宽度
EEE=4.88e9;
%地基弹模
psb=0.22;
%地基泊松比
GGG=EEE/2/(1+psb);
%地基剪切模量
DENSITY=2000;
%地基密度
CP=sqrt((1-psb)*EEE/(1+psb)/DENSITY/(1-2*psb));
%计算地基中纵波波速
CS=sqrt(EEE/2/(1+psb)/DENSITY);
%计算地基中横波波速
lmt=(CP*CP-2*CS*CS)*DENSITY;
%拉梅常数中的 lmt
alpha_N=1.33;
%弹簧阻尼器系数
alpha_T=0.67;
%弹簧阻尼器系数
RRR(1)=LLLX/2;
RRR(2)=LLLY/2;
RRR(3)=H;
%%%%---------------------------------------------------------------------
%%%%
%%%%--------------------整理节点数据-------------------------------------
%%%%
%%%%加载文件数据
coord_point=load('coord_point.rpt')
boundary_Z0=load('boundary1.rpt')
boundary_X1=load('boundary2.rpt')
boundary_X0=load('boundary3.rpt')
boundary_Y1=load('boundary4.rpt')
boundary_Y0=load('boundary5.rpt')
% 加载地震波数据，默认位移和速度数据是一样长的
dis=load('DIS.txt');
vel=load('VEL.txt');
m=length(dis);
n(1)=length(boundary_Z0);%底面点数
n(2)=length(boundary_X1);%X 正向面的点数
n(3)=length(boundary_X0);%X 负向面的点数
n(4)=length(boundary_Y1);%y 正向面的点数
n(5)=length(boundary_Y0);%y 负向面的点数
n(6)=length(coord_point);%总节点数，但不等于 5 个面的和，因为有重复的点
error_point_num=0;
%coord_point 第 1 列节点号，第 2-4 列 XYZ 坐标，第 5 列归属的边界面号（1/2/3/4/5），第 6 列 R 值，第 7
列面积
%循环依据坐标确定归属号和 R 值
for i=1:n(6)
if(coord_point(i,4)==Z0)%底面点
coord_point(i,5)=1;
coord_point(i,6)=RRR(3);
elseif(coord_point(i,2)==X1)%X 正向面的点
coord_point(i,5)=2;
coord_point(i,6)=RRR(1);
elseif(coord_point(i,2)==X0)%X 负向面的点
coord_point(i,5)=3;
coord_point(i,6)=RRR(1);
elseif(coord_point(i,3)==Y1)%y 正向面的点
coord_point(i,5)=4;
coord_point(i,6)=RRR(2);
elseif(coord_point(i,3)==Y0)%y 负向面的点
coord_point(i,5)=5;
coord_point(i,6)=RRR(2);
else%理论上没有其他的点了，如果 error_point_num 非 0，说明节点坐标文件不正确
coord_point(i,5)=0;
coord_point(i,6)=0;
error_point_num=error_point_num+1;
end
end
%循环确定节点控制面积
for i=1:n(6)
if(coord_point(i,5)==1)
for j=1:n(1)
if(coord_point(i,1)==boundary_Z0(j,1))
coord_point(i,7)=abs(boundary_Z0(j,4));
end
end
elseif (coord_point(i,5)==2)
for j=1:n(2)
if(coord_point(i,1)==boundary_X1(j,1))
coord_point(i,7)=abs(boundary_X1(j,2));
end
end
elseif (coord_point(i,5)==3)
for j=1:n(3)
if(coord_point(i,1)==boundary_X0(j,1))
coord_point(i,7)=abs(boundary_X0(j,2));
end
end
elseif (coord_point(i,5)==4)
for j=1:n(4)
if(coord_point(i,1)==boundary_Y1(j,1))
coord_point(i,7)=abs(boundary_Y1(j,3));
end
end
elseif (coord_point(i,5)==5)
for j=1:n(5)
if(coord_point(i,1)==boundary_Y0(j,1))
coord_point(i,7)=abs(boundary_Y0(j,3));
end
end
else
error_point_num=error_point_num+1;
end
end
%%%%---------------------------------------------------------------------
%%%%
%%%%--------------------生成弹簧阻尼器和地震荷载的 input 文件--------------------------
%%%%
fid=fopen('springs&dashpot.inp','w')
fin_amp=fopen('Amplitude.inp','w')
fin_load=fopen('load.inp','w')
k_sum=0;
m_sum=0;
% 方向向量
for i=1:3
dri(i)=i;
end
%%%%
%循环每个节点
for i=1:n(6)
%依据节点属于哪个面，确定 alpha 和 C_vel 的顺序
if (coord_point(i,5)==1)%如果是底面（法向方向平行于 Z 面）的话
alpha(1)=alpha_T;
alpha(2)=alpha_T;
alpha(3)=alpha_N;
C_vel(1)=CS;
C_vel(2)=CS;
C_vel(3)=CP;
R=RRR(3);
elseif (coord_point(i,5)==2)%如果法向方向平行于 X 方向的话
alpha(1)=alpha_N;
alpha(2)=alpha_T;
alpha(3)=alpha_T;
C_vel(1)=CP;
C_vel(2)=CS;
C_vel(3)=CS;
R=RRR(1);
elseif (coord_point(i,5)==3)%如果法向方向平行于 X 方向的话
alpha(1)=alpha_N;
alpha(2)=alpha_T;
alpha(3)=alpha_T;
C_vel(1)=CP;
C_vel(2)=CS;
C_vel(3)=CS;
R=RRR(1)
elseif (coord_point(i,5)==4)%如果法向方向平行于 Y 方向的话
alpha(1)=alpha_T;
alpha(2)=alpha_N;
alpha(3)=alpha_T;
C_vel(1)=CS;
C_vel(2)=CP;
C_vel(3)=CS;
R=RRR(2);
elseif (coord_point(i,5)==5)%如果法向方向平行于 Y 方向的话
alpha(1)=alpha_T;
alpha(2)=alpha_N;
alpha(3)=alpha_T;
C_vel(1)=CS;
C_vel(2)=CP;
C_vel(3)=CS;
R=RRR(2);
else
error_point_num=error_point_num+1;
end
%%%% 计算每个节点的弹簧阻尼器参数，并写入文件中
%%%% X 方向弹簧阻尼器
k_sum=k_sum+1;
fprintf(fid,'%s%d\r\n','*Element, type=Spring1, elset=Springs-',k_sum)
m_sum=m_sum+1;
fprintf(fid,'%d%s%d\r\n',m_sum,', PART-1-1.',coord_point(i,1))
fprintf(fid,'%s%d\r\n','*Element, type=Dashpot1, elset=Dashpots-',k_sum)
m_sum=m_sum+1;
fprintf(fid,'%d%s%d\r\n',m_sum,', PART-1-1.',coord_point(i,1))
fprintf(fid,'%s%d\r\n','*Spring, elset=Springs-',k_sum)
fprintf(fid,'%d\r\n',dri(1))
KKK(1)=alpha(1)*GGG/R*coord_point(i,7);
fprintf(fid,'%d\r\n',KKK(1))
fprintf(fid,'%s%d\r\n','*Dashpot, elset=Dashpots-',k_sum)
fprintf(fid,'%d\r\n',dri(1))
CCC(1)=DENSITY*C_vel(1)*coord_point(i,7);
fprintf(fid,'%d\r\n',CCC(1))
%%%% Y 方向弹簧阻尼器
k_sum=k_sum+1;
fprintf(fid,'%s%d\r\n','*Element, type=Spring1, elset=Springs-',k_sum)
m_sum=m_sum+1;
fprintf(fid,'%d%s%d\r\n',m_sum,', PART-1-1.',coord_point(i,1))
fprintf(fid,'%s%d\r\n','*Element, type=Dashpot1, elset=Dashpots-',k_sum)
m_sum=m_sum+1;
fprintf(fid,'%d%s%d\r\n',m_sum,', PART-1-1.',coord_point(i,1))
fprintf(fid,'%s%d\r\n','*Spring, elset=Springs-',k_sum)
fprintf(fid,'%d\r\n',dri(2))
KKK(2)=alpha(2)*GGG/R*coord_point(i,7);
fprintf(fid,'%d\r\n',KKK(2))
fprintf(fid,'%s%d\r\n','*Dashpot, elset=Dashpots-',k_sum)
fprintf(fid,'%d\r\n',dri(2))
CCC(2)=DENSITY*C_vel(2)*coord_point(i,7);
fprintf(fid,'%d\r\n',CCC(2))
%%%% Z 方向弹簧阻尼器
k_sum=k_sum+1;
fprintf(fid,'%s%d\r\n','*Element, type=Spring1, elset=Springs-',k_sum)
m_sum=m_sum+1;
fprintf(fid,'%d%s%d\r\n',m_sum,', PART-1-1.',coord_point(i,1))
fprintf(fid,'%s%d\r\n','*Element, type=Dashpot1, elset=Dashpots-',k_sum)
m_sum=m_sum+1;
fprintf(fid,'%d%s%d\r\n',m_sum,', PART-1-1.',coord_point(i,1))
fprintf(fid,'%s%d\r\n','*Spring, elset=Springs-',k_sum)
fprintf(fid,'%d\r\n',dri(3))
KKK(3)=alpha(3)*GGG/R*coord_point(i,7);
fprintf(fid,'%d\r\n',KKK(3))
fprintf(fid,'%s%d\r\n','*Dashpot, elset=Dashpots-',k_sum)
fprintf(fid,'%d\r\n',dri(3))
CCC(3)=DENSITY*C_vel(3)*coord_point(i,7);
fprintf(fid,'%d\r\n',CCC(3))
%%%% 计算每个节点的集中力时程并写入文件中
Z1=coord_point(i,4)-Z0;
%%%% 时间循环
for j=1:m
time=(j-1)*d_time; %第一个时刻为 0
UNUM(1)=(time-Z1/CS)/d_time;
UNUM(2)=(time-(2*H-Z1)/CS)/d_time;
UNUM(3)=(time-Z1/CS)/d_time;
UNUM(4)=(time-(2*H-Z1)/CS)/d_time;
UNUM(5)=(time-Z1/CP)/d_time;
UNUM(6)=(time-(2*H-Z1)/CP)/d_time;
if(UNUM(1)>0)
K_NUM(1)=round(UNUM(1))+1;
U(1)=dis(K_NUM(1),1);
V(1)=vel(K_NUM(1),1);
else
K_NUM(1)=0;
U(1)=0;
V(1)=0;
end
if(UNUM(2)>0)
K_NUM(2)=round(UNUM(2))+1;
U(2)=dis(K_NUM(2),1);
V(2)=vel(K_NUM(2),1);
else
K_NUM(2)=0;
U(2)=0;
相关推荐
2023年江西萍乡中考道德与法治真题及答案.doc
2012年重庆南川中考生物真题及答案.doc
2013年江西师范大学地理学综合及文艺理论基础考研真题.doc
2020年四川甘孜小升初语文真题及答案I卷.doc
2020年注册岩土工程师专业基础考试真题及答案.doc
2023-2024学年福建省厦门市九年级上学期数学月考试题及答案.doc
2021-2022学年辽宁省沈阳市大东区九年级上学期语文期末试题及答案.doc
2022-2023学年北京东城区初三第一学期物理期末试卷及答案.doc
2018上半年江西教师资格初中地理学科知识与教学能力真题及答案.doc
2012年河北国家公务员申论考试真题及答案-省级.doc
2020-2021学年江苏省扬州市江都区邵樊片九年级上学期数学第一次质量检测试题及答案.doc
2022下半年黑龙江教师资格证中学综合素质真题及答案.doc
