ABAQUS粘弹性边界及地震荷载施加的简单实现.docx

1、ABAQUS粘弹性边界及地震荷载施加的简单实现ABAQUS粘弹性边界及地震荷载施加的简单实现(Matlab生成input文件)思路粘弹性边界因为能够考虑地基辐射阻尼而使得结构抗震的计算结果更趋于合理，所以在需要考虑结构地基相互作用的结构抗震计算时，是较为常用的地基边界处理和地震荷载施加方法。而ABAQUS软件是经常用来进行结构响应分析的有限元软件。下面介绍一种在ABAQUS中实现粘弹性边界及地震荷载施加的方法。粘弹性边界是通过在有限元模型的地基边界节点上施加弹簧阻尼器实现的，在ABAQUS中的实现有以下几种方法：第一种，通过ABAQUS自有的弹簧单元spring单元和阻尼单元dashpot实现

2、，具体的单元参数可以参考文献1，这种较为精确；第二种是通过ABAQUS的UEL子程序实现，可以看下文献2；还有一种是等效单元替代的方法，就是在地基周围加一层单元，然后设置近似的材料参数，参考文献3，这一种精度较差，但实现起来较为简单。我采用的是第一种方法，但操作起来较为繁琐，具体程序及过程后面介绍。采用粘弹性边界，其配套的地震荷载输入方法就是在已知输入地震位移和速度的情况下，计算各个时刻地基边界各个结点上应当施加的集中力荷载，然后施加荷载，一步一步的进行计算。地震荷载的施加在ABAQUS中也有两种不同的思路，文献2中的方法是通过ABAQUS的DLOAD和UTRACLOAD两个子程序实现。DLO

3、AD子程序用于施加边界面的法向荷载，UTRACLOAD用于施加边界面的切向荷载。而文献1中则是将边界上每一个节点每个时刻的力都计算出来，然后导入到ABAQUS中作为幅值数据，对每个对应节点施加。我最初的想法是两篇文章的思路各取一半，用文献1的方法实现粘弹性边界，用文献2的方法施加地震荷载。然而尝试了很久，发现这样做的效果并不是太好，可能我编的程序哪儿还是有问题吧。最后放弃了，统一采用文献1的方法实现，具体实现采用MATLAB语言生成ABAQUS的input文件，然后将生成的input文件在模型文件的指定位置插入，用ABAQUS运行即可。MATLAB程序与运行前的准备工作首先需要准备一些必要的数

4、据文件（上图中红色框内的文件），其余黄色框内为模型文件，蓝色框内的文件为程序运行后的生成文件，boundary15.rpt是从ABAQUS反力文件中提取的反力文件，其值代表某一节点的控制面积，可以通过在地基边界（5个面）施加值为1的压力荷载，即可提取得到这些反力文件；coord_point.rpt为5个边界面上节点的坐标文件，提取方法可以在很容易的XX到；DIS.txt是三个方向地震波的位移文件；VEL.txt是三个方向地震波的速度文件；job-996.inp是模型的inp文件；Amplitude.inp里面是计算过程中边界节点上随时间要施加的所有集中力荷载，文件较大；load.inp是将Am

5、plitude.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, Inpu

6、t=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，保存节点编号、

7、节点坐标、控制面积% 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);%地基剪切模量D

8、ENSITY=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;%拉梅常数中的lmtalpha_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(bo

9、undary1.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

10、负向面的点数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);el

11、seif(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)=R

12、RR(2);else%理论上没有其他的点了，如果error_point_num非0，说明节点坐标文件不正确coord_point(i,5)=0;coord_point(i,6)=0;error_point_num=error_point_num+1;endend%循环确定节点控制面积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);endendelseif (coord_point(i,5)=2)for j

