一、前言

** 2d-visco-plastic-boundary&pressure test**

作為學(xué)習(xí)三維黏彈性邊界的前期探索，先從二維的簡(jiǎn)單模型入手。以《_人工邊界及地震動(dòng)輸入在有限元軟件中的實(shí)現(xiàn)》一文中給出的算例為參考，取其參數(shù)來進(jìn)行驗(yàn)證黏彈性邊界設(shè)置的正確性。

二、整體思路

建立模型（部件，材料，裝配，分析步，荷載），導(dǎo)出inp文件

編寫for程序；

用vs2013運(yùn)行for程序，生成彈簧-阻尼器的語句；

再將語句放入模型inp文件對(duì)應(yīng)的位置；

提交計(jì)算inp文件，后處理。

三、模型建立過程

1. 模型幾何以及物理參數(shù)

模型范圍.jpg

模型范圍為800x400，網(wǎng)格尺寸為20x20

物理參數(shù)如下

參數(shù).jpg

2. 分析步設(shè)置

總分析時(shí)間為20s，采用固定分析步長(zhǎng)，取為0.005

3 荷載設(shè)置

荷載取原文，如下圖

荷載.jpg

本次模擬的荷載取為分布荷載，荷載的位置如圖，作用在自由面中間節(jié)點(diǎn)兩邊，4個(gè)單元網(wǎng)格長(zhǎng)度。

荷載加載位置.jpg

四、彈簧-阻尼器設(shè)置

彈簧剛度和阻尼系數(shù)的計(jì)算，整體思路是

通過編寫for程序；

用vs2013運(yùn)行出彈簧-阻尼器的語句；

再將語句放入模型inp文件對(duì)應(yīng)的位置；

最后提交計(jì)算inp，進(jìn)行后處理。

1.編寫for程序

01參數(shù)取值

• 散射源到人工邊界的距離R的取值

到側(cè)面的距離R取值為400

到底部面的距離R取值為400

• 邊界修正系數(shù)取值

彈簧剛度和阻尼系數(shù)的計(jì)算公式不按原文計(jì)算的取值，按照文獻(xiàn)《三維黏彈性靜一動(dòng)力統(tǒng)一人工邊界》中推薦的公式，以及文獻(xiàn)《黏彈性人工邊界及地震動(dòng)輸入在通用有限元軟件中的實(shí)現(xiàn)》中推薦的邊界修正系數(shù)。

法向的修正系數(shù)取值為1，切向的修正系數(shù)取值為0.5

彈簧-阻尼示意圖.jpg

修正參數(shù)取值.jpg

• 其他參數(shù)根據(jù)文獻(xiàn)中物理參數(shù)表所給取值即可

節(jié)點(diǎn)控制面積單獨(dú)說明

02 節(jié)點(diǎn)控制面積計(jì)算

1) 原理：邊界節(jié)點(diǎn)等效面積的取值

根據(jù)文獻(xiàn)《波動(dòng)問題中的三維時(shí)域粘彈性人工邊界》指出，

不同位置處的節(jié)點(diǎn)面積取值（邊角1/4，邊線1/2，中間1/1）

二維的模型，其節(jié)點(diǎn)的等效面積取值按如圖規(guī)則：

節(jié)點(diǎn)等效面積取值.jpg

2）實(shí)現(xiàn)思路

01 提取某一個(gè)面的節(jié)點(diǎn)面積，該面則采用固定約束，并施加1的壓強(qiáng)，之后提取節(jié)點(diǎn)反力，就是節(jié)點(diǎn)控制面積。

02 再用語句** open(unit=1,file='2dx+.txt')、** 引用到 for程序中。

本模型是提取兩個(gè)側(cè)面、和底部面的節(jié)點(diǎn)反力。分別對(duì)每一個(gè)面采用固定約束，施加1的壓強(qiáng)，提交計(jì)算，在odb里面report 反力（RF Magnitude），這樣得到的就是節(jié)點(diǎn)面積。

3) 其他說明

open(unit=10,file='2dx+result.txt')

代表儲(chǔ)存生成的彈簧-阻尼器語句的文件

open(unit=1,file='2dx+.txt')

代表節(jié)點(diǎn)控制面積的文件

02 使用vs編譯、鏈接、運(yùn)行fortan程序

可參考一下操作

使用vs編譯、鏈接、運(yùn)行fortan程序

03 for程序完整代碼(以兩個(gè)側(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

四、提交計(jì)算&后處理

1.提交計(jì)算

• 將for程序生成的文件【2dx+result】【2dx-result】【2dyresult】（分別代表x兩個(gè)方向邊界和底部邊界的彈簧-阻尼系數(shù)語句）同inp文件放入同一目錄下
• 打開inp文件，定位end assembly,在之前插入語句，即能引用生成的黏彈性邊界彈簧-阻尼系數(shù)語句。

*include, input=2dspringdashpotresults.txt

• 提交計(jì)算inp即可。

  
  
  inp中引用彈簧-阻尼語句.png

2.后處理

（1）自由面中點(diǎn) u2方向位移時(shí)程曲線

自由面中點(diǎn)時(shí)程.jpg

（2）各個(gè)邊界面 u2方向位移時(shí)程曲線

各個(gè)面中點(diǎn)位移時(shí)程.jpg

結(jié)果顯示，算例設(shè)置的黏彈性邊界的效果較好。和原文基本一致。

參考文獻(xiàn)

• 《波動(dòng)問題中的三維時(shí)域黏彈性人工邊界》

• 《_人工邊界及地震動(dòng)輸入在有限元軟件中的實(shí)現(xiàn)》

• 《三維黏彈性靜一動(dòng)力統(tǒng)一人工邊界》

• 《黏彈性人工邊界及地震動(dòng)輸入在通用有限元軟件中的實(shí)現(xiàn)》