一晕鹊、前言
** 2d-visco-plastic-boundary&pressure test**
作為學(xué)習三維黏彈性邊界的前期探索玛臂，先從二維的簡單模型入手画舌。以《_人工邊界及地震動輸入在有限元軟件中的實現(xiàn)》一文中給出的算例為參考兵拢，取其參數(shù)來進行驗證黏彈性邊界設(shè)置的正確性孕锄。
二吮廉、整體思路
建立模型（部件，材料畸肆，裝配宦芦，分析步，荷載）轴脐，導(dǎo)出inp文件
編寫for程序调卑；
用vs2013運行for程序，生成彈簧-阻尼器的語句大咱；
再將語句放入模型inp文件對應(yīng)的位置恬涧；
提交計算inp文件，后處理碴巾。
三溯捆、模型建立過程
1. 模型幾何以及物理參數(shù)
模型范圍為800x400，網(wǎng)格尺寸為20x20
物理參數(shù)如下
2. 分析步設(shè)置
總分析時間為20s厦瓢，采用固定分析步長提揍，取為0.005
3 荷載設(shè)置
荷載取原文啤月，如下圖
本次模擬的荷載取為分布荷載，荷載的位置如圖碳锈，作用在自由面中間節(jié)點兩邊顽冶，4個單元網(wǎng)格長度欺抗。
四售碳、彈簧-阻尼器設(shè)置
彈簧剛度和阻尼系數(shù)的計算，整體思路是
通過編寫for程序绞呈；
用vs2013運行出彈簧-阻尼器的語句贸人；
再將語句放入模型inp文件對應(yīng)的位置；
最后提交計算inp佃声，進行后處理艺智。
1.編寫for程序
01參數(shù)取值
- 散射源到人工邊界的距離R的取值
到側(cè)面的距離R取值為400
到底部面的距離R取值為400
- 邊界修正系數(shù)取值
彈簧剛度和阻尼系數(shù)的計算公式不按原文計算的取值，按照文獻《三維黏彈性靜一動力統(tǒng)一人工邊界》中推薦的公式圾亏，以及文獻《黏彈性人工邊界及地震動輸入在通用有限元軟件中的實現(xiàn)》中推薦的邊界修正系數(shù)十拣。
法向的修正系數(shù)取值為1，切向的修正系數(shù)取值為0.5
- 其他參數(shù)根據(jù)文獻中物理參數(shù)表所給取值即可
節(jié)點控制面積單獨說明
02 節(jié)點控制面積計算
1) 原理：邊界節(jié)點等效面積的取值
根據(jù)文獻《波動問題中的三維時域粘彈性人工邊界》指出志鹃，
不同位置處的節(jié)點面積取值（邊角1/4夭问，邊線1/2，中間1/1）
二維的模型曹铃，其節(jié)點的等效面積取值按如圖規(guī)則：
2）實現(xiàn)思路
01 提取某一個面的節(jié)點面積缰趋，該面則采用固定約束，并施加1的壓強陕见，之后提取節(jié)點反力秘血，就是節(jié)點控制面積。
02 再用語句** open(unit=1,file='2dx+.txt')评甜、** 引用到 for程序中灰粮。
本模型是提取兩個側(cè)面、和底部面的節(jié)點反力忍坷。分別對每一個面采用固定約束粘舟，施加1的壓強，提交計算承匣，在odb里面report 反力（RF Magnitude）蓖乘，這樣得到的就是節(jié)點面積。
3) 其他說明
open(unit=10,file='2dx+result.txt')
代表儲存生成的彈簧-阻尼器語句的文件
open(unit=1,file='2dx+.txt')
代表節(jié)點控制面積的文件
02 使用vs編譯韧骗、鏈接嘉抒、運行fortan程序
可參考一下操作
03 for程序完整代碼(以兩個側(cè)面x-些侍、x+為例)
program main
implicit none
integer::i,s,node
real::density,G,R,alfaN,alfaT,mu,E,lamta,Cp,Cs&
&,KBN,KBT,CBN,CBT,A,y
!write(*,*)'輸入密度'
!read(*,*)density
!write(*,*)'輸入剪切模量G和泊松比'
!read(*,*)G,mu
!write(*,*)'輸入散射源到人工邊界的距離'
!read(*,*)R
!write(*,*)'輸入alfaN和alfaT'
!read(*,*)alfaN,alfaT
R=400
density=2000
E=1250000000
mu=0.25
G=0.5*E/(1+mu)
lamta=E*mu/((1+mu)*(1-2*mu))
Cp=sqrt((lamta+2*G)/density)
Cs=sqrt(G/density)
KBN=1.0*G/R
KBT=0.5*G/R
CBN=density*Cp
CBT=density*Cs
write(*,*)E,lamta,Cp,Cs,KBN,KBT,CBN,CBT
!1-10為原始數(shù)據(jù)
!10-20為輸出數(shù)據(jù)
!X左側(cè)
open(unit=10,file='2dx+result.txt')
open(unit=1,file='2dx+.txt')
s=0
do i=1,10000
read(1,*,end=500)node,a
!法向x
write(10,50)node
50 format('*Spring, elset="X+N',i0,'-spring"')
write(99,150)node
150 format('X+N',i0,'-spring')
write(10,51)
51 format('1')
write(10,52)KBN*a
52 format(e10.3)
write(10,53)node
53 format('*Dashpot, elset="X+N',i0,'-dashpot"')
write(99,151)node
151 format('X+N',i0,'-dashpot')
write(10,54)
54 format('1')
write(10,55)CBN*a
55 format(e10.3)
write(10,56)node
56 format('*Element, type=Spring1, elset="X+N',i0,'-spring"')
s=s+1
write(10,57)s,node
57 format(i0,',',' Part-1-1','.',i0)
write(10,58)node
58 format('*Element, type=Dashpot1, elset="X+N',i0,'-dashpot"')
s=s+1
write(10,59)s,node
59 format(i0,',',' Part-1-1','.',i0)
!切向y
write(10,60)node
60 format('*Spring, elset="X+YT',i0,'-spring"')
write(99,152)node
152 format('X+YT',i0,'-spring')
write(10,61)
61 format('2')
write(10,62)KBT*a
62 format(e10.3)
write(10,63)node
63 format('*Dashpot, elset="X+YT',i0,'-dashpot"')
write(99,153)node
153 format('X+YT',i0,'-dashpot')
write(10,64)
64 format('2')
write(10,65)CBT*a
65 format(e10.3)
write(10,66)node
66 format('*Element, type=Spring1, elset="X+YT',i0,'-spring"')
s=s+1
write(10,67)s,node
67 format(i0,',',' Part-1-1','.',i0)
write(10,68)node
68 format('*Element, type=Dashpot1, elset="X+YT',i0,'-dashpot"')
s=s+1
write(10,69)s,node
69 format(i0,',',' Part-1-1','.',i0)
end do
500 continue
!================X -（側(cè)面）=============
open(unit=12,file='2dx-result.txt')
open(unit=2,file='2dx-.txt')
do i=1,10000
read(2,*,end=501)node,a
!法向x
write(12,80)node
80 format('*Spring, elset="X-N',i0,'-spring"')
write(99,154)node
154 format('X-N',i0,'-spring')
write(12,81)
81 format('1')
write(12,82)KBN*a
82 format(e10.3)
write(12,83)node
83 format('*Dashpot, elset="X-N',i0,'-dashpot"')
write(99,155)node
155 format('X-N',i0,'-dashpot')
write(12,84)
84 format('1')
write(12,85)CBN*a
85 format(e10.3)
write(12,86)node
86 format('*Element, type=Spring1, elset="X-N',i0,'-spring"')
s=s+1
write(12,87)s,node
87 format(i0,',',' Part-1-1','.',i0)
write(12,88)node
88 format('*Element, type=Dashpot1, elset="X-N',i0,'-dashpot"')
s=s+1
write(12,89)s,node
89 format(i0,',',' Part-1-1','.',i0)
!切向y
write(12,90)node
90 format('*Spring, elset="X-YT',i0,'-spring"')
write(99,156)node
156 format('X-YT',i0,'-spring')
write(12,91)
91 format('2')
write(12,92)KBT*a
92 format(e10.3)
write(12,93)node
93 format('*Dashpot, elset="X-YT',i0,'-dashpot"')
write(99,157)node
157 format('X-YT',i0,'-dashpot')
write(12,94)
94 format('2')
write(12,95)CBT*a
95 format(e10.3)
write(12,96)node
96 format('*Element, type=Spring1, elset="X-YT',i0,'-spring"')
s=s+1
write(12,97)s,node
97 format(i0,',',' Part-1-1','.',i0)
write(12,98)node
98 format('*Element, type=Dashpot1, elset="X-YT',i0,'-dashpot"')
s=s+1
write(12,99)s,node
99 format(i0,',',' Part-1-1','.',i0)
end do
501 continue
end
四隶症、提交計算&后處理
1.提交計算
- 將for程序生成的文件【2dx+result】【2dx-result】【2dyresult】（分別代表x兩個方向邊界和底部邊界的彈簧-阻尼系數(shù)語句）同inp文件放入同一目錄下
- 打開inp文件，定位end assembly,在之前插入語句岗宣，即能引用生成的黏彈性邊界彈簧-阻尼系數(shù)語句蚂会。
*include, input=2dspringdashpotresults.txt
-
提交計算inp即可。
inp中引用彈簧-阻尼語句.png
2.后處理
（1）自由面中點 u2方向位移時程曲線
（2）各個邊界面 u2方向位移時程曲線
結(jié)果顯示耗式，算例設(shè)置的黏彈性邊界的效果較好胁住。和原文基本一致。
參考文獻
《波動問題中的三維時域黏彈性人工邊界》
《_人工邊界及地震動輸入在有限元軟件中的實現(xiàn)》
《三維黏彈性靜一動力統(tǒng)一人工邊界》
《黏彈性人工邊界及地震動輸入在通用有限元軟件中的實現(xiàn)》