13、=1:n(2)if(coord_point(i,1)=boundary_X1(j,1)coord_point(i,7)=abs(boundary_X1(j,2);endendelseif (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);endendelseif (coord_point(i,5)=4)for j=1:n(4)if(coord_point(i,1)=boundary_Y1(j,1)coord_point(i,7)=ab

14、s(boundary_Y1(j,3);endendelseif (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);endendelseerror_point_num=error_point_num+1;endend%-%-生成弹簧阻尼器和地震荷载的input文件-%fid=fopen(springs&dashpot.inp,w)fin_amp=fopen(Amplitude.inp,w)fin_load=fopen(load.inp,

15、w)k_sum=0;m_sum=0;% 方向向量for i=1:3dri(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)=alp

16、ha_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_v

17、el(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);elseerror_point_num=error_point_num+1;end% 计算每个节点的弹簧阻尼器参数，并写入文件中% X方向弹簧阻尼器k_sum=k_sum+1;fprintf(fid,%s%drn,*Element, typ

18、e=Spring1, elset=Springs-,k_sum)m_sum=m_sum+1;fprintf(fid,%d%s%drn,m_sum, PART-1-1.,coord_point(i,1)fprintf(fid,%s%drn,*Element, type=Dashpot1, elset=Dashpots-,k_sum)m_sum=m_sum+1;fprintf(fid,%d%s%drn,m_sum, PART-1-1.,coord_point(i,1)fprintf(fid,%s%drn,*Spring, elset=Springs-,k_sum)fprintf(fid,%drn,

19、dri(1)KKK(1)=alpha(1)*GGG/R*coord_point(i,7);fprintf(fid,%drn,KKK(1)fprintf(fid,%s%drn,*Dashpot, elset=Dashpots-,k_sum)fprintf(fid,%drn,dri(1)CCC(1)=DENSITY*C_vel(1)*coord_point(i,7);fprintf(fid,%drn,CCC(1)% Y方向弹簧阻尼器k_sum=k_sum+1;fprintf(fid,%s%drn,*Element, type=Spring1, elset=Springs-,k_sum)m_sum=

20、m_sum+1;fprintf(fid,%d%s%drn,m_sum, PART-1-1.,coord_point(i,1)fprintf(fid,%s%drn,*Element, type=Dashpot1, elset=Dashpots-,k_sum)m_sum=m_sum+1;fprintf(fid,%d%s%drn,m_sum, PART-1-1.,coord_point(i,1)fprintf(fid,%s%drn,*Spring, elset=Springs-,k_sum)fprintf(fid,%drn,dri(2)KKK(2)=alpha(2)*GGG/R*coord_poin

21、t(i,7);fprintf(fid,%drn,KKK(2)fprintf(fid,%s%drn,*Dashpot, elset=Dashpots-,k_sum)fprintf(fid,%drn,dri(2)CCC(2)=DENSITY*C_vel(2)*coord_point(i,7);fprintf(fid,%drn,CCC(2)% Z方向弹簧阻尼器k_sum=k_sum+1;fprintf(fid,%s%drn,*Element, type=Spring1, elset=Springs-,k_sum)m_sum=m_sum+1;fprintf(fid,%d%s%drn,m_sum, PA

22、RT-1-1.,coord_point(i,1)fprintf(fid,%s%drn,*Element, type=Dashpot1, elset=Dashpots-,k_sum)m_sum=m_sum+1;fprintf(fid,%d%s%drn,m_sum, PART-1-1.,coord_point(i,1)fprintf(fid,%s%drn,*Spring, elset=Springs-,k_sum)fprintf(fid,%drn,dri(3)KKK(3)=alpha(3)*GGG/R*coord_point(i,7);fprintf(fid,%drn,KKK(3)fprintf(

23、fid,%s%drn,*Dashpot, elset=Dashpots-,k_sum)fprintf(fid,%drn,dri(3)CCC(3)=DENSITY*C_vel(3)*coord_point(i,7);fprintf(fid,%drn,CCC(3)% 计算每个节点的集中力时程并写入文件中Z1=coord_point(i,4)-Z0;% 时间循环for j=1:mtime=(j-1)*d_time; %第一个时刻为0UNUM(1)=(time-Z1/CS)/d_time;UNUM(2)=(time-(2*H-Z1)/CS)/d_time;UNUM(3)=(time-Z1/CS)/d_

24、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);elseK_NUM(1)=0;U(1)=0;V(1)=0;endif(UNUM(2)0)K_NUM(2)=round(UNUM(2)+1;U(2)=dis(K_NUM(2),1);V(2)=vel(K_NUM(2),1);elseK_NUM

25、(2)=0;U(2)=0;V(2)=0;endif(UNUM(3)0)K_NUM(3)=round(UNUM(3)+1;U(3)=dis(K_NUM(3),2);V(3)=vel(K_NUM(3),2);elseK_NUM(3)=0;U(3)=0;V(3)=0;endif(UNUM(4)0)K_NUM(4)=round(UNUM(4)+1;U(4)=dis(K_NUM(4),2);V(4)=vel(K_NUM(4),2);elseK_NUM(4)=0;U(4)=0;V(4)=0;endif(UNUM(5)0)K_NUM(5)=round(UNUM(5)+1;U(5)=dis(K_NUM(5),

26、3);V(5)=vel(K_NUM(5),3);elseK_NUM(5)=0;U(5)=0;V(5)=0;endif(UNUM(6)0)K_NUM(6)=round(UNUM(6)+1;U(6)=dis(K_NUM(6),3);V(6)=vel(K_NUM(6),3);elseK_NUM(6)=0;U(6)=0;V(6)=0;end%依据节点属于哪个面，确定集中力荷载中的自由场应力项sigmaif (coord_point(i,5)=1)%如果是底面（法向方向平行于Z面）的话sigma(1)=DENSITY*CS*(V(1)-V(2)*coord_point(i,7);sigma(2)=DEN

27、SITY*CS*(V(3)-V(4)*coord_point(i,7);sigma(3)=DENSITY*CP*(V(5)-V(6)*coord_point(i,7);elseif (coord_point(i,5)=2)%如果法向方向平行于X方向的话sigma(1)=-lmt/CP*(V(5)-V(6)*coord_point(i,7);sigma(2)=0;sigma(3)=-DENSITY*CS*(V(1)-V(2)*coord_point(i,7);elseif (coord_point(i,5)=3)%如果法向方向平行于X方向的话sigma(1)=lmt/CP*(V(5)-V(6)*coord_point(i,7);sigma(2)=0;sigma(3)=DENSITY*