一、前言

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

作为学习三维黏弹性边界的前期探索，先从二维的简单模型入手。以《_人工边界及地震动输入在有限元软件中的实现》一文中给出的算例为参考，取其参数来进行验证黏弹性边界设置的正确性。

二、整体思路

建立模型（部件，材料，装配，分析步，荷载），导出inp文件

编写for程序；

用vs2013运行for程序，生成弹簧-阻尼器的语句；

再将语句放入模型inp文件对应的位置；

提交计算inp文件，后处理。

三、模型建立过程

1. 模型几何以及物理参数

模型范围.jpg

模型范围为800x400，网格尺寸为20x20

物理参数如下

参数.jpg

2. 分析步设置

总分析时间为20s，采用固定分析步长，取为0.005

3 荷载设置

荷载取原文，如下图

荷载.jpg

本次模拟的荷载取为分布荷载，荷载的位置如图，作用在自由面中间节点两边，4个单元网格长度。

荷载加载位置.jpg

四、弹簧-阻尼器设置

弹簧刚度和阻尼系数的计算，整体思路是

通过编写for程序；

用vs2013运行出弹簧-阻尼器的语句；

再将语句放入模型inp文件对应的位置；

最后提交计算inp，进行后处理。

1.编写for程序

01参数取值

• 散射源到人工边界的距离R的取值

到侧面的距离R取值为400

到底部面的距离R取值为400

• 边界修正系数取值

弹簧刚度和阻尼系数的计算公式不按原文计算的取值，按照文献《三维黏弹性静一动力统一人工边界》中推荐的公式，以及文献《黏弹性人工边界及地震动输入在通用有限元软件中的实现》中推荐的边界修正系数。

法向的修正系数取值为1，切向的修正系数取值为0.5

弹簧-阻尼示意图.jpg 修正参数取值.jpg

• 其他参数根据文献中物理参数表所给取值即可

节点控制面积单独说明

02 节点控制面积计算

1) 原理：边界节点等效面积的取值

根据文献《波动问题中的三维时域粘弹性人工边界》指出，

不同位置处的节点面积取值（边角1/4，边线1/2，中间1/1）

二维的模型，其节点的等效面积取值按如图规则：

节点等效面积取值.jpg

2）实现思路

01 提取某一个面的节点面积，该面则采用固定约束，并施加1的压强，之后提取节点反力，就是节点控制面积。

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

本模型是提取两个侧面、和底部面的节点反力。分别对每一个面采用固定约束，施加1的压强，提交计算，在odb里面report 反力（RF Magnitude），这样得到的就是节点面积。

3) 其他说明

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

代表储存生成的弹簧-阻尼器语句的文件

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

代表节点控制面积的文件

02 使用vs编译、链接、运行fortan程序

可参考一下操作

使用vs编译、链接、运行fortan程序

03 for程序完整代码(以两个侧面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为原始数据

!10-20为输出数据

!X左侧

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 -（侧面）=============

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两个方向边界和底部边界的弹簧-阻尼系数语句）同inp文件放入同一目录下
• 打开inp文件，定位end assembly,在之前插入语句，即能引用生成的黏弹性边界弹簧-阻尼系数语句。

*include, input=2dspringdashpotresults.txt

• 提交计算inp即可。

  
  inp中引用弹簧-阻尼语句.png

2.后处理

（1）自由面中点 u2方向位移时程曲线

自由面中点时程.jpg

（2）各个边界面 u2方向位移时程曲线

各个面中点位移时程.jpg

结果显示，算例设置的黏弹性边界的效果较好。和原文基本一致。

参考文献

• 《波动问题中的三维时域黏弹性人工边界》

• 《_人工边界及地震动输入在有限元软件中的实现》

• 《三维黏弹性静一动力统一人工边界》

• 《黏弹性人工边界及地震动输入在通用有限元软件中的实现》