思路
粘彈性邊界因為能夠考慮地基輻射阻尼而使得結構抗震的計算結果更趨於合理，所以在需要考慮結構地基相互作用的結構抗震計算時，是較為常用的地基邊界處理和地震荷載施加方法。而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.