ABAQUS粘彈性邊界及地震荷載施加的簡單實現(Matlab生成input文件)

 思路

粘彈性邊界由於可以考慮地基輻射阻尼而使得結構抗震的計算結果更趨於合理，因此在須要考慮結構地基相互做用的結構抗震計算時，是較爲經常使用的地基邊界處理和地震荷載施加方法。而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;

V(2)=0;

end

if(UNUM(3)>0)

K_NUM(3)=round(UNUM(3))+1;

U(3)=dis(K_NUM(3),2);

V(3)=vel(K_NUM(3),2);

else

K_NUM(3)=0;

U(3)=0;

V(3)=0;

end

if(UNUM(4)>0)

K_NUM(4)=round(UNUM(4))+1;

U(4)=dis(K_NUM(4),2);

V(4)=vel(K_NUM(4),2);

else

K_NUM(4)=0;

U(4)=0;

V(4)=0;

end

if(UNUM(5)>0)

K_NUM(5)=round(UNUM(5))+1;

U(5)=dis(K_NUM(5),3);

V(5)=vel(K_NUM(5),3);

else

K_NUM(5)=0;

U(5)=0;

V(5)=0;

end

if(UNUM(6)>0)

K_NUM(6)=round(UNUM(6))+1;

U(6)=dis(K_NUM(6),3);

V(6)=vel(K_NUM(6),3);

else

K_NUM(6)=0;

U(6)=0;

V(6)=0;

end

 

%依據節點屬於哪一個面，肯定集中力荷載中的自由場應力項sigma

if (coord_point(i,5)==1)%若是是底面（法向方向平行於Z面）的話

sigma(1)=DENSITY*CS*(V(1)-V(2))*coord_point(i,7);

sigma(2)=DENSITY*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*CS*(V(1)-V(2))*coord_point(i,7);

elseif (coord_point(i,5)==4)%若是法向方向平行於Y方向的話

sigma(1)=0;

sigma(2)=-lmt/CP*(V(5)-V(6))*coord_point(i,7);

sigma(3)=-DENSITY*CS*(V(3)-V(4))*coord_point(i,7);

elseif (coord_point(i,5)==5)%若是法向方向平行於Y方向的話

sigma(1)=0;

sigma(2)=lmt/CP*(V(5)-V(6))*coord_point(i,7);

sigma(3)=DENSITY*CS*(V(3)-V(4))*coord_point(i,7);

else

error_point_num=error_point_num+1;

end

 

%節點3個方向的集中力時程

amp(j,1)=KKK(1)*(U(1)+U(2))+CCC(1)*(V(1)+V(2))+sigma(1);

amp(j,2)=KKK(2)*(U(3)+U(4))+CCC(2)*(V(3)+V(4))+sigma(2);

amp(j,3)=KKK(3)*(U(5)+U(6))+CCC(3)*(V(5)+V(6))+sigma(3);

end

 

for j=1:3

%將節點3個方向的集中力時程寫入文件中

fprintf(fin_amp,'%s%d%s%d\r\n','*Amplitude, name=Amp-',coord_point(i,1),'-',dri(j))

for k=1:m

fprintf(fin_amp,'%f%s%f\r\n',k*d_time,',   ',amp(k,j))

end

end

 

for j=1:3

%生成集中力荷載施加的文件

fprintf(fin_load,'%s%d%s%d\r\n','*Cload, amplitude=Amp-',coord_point(i,1),'-',dri(j))

fprintf(fin_load,'%s%d%s%d%s\r\n','PART-1-1.',coord_point(i,1),', ',dri(j),', 1.')

end

 

end

status=fclose(fid);

status=fclose(fin_amp);

status=fclose(fin_load);

算例驗證

以文獻[2]中的模型做爲驗證模型，具體參數見下圖：

 

模型以下圖所示：

驗證模型400m×400m×600m，紫色爲彈簧阻尼器

計算結果以下所示

Z方向地基地面點和地基上部點的位移對比

X方向地基地面點和地基上部點的位移對比

Y方向地基地面點和地基上部點的位移對比

 

入射波幅值爲0.5m，地面放大2倍接近1m，反射波0.5m

 

參考文獻

[1]何建濤, 馬懷發, 張伯豔,等. 黏彈性人工邊界地震動輸入方法及實現[J]. 水利學報, 2010, 39(8):960-969.

[2]苑舉衛, 杜成斌, 陳燈紅. 基於ABAQUS的三維粘彈性邊界單元及地震動輸入方法研究[J]. 三峽大學學報(天然科學版), 2010, 32(3):9-13.

[3]谷音, 劉晶波, 杜義欣. 三維一致粘彈性人工邊界及等效粘彈性邊界單元[J]. 工程力學, 2007, 24(12):31-37.

[4]潘堅文. ABAQUS水利工程應用實例教程[M]. 中國建築工業出版社, 2015.