粘弹性人工边界在ANSYS中的实现

粘弹性人工边界的有限元分析徐浩【摘要】采用数值方法集中比较了工程中广泛应用的几种人工边界在SSI体系分析中的动力反应,研究结果表明:粘性边界能较好的模拟土的边界但计算的位移会发生整体的飘逸,精度也不如粘弹性边界;粘弹性边界能很好的模拟土一结构相互作用体系中土的边界问题,但实现过程比较麻烦.【期刊名称】《山西建筑》【年(卷),期】2011(037)009【总页数】2页(P72-73)【关键词】人工边界;SSI体系;粘弹性人工边界;有限元法【作者】徐浩【作者单位】同济大学土木工程学院建筑工程系,上海,200092【正文语种】中文【中图分类】TU411土木工程中有许多涉及弹性半无限地基的振动及波动问题,诸如土—结构动力相互作用(SSI)问题、地震波的传播问题、动力机器基础的振动问题、打桩及道路交通引起的振动问题等。

对这类问题用有限元法求解与用其他计算方法相比具有可灵活地适用于地基物性的不均匀性并可考虑地基的非线性特性等优点。

但由于有限元法必须对求解对象的全领域进行离散,故在分析弹性半无限地基问题时必须截取一定范围的计算模型,这就要求在切取的边界上建立人工边界,用于模拟切除的无限域影响。

目前人工边界主要分为两类：一类是全局人工边界条件,如边界元法等;另一类是局部人工边界条件,如旁轴近似人工边界、透射人工边界、粘弹性人工边界等。

局部人工边界具有时空解耦的特点,得到了广泛的应用。

其中粘弹性人工边界是通过沿人工边界设置一系列由弹簧和阻尼器组成的简单物理元件来吸收射向人工边界的波动能量和反射波的散射,其模型简单,物理意义清晰,便于在大型通用有限元计算软件中得到实现。

ANSYS是一功能强大的有限元计算软件,其中的 ANSYS中的Combin14单元和 LSDYNA中的Combin165单元,是弹簧与阻尼器的元件,易于实现粘弹性人工边界。

本文在 ANSYS中实现了粘弹性边界并进行了验证和对比。

粘弹性人工边界作为一种应力边界条件,该应力是边界结点位移和速度的函数,一般形式写为：由式(1)可以看出,粘弹性边界相当于在边界结点每个方向施加一个一段固定的单向弹簧—阻尼元件,而且该弹簧—阻尼元件的刚度和阻尼系数仅与该边界结点在该方向该时刻的反应相关,从而通过粘性阻尼的吸能作用和弹簧的刚性恢复作用模拟无限域对广义结构的影响。

应用多种超弹性本构的案例：膝盖韧带分析
Cartilage(软骨) Isotropic Mooney-Rivlin
Ligament(韧带) Transversely isotropic Neo-Hookean
Meniscus(半月板)
Fung
W
2
I1
3
F2 (I4 )
Orthotropic Hyperelastic:
i 1
2i
3 i
3
第三类分子统计学本构模型，Neo-Hookean、Gent。
W
2
I1
3
5
超弹性模型 介绍
6
超弹性模型 选择
• 2阶项的 Mooney-Rivlin模型（ANSYS提供2,3,5,9参数的 模型）
W C1I1 3 C2 I2 3
W C1(I1 3) C2 (I2 3) C3 (I1 3)(I2 3)
Cauchy应力 剪切松弛模量 Green应变
t
0
2G t
de
d
d
t
I 0
K
t
d
d
d
式中σ——Cauchy应力(真实应力)
G(t)——为剪切松弛核函数
K(t)——为体积松弛核函数(由于体积几乎不变，可不考虑)
e——为应变偏量部分(剪切变形)
Δ——为应变体积部分(体积变形)
t——当前时间
τ——过去时间
粘弹性模型 静态 时温等效
• 时温等效原理 升高温度与延长时间对分子运动或高聚物的粘弹行为都
是等效的，这个等效性可以借助转换因子aT，将在某一温度 下测定的力学数据转换成另一温度下的数据.
静态下，升高温度与延长时间等效 降低温度与缩短时间等效

从半空间无限域取一4X2的矩形平面结构，顶部中间一定范围内受随时间变化的均布荷载，荷载如下p(t)=t 当0< DIV>p(t)=2-t 当1<=t<=2时p(t)=0 当t>2时材料弹性模量E=2.5，泊松比0.25，密度1网格尺寸0.1X0.1，在网格边界上所有结点加法向和切向combin14号单元用以模拟粘弹性人工边界（有关理论可参考刘晶波老师的相关文章）。

combine14单元的两个结点，其中一个与实体单元相连，另一个结点固定。

网格图如图1所示时程分析的时间步长为0.02秒，共计算16秒。

计算得到四个控制点位移时程图如图2所示，控制点坐标A(0,2)、B(0,1)、C(0,0)、D(2,2).计算所用命令流如下：/PREP7L=4 !水平长度H=2 !竖起深度E=2.5 !弹性模量density=1 !密度nu=0.25 !泊松比dxyz=0.1 !网格尺寸G = E/(2.*(1.+nu)) !剪切模量alfa = E*(1-nu)/((1.+nu)*(1.-2.*nu)) !若计算平面应力，此式需要修改Cp=sqrt(alfa/density) !压缩波速Cs=sqrt(g/density) !剪切波速R=sqrt(L*L/4.+H*H/4.) !波源到边界点等效长度KbT=0.5*G/R*dxyzKbN=1.0*G/R*dxyzCbT=density*Cs*dxyzCbN=density*Cp*dxyzET, 1, plane42,,,2 !按平面应变计算et, 2, combin14, ,, 2 !切向et, 3, combin14, ,, 2 !法向r, 2, KbT, CbTr, 3, KbN, CbNMP, EX, 1, EMP, PRXY, 1, nuMP, DENS, 1, densityrectng,-L/2.,L/2,0.,Hasel, allaesize, all, dxyzmshape,0,2Dmshkey,1amesh, all!以下建立底边界法向和切向弹簧阻尼单元nsel,s,loc,y,0.*get,np,node,,count !得到选中的结点数，存入np*get,npmax,node,,num,maxd !得到已经定义的最大结点数，存入npmax*do,ip,1,npnpnum=node((ip-1)*dxyz-L/2.,0.,0.)x=nx(npnum)y=ny(npnum)z=nz(npnum)npmax=npmax+1n,npmax,x.,y-dxyz/2,z !定义底边界法向结点以便与边界点形成法向单元type,3real,3e,npnum,npmaxd,npmax,all,0. !约束新生成的点npmax=npmax+1n,npmax,x-dxyz/2.,y,z !定义底边界切向结点以便与边界点形成切向单元type,2real,2e,npnum,npmaxd,npmax,all,0. !约束新生成的点*enddo!以下建立左边界法向和切向弹簧阻尼单元nsel,s,loc,x,-L/2*get,np,node,,count !得到选中的结点数，存入np*get,npmax,node,,num,maxd !得到已经定义的最大结点数，存入npmax*do,ip,2,np !侧边界最下面一个点按底边界上处理npnum=node(-L/2,(ip-1)*dxyz,0.)x=nx(npnum)y=ny(npnum)z=nz(npnum)npmax=npmax+1n,npmax,x-dxyz/2.,y,z !定义左边界法向结点以便与边界点形成法向单元type,3real,3e,npnum,npmaxd,npmax,all,0. !约束新生成的点npmax=npmax+1n,npmax,x,y-dxyz/2.,z !定义左边界切向结点以便与边界点形成切向单元type,2real,2e,npnum,npmaxd,npmax,all,0. !约束新生成的点*enddo!以下建立右边界法向和切向弹簧阻尼单元nsel,s,loc,x,L/2*get,np,node,,count !得到选中的结点数，存入np*get,npmax,node,,num,maxd !得到已经定义的最大结点数，存入npmax*do,ip,2,np !侧边界最下面一个点按底边界上处理npnum=node(L/2,(ip-1)*dxyz,0.)x=nx(npnum)y=ny(npnum)z=nz(npnum)npmax=npmax+1n,npmax,x+dxyz/2.,y,z !定义右边界法向结点以便与边界点形成法向单元type,3real,3e,npnum,npmaxd,npmax,all,0. !约束新生成的点npmax=npmax+1n,npmax,x,y-dxyz/2.,z !定义右边界切向结点以便与边界点形成切向单元type,2real,2e,npnum,npmaxd,npmax,all,0. !约束新生成的点*enddoallsel,all/pnum,type,1/number,1eplotfinish/soluANTYPE,trans!*TRNOPT,FULLLUMPM,0btime=0.02etime=16.00dtime=0.02*DO,itime,btime,etime,dtimeTIME,itimensel,s,loc,y,H !选中需要加荷载的点nsel,r,loc,x,-L/4,L/4*if,itime,lt,1.,thenf,all,fy,1*itime*elseif,itime,ge,1.0,and,itime,le,2.0f,all,fy,1*(2-itime)*elsef,all,fy,0.0*endifallsel,allSOLVE*ENDDO另外，还用自己编写的有限元程序计算了一下这个例子，并与ANSYS得到的结果进行了比较，结果非常吻合，这里给出A点的比较结果。

粘塑性
A. 粘塑性背景
• 与前述的率无关塑性相反, 率相关塑性与应变速率或时间有关。从 材料的观点看, 粘塑性 和蠕变是相同的。
– 通常, 对工程用途, 蠕变用于描述在恒载荷下应变的变化。通常当温度 达到材料熔点的30-60%时, 蠕变变形就很重要，而且, 时间更长，蠕 变和塑性应变解耦。
– ANSYS 中, Perzyna 和 Peirce 模型(TB,RATE) 意味着高应变率(即冲 击)载荷状态，非弹性应变是不解耦的。 – Anand 模型与 ANSYS 蠕变法则类似, 主要区别是使用了一个内部变 量，变形抗力, 来表示对非弹性材料流动的各向同性抗力。
1 o g in
n
Peirce 模型 n 1 o in g
n
g 1 in o in g
Strain Rates Time scale Temperature Effects
Supported Element Types 3
BISO, MISO, NLISO, BKIN, HILL BISO, MISO, NLISO, HILL None Suitable for small strain rates Suitable for large strain rates Suitable for small strain rates Long periods, creep and plasticity Short periods, usually for impactShort/medium periods have different time scales type problems Temperature effects included as Can input temperature-dependent Anand's equation considers part of equation (or material material constants, but equations do temperature effects directly. No constants can be temperaturenot consider temp effects directly need to input temperature2 dependent material constants dependent) Implicit - core and 18x Core and 18x VISCO106-108 Explicit - core and misc

& σ = SR (εin )σo
等效有效应力 应变率强化 [应力比率] (TB,RATE) 静态屈服应力 (TB,BISO/MISO/NLISO)
3
应力比率
ε &in m & SR (εin ) = 1+ γ & &in ) = 1+ εin SR (ε γ
m
PERZYNA
Perzyna 模型
2
1
PEIRCE
Peirce 模型
0 0 1 2 3 4 5
正则化应变率
September 30, 2001 Inventory #001491 5-8
粘塑性
... RATE 粘塑性选项
• 这意味着应力比率作为应变率的函数而改变。 这意味着应力比率作为应变率的函数而改变。
用 Perzyna 模型的例子
September 30, 2001 Inventory #001491 5-9
粘塑性
... RATE 粘塑性选项
• 关于 关于RATE模型的一些假设 模型的一些假设: 模型的一些假设
Training Manual
Advanced Structural Nonlinearities 6.0
September 30, 2001 Inventory #001491 5-3
粘塑性
... 粘塑性理论的背景
• ANSYS中率相关塑性模型的总结 中率相关塑性模型的总结: 中率相关塑性模型的总结
CREEP Behavior Isotropic or anisotropic creep (see HILL below) No explicitly defined yield surface RATE 1 Strain rate- or time-dependent Isotropic or anisotropic viscoplasticity (see HILL below) Includes yield surface ANAND Isotropic

防灾减灾工程学报第２７卷簧一阻尼器系统，图１为三维粘弹性人工边界示意图，人工边界上法向与切向的弹簧刚度和阻尼系数按照公式（１）和（２）取值‘８’９］。

，结构图１人工边界等效弹簧一阻尼器系统Ｋ删＝口Ⅳ号，ＣＢ．Ⅳ＝ｐｃｐ（１）广Ｋ盯＝奸百ｋＪ＂，Ｃ盯＝ｐｃ，（２）式中ＫｓＮ、Ｋ胛分别为法向与切向弹簧刚度；ＧⅣ、Ｃ胛分别为法向与切向阻尼器的阻尼系数；Ｒ为波源至人工边界点的距离；ｏ和以分别为介质的Ｐ波和Ｓ波波速；Ｇ为介质剪切模量；Ｐ为介质质量密度；ａＮ与即分别为法向与切向粘弹性人工边界修正系数。

大量数值计算表明，粘弹性人工边界具有良好的鲁棒性，人工边界参数口Ⅳ与奸在一定范围内取值均可以给出良好的计算结果，经过大量算例分析，推荐使用表１中的数据。

表ｌ粘弹性人工边界参数蛳与昕的取值范围粘弹性人工边界所模拟的是人工边界上的应力条件，因此这是１种连续分布的人工边界条件。

当采用有限元法或其它离散方法将人工边界所包围的计算区离散化时，人工边界面也将随之离散化，此时可以采用有限元法的形函数将连续分布的人工边界物理元件化为耦联的人工边界，可称为一致粘弹性人工边界；也可以简单地采用集中处理方法，形成解耦的人工边界，称为集中粘弹性人工边界。

三维粘弹性人工边界有限元离散系统如图２（ａ）所示，由三维单元边界面上的节点１、２、３、４所确定的平面上均匀分布着法向和切向粘弹性人工边界的阻尼器和弹簧。

（ａ）一致粘弹性人工边界（ｂ）一致粘弹性人工边界单元图２三维一致粘弹性人Ｔ边界及粘弹性人工边界单元示意１．２一致粘弹性人工边界单元为了更为简便地应用粘弹性人工边界，在有限元中使用等效实体单元来替换空间分布的弹簧一阻尼器元件，即在已有有限元模型的边界上沿边界面法向延伸一层一定厚度的实体单元。

并将外层边界固定，如图２（ｂ）所示‘１“。

如果用与计算区域相同的单元来模拟一致粘弹性边界，单元的等效剪切模量０和弹性模量雷为ｆｏ—ｈＫｖｒ＝的＾簧｛卺＝ｇ二≮｝铲＾Ｋ肼＝ｃ３）【蛳ＪＩｌ页Ｇ·型导％型式中Ｋ删、Ｋ朋．为人工边界弹簧法向与切向刚度；ｈ为边界单元的厚度；孝为等效泊松比，按下式取值：移一』者各口≥２移一＿２（口一１）一７。

粘塑性
A. 粘塑性背景
Training Manual
• 与前述的率无关塑性相反, 率相关塑性与应变速率或时间有关。从
材料的观点看, 粘塑性 和蠕变是相同的。
– 通常, 对工程用途, 蠕变用于描述在恒载荷下应变的变化。通常当温度 达到材料熔点的30-60%时, 蠕变变形就很重要，而且, 时间更长，蠕 变和塑性应变解耦。
Training Manual
记住首先定义线弹性材 料(EX 和 PRXY)。
否则, 当定义粘塑性时, 材料GUI将提示该信息。
第112页/共34页
Advanced Structural Nonlinearities 6.0
粘塑性
... 定义 RATE 材料特性
• 本例中选择了双线性等向强化的各向异性粘塑性。
• 当输入RATE常数后, 会出现率无关塑性参数对话框。
– 该例中, 输入双线性等向强化模型( BISO)的值, 即屈服应力和剪切模 量。
第143页/共34页
粘塑性
... 定义 RATE 材料特性
• 最后, 需要输入HILL参数来指定各向异性。
– 下面对话框中需要输入6个应力比率。
Training Manual
have different time scales
type problems
Temperature effects included as Can input temperature-dependent Anand's equation considers
part of equation (or material material constants, but equations do temperature effects directly. No

第 １７卷第 ５期 ２０１９年 １０月
水利与建筑工程学报
ＪｏｕｒｎａｌｏｆＷａｔｅｒＲｅｓｏｕｒｃｅｓａｎｄＡｒｃｈｉｔｅｃｔｕｒａｌＥｎｇｉｎｅｅｒｉｎｇ
ＤＯＩ：１０．３９６９／ｊ．ｉｓｓｎ．１６７２－１１４４．２０１９．０５．０１０
Ｖｏｌ．１７Ｎｏ．５ Ｏｃｔ．，２０１９
二维黏弹性人工边界单元及地震波输入 在 ＡＮＳＹＳ中的实现
对无限域或半无限域进行地震响应分析时，往 往采用有限元法。该方法通常采取的处理措施是在 无限域或半无限域空间中人为地截取一定范围的有 限域，进而可以将无限区域转化为有限区域，然后在 有限域的基础上，对其进行离散法处理，进而无限域
或半无限域最终被转化为有限个单元。同时为了达 到近似处理的效果，通常在截取的有限域的边界上 人为的施加约束，但是这一处理的缺点是在处理过 程中，不仅没有考虑远场无限地基的辐射阻尼效应， 同时半空间无限地基弹卷
截取的有限域的范围存在不确定性，随着有限域截 取范围的增加，自由度也随之增加。基于此离散范 围也更加宽广，进而网格划分的单元数及结点数成 倍增加，因此，在地震动的动力计算中工作量也进一 步被增加，对计算机的要求也更高，进而，在动力学 中非线性问题的求解难度进一步增大。目前，针对 该关键性问题处理的诸多措施中，较好的处理方法 是在截取的有限域边界处设置局部人工边界。目前 在工程中较常用的局部边界主要有黏性边界［１］、透 射边界［２］、黏弹性边界［３］等。黏性边界的优点是物 理概念清晰，在实际动力计算程序中实现起来相对 比较容易，对于 １阶动力响应精确度较高。基于此， 黏性边界单元被嵌入到了很多大型商用软件中以便 于工程的实际应用，如 ＬＳＤＹＮＡ、ＦＩＡＣ、ＡＢＡＱＵＳ［４］ 等。它的缺点是在施加人工边界过程中只着重考虑 了人工边界辐射阻尼效应，忽略了人工边界处半无 限空间介质的弹性恢复性能。因此在实际工程应用 中容易发生整体漂移的现象，进而使计算结果精度 不高。透射边界的优点是对于 ２阶动力响应计算结 果存在较高的精度，其缺点是动力计算分析程序的 代码编写较复杂、在实际工程应用中可能存在高频 失稳的现象。黏弹性边界不仅能有效的模拟半无限 空间介质的辐射阻尼效应，同时还能较好的模拟半 无限地基的弹性恢复性能。此外还具有良好的频率 稳定性和较高的精度。因此在实际工程中得到了广 泛的应用［５－１０］。

September 30, 2001 Inventory #001491 7-4
粘弹性
... 粘弹性理论背景
Training Manual
Advanced Structural Nonlinearities 6.0
• 粘弹性是率相关行为 材料特性可能与时间和温度都有关，粘弹性 粘弹性是率相关行为, 响应可看作由弹性和粘性部分组成。 响应可看作由弹性和粘性部分组成。
响应
t
ε = constant
(plus delta function)
ε
←δ(t-to)
加载
t
响应
t
September 30, 2001 Inventory #001491 7-14
粘弹性
... Kelvin-Voigt 模型
• Kelvin-Voigt 模型有以下特征 模型有以下特征:
– 对蠕变 在外载作用下 应变是指数的并渐近至 σo/E 对蠕变, 在外载作用下,
– 首先定义一些常用的术语 – 通过使用一维流变模型来解释线性粘弹性行为，这将有助于说明广义 通过使用一维流变模型来解释线性粘弹性行为， Maxwell 模型的基本特征。 模型的基本特征。 – 使用所有粘弹性材料通用的一般输入要求。 使用所有粘弹性材料通用的一般输入要求。 – 焦点将转至 WLF 偏移函数和聚合物。 偏移函数和聚合物。 – 将他讨论 TN偏移函数及其对玻璃材料的适用性。 偏移函数及其对玻璃材料的适用性。 偏移函数及其对玻璃材料的适用性
η
September 30, 2001 Inventory #001491 7-11
粘弹性
... Maxwell 模型
• Maxwell模型是一个串联的弹簧和阻尼器。 模型是一个串联的弹簧和阻尼器。 模型是一个串联的弹簧和阻尼器 σ spring = σ dashpot

粘弹性人工边界应用中的几个关键问题及其在ANSYS 中的实现蒋伟河海大学土木工程学院，江苏南京 (210098)E-mail: jw800403@摘 要：粘弹性人工边界能同时模拟半无限地基的能量辐射效应和弹性恢复能力，精度较高，计算结果稳定，在工程中受到越来越广泛的应用。

本文通过粘弹性人工边界理论，比较全面地介绍了粘弹性人工边界应用中人工边界的设置、参数选取、波动输入方法等几个关键问题以及在通用有限元分析软件ANSYS 中的实现，并结合平面问题算例，验证了该方法的有效性和准确性。

关键词：粘弹性人工边界；结构-地基动力相互作用；ANSYS ；波动输入1. 引言半无限地基的模拟问题是结构-地基动力相互作用分析中的一个关键问题。

目前通常的做法是在截取的有限域截断面上设置人工边界，合理地设置人工边界对于正确反映结构-地基的整体动力特性很重要。

人工边界大致可分为全局人工边界和局部人工边界两大类。

局部人工边界与全局人工边界相比，具有所需计算机存储量小、计算时间短、实用性强等优点，因此在实际工程中得到了比较广泛的应用。

局部人工边界中，工程上目前较常用的有廖振鹏等提出的透射边界[1]、Lysmer 等提出的粘性边界[2]，以及Deeks 在粘性边界的基础上提出了粘弹性人工边界[3]等。

透射边界虽具有较高精度，但在实际应用中一般仅限于二阶精度以内，并且存在编程较复杂、计算中可能引起高频失稳等问题。

粘性边界虽只有一阶精度，但概念清楚，易于程序实现，所以应用比较广泛，但其仅考虑了对散射波的吸收，不能模拟半无限地基的弹性恢复能力。

粘弹性边界具有能同时模拟散射波辐射和半无限地基的弹性恢复能力的优点，且能克服粘性边界引起的低频漂移问题，稳定性好。

目前，粘弹性人工边界已经开始应用到实际工程中，并越来越受到工程界的重视。

本文将以二维平面问题结合大型通用有限元计算软件ANSYS ，就粘弹性人工边界如何实现的几个问题做一简要的介绍。

粘弹性人工边界在有限元分析中的应用杜兴华;高扬【摘要】主要介绍了粘弹性人工边界的相关理论,牯弹性人工边界单元在ANSYS 中的具体实现方法,以及相关参数的计算公式.并且通过一个算例验证了粘弹性人工边界具有良好的计算精度和稳定性.【期刊名称】《低温建筑技术》【年(卷),期】2012(034)001【总页数】3页(P76-78)【关键词】粘弹性人工边界;ANSYS;地下结构【作者】杜兴华;高扬【作者单位】天津大学建筑工程学院,天津300072;天津大学建筑工程学院,天津300072【正文语种】中文【中图分类】TU471.2半无限地基的模拟是地下结构数值分析的一个关键问题。

目前普遍采用有限域来模拟半无限域，所以人工边界选取得是否合理，直接关系到数值分析的准确性。

目前较常用的有透射边界、粘性边界，粘弹性人工边界等。

透射边界虽具有较高精度，但编程较复杂、计算中可能引起高频失稳等问题［1］。

粘性边界概念清楚，易于程序实现，所以应用比较广泛，但其仅考虑了对散射波的吸收，不能模拟半无限地基的弹性恢复能力。

粘弹性边界具同时模拟散射波辐射和半无限地基的弹性恢复能力的优点，且能克服粘性边界引起的低频漂移问题，稳定性好。

目前，粘弹性人工边界已应用到实际工程中，并越来越受到工程界的重视。

1 粘弹性人工边界相关问题1.1 粘弹性人工边界理论粘弹性人工边界的推导过程同粘性边界相类似，在假设边界上不存在能量反射前提下，基于二维散射波为柱面波的情形可推导出任一半径γb处，以γb为外法线的微元面上应力同该处速度和位移的关系为：其中，G为剪切模量;ρ为介质密度;cs为介质中的剪切波速。

由式(1)可以看出，如果在半径rb处截断介质，并且在截断边界处施加等效的物理元件就可以消除波在人工边界处的反射。

由公式可知施加的物理元件为一个弹簧和一个阻尼器。

对于平面内波动问题，在人工边界的切线和法线两个方向上均需施加弹簧阻尼器，法线方向上的弹簧阻尼器值应从理论上重新推导，但可以将G和cs 简单地用E和cp替换。

粘弹性人工边界在ANSYS中的实现作者：河海水妖 2007-11-07 00:25:58标签：知识/探索ansys粘弹性人工边界动力边界条件粘弹性人工边界在ANSYS中的实现从半空间无限域取一4X2的矩形平面结构，顶部中间一定范围内受随时间变化的均布荷载，荷载如下p(t)=t 当0< DIV>p(t)=2-t 当1<=t<=2时p(t)=0 当t>2时材料弹性模量E=2.5，泊松比0.25，密度1网格尺寸0.1X0.1，在网格边界上所有结点加法向和切向combin14号单元用以模拟粘弹性人工边界（有关理论可参考刘晶波老师的相关文章）。

combine14单元的两个结点，其中一个与实体单元相连，另一个结点固定。

网格图如图1所示时程分析的时间步长为0.02秒，共计算16秒。

计算得到四个控制点位移时程图如图2所示，控制点坐标A(0,2)、B(0,1)、C(0,0)、D(2,2).计算所用命令流如下：/PREP7L=4 !水平长度H=2 !竖起深度E=2.5 !弹性模量density=1 !密度nu=0.25 !泊松比dxyz=0.1 !网格尺寸G = E/(2.*(1.+nu)) !剪切模量alfa = E*(1-nu)/((1.+nu)*(1.-2.*nu)) !若计算平面应力，此式需要修改 Cp=sqrt(alfa/density) !压缩波速Cs=sqrt(g/density) !剪切波速R=sqrt(L*L/4.+H*H/4.) !波源到边界点等效长度KbT=0.5*G/R*dxyzKbN=1.0*G/R*dxyzCbT=density*Cs*dxyzCbN=density*Cp*dxyzET, 1, plane42,,,2 !按平面应变计算et, 2, combin14, ,, 2 !切向et, 3, combin14, ,, 2 !法向r, 2, KbT, CbTr, 3, KbN, CbNMP, EX, 1, EMP, PRXY, 1, nuMP, DENS, 1, densityrectng,-L/2.,L/2,0.,Hasel, allaesize, all, dxyzmshape,0,2Dmshkey,1amesh, all!以下建立底边界法向和切向弹簧阻尼单元nsel,s,loc,y,0.*get,np,node,,count !得到选中的结点数，存入np*get,npmax,node,,num,maxd !得到已经定义的最大结点数，存入npmax*do,ip,1,npnpnum=node((ip-1)*dxyz-L/2.,0.,0.)x=nx(npnum)y=ny(npnum)z=nz(npnum)npmax=npmax+1n,npmax,x.,y-dxyz/2,z !定义底边界法向结点以便与边界点形成法向单元 type,3real,3e,npnum,npmaxd,npmax,all,0. !约束新生成的点npmax=npmax+1n,npmax,x-dxyz/2.,y,z !定义底边界切向结点以便与边界点形成切向单元 type,2real,2e,npnum,npmaxd,npmax,all,0. !约束新生成的点*enddo!以下建立左边界法向和切向弹簧阻尼单元nsel,s,loc,x,-L/2*get,np,node,,count !得到选中的结点数，存入np*get,npmax,node,,num,maxd !得到已经定义的最大结点数，存入npmax*do,ip,2,np !侧边界最下面一个点按底边界上处理npnum=node(-L/2,(ip-1)*dxyz,0.)x=nx(npnum)y=ny(npnum)z=nz(npnum)npmax=npmax+1n,npmax,x-dxyz/2.,y,z !定义左边界法向结点以便与边界点形成法向单元 type,3real,3e,npnum,npmaxd,npmax,all,0. !约束新生成的点npmax=npmax+1n,npmax,x,y-dxyz/2.,z !定义左边界切向结点以便与边界点形成切向单元 type,2real,2e,npnum,npmaxd,npmax,all,0. !约束新生成的点*enddo!以下建立右边界法向和切向弹簧阻尼单元nsel,s,loc,x,L/2*get,np,node,,count !得到选中的结点数，存入np*get,npmax,node,,num,maxd !得到已经定义的最大结点数，存入npmax*do,ip,2,np !侧边界最下面一个点按底边界上处理npnum=node(L/2,(ip-1)*dxyz,0.)x=nx(npnum)y=ny(npnum)z=nz(npnum)npmax=npmax+1n,npmax,x+dxyz/2.,y,z !定义右边界法向结点以便与边界点形成法向单元 type,3real,3e,npnum,npmaxd,npmax,all,0. !约束新生成的点npmax=npmax+1n,npmax,x,y-dxyz/2.,z !定义右边界切向结点以便与边界点形成切向单元 type,2real,2e,npnum,npmaxd,npmax,all,0. !约束新生成的点*enddoallsel,all/pnum,type,1/number,1eplotfinish/soluANTYPE,trans!*TRNOPT,FULLLUMPM,0btime=0.02etime=16.00dtime=0.02*DO,itime,btime,etime,dtimeTIME,itimensel,s,loc,y,H !选中需要加荷载的点nsel,r,loc,x,-L/4,L/4*if,itime,lt,1.,thenf,all,fy,1*itime*elseif,itime,ge,1.0,and,itime,le,2.0f,all,fy,1*(2-itime)*elsef,all,fy,0.0*endifallsel,allSOLVE*ENDDO另外，还用自己编写的有限元程序计算了一下这个例子，并与ANSYS得到的结果进行了比较，结果非常吻合，这里给出A点的比较结果。

粘弹性
B. 流变模型
Training Manual
• 在详细讨论 ANSYS 中可用的粘弹性选项之前, 对一些常用的流变 模型(1-D)进行总结有助于理解粘弹性行为。
• 下面的讨论依赖于两个基本模型, 弹簧和缓冲器
– 弹簧采用 Hooke 定律作为应力和应变的关系，这代表弹性固体。
s
E
– 缓冲器(或阻尼器)定义应力和应变率间的行为, 用于表示粘性流体。
– Maxwell 模型充分地描述了应力松弛而不是蠕变。 – Kelvin-Voigt 模型充分地描述了蠕变而不是应力松弛。 – SLM 提供了一个简单的蠕变和松弛行为的表示法。
• 速率效应如下:
– 快速加载极限: sE0E1e – 缓慢加载极限: sE0e
September 30, 2001 Inventory #001491 7-17
– 这是一些高阶单元(能使用退化形式)，虽然可能, 但不推荐用它们作为 低阶单元(减去中间节点)。
– 目前 2D 平面应力、壳和梁单元对粘弹性不可用。 – VISCO88/89 单元有应力刚化能力, 但假设为小应变、小位移行为。
September 30, 2001 Inventory #001491 7-19
Advanced Structural Nonlinearities 6.0
粘弹性
... 广义 Maxwell 模型
Training Manual
• 广义 Maxwell 模型是由k 个并联的弹簧 和缓冲器组成.
G0
G1
h1
k
s si i 0
s i

si i

G
i e
i

hi Gi

1.粘弹性：ANSYS中的粘弹性模型是Maxwell模型的通用积分形式，其松弛函数由Prony级数表示。

该模型功能全面，Maxwell、Kevin和标准线性实体都是其特殊形式，全面支持亚粘弹性和大应变超粘弹性。

大应变超粘弹性基于Simo建议的列式，粘弹性行为的定义分为超弹性和松弛两个部分，所有的ANSYS超弹性材料模型都可采用粘弹性选项（PRONY）。

2. 粘弹性是率相关行为, 材料特性可能与时间和温度都有关，粘弹性响应可看作由弹性和粘性部分组成。

–弹性部分是可恢复的, 且是瞬时的。

–粘性部分是不可恢复的, 且在整个时间范围内发生。

ANSYS 中能模拟线性粘弹性，这导致如下假设:
–应变率与瞬态应力成比例
–瞬态应变与瞬态应力也成比例
–限于小应变、小变形行为(NLGEOM,OFF)
–
C5=1
FICT TEMP可以从帮助文件里找到
注意密度。

图 6 时间-荷载示意图
图 7 变形曲线图
5.4 模拟结果
荷载作用 20 次后永久变形见图 7，横坐标为水平方向，纵坐标为深度方向，可以看到其最大永 久变形为 6.85mm，并车辙现象形成相似成 U 型。
图 8 水平方向应力图示
图 10 Y 方向轮边缘应力分布
图 9 Y 方向轮中心应力分布
图 11 荷载中心的竖向位移时间曲线图
∑ G (t)
=
G∞
+
nG i =1
Gi
exp
⎛ ⎜ ⎝
−
t
τ
G i
⎞ ⎟ ⎠
∑ （2）
K
(t)
=
K∞
+
nK i =1
Ki
exp
⎛ ⎜ ⎝
−
t
τ
K i
⎞ ⎟ ⎠
（3）
式中，G∞
和 Gi
是剪切模量， K∞
和
Ki
是体积模量，τ
G i
和τ
K i
是各
Prony
级数分量的松弛时间。
再定义下面相对模量：
αiG = Gi G0
（4）
α
K i
=
Ki
K0
（5）
式中， G0 ， K0 分别为粘弹性材质的瞬态模量，并定义式如下：
nG
∑ G0 = G (t = 0) = G∞ + Gi i =1
nK
∑ （6） K0 = K (t = 0) = K∞ + Ki i =1
（7）
4. 材料参数的实验数据
在 ANSYS 里有这么一组材料数据：剪切模量与时间的关系如表 1：
时间（s） 0.01 1 2 4 6

ANSYS 粘弹性材料1.1 ANSYS 中表征粘弹性属性问题粘弹性材料的应力响应包括弹性部分和粘性部分，在载荷作用下弹性部分是即时响应的，而粘性部分需要经过一段时间才能表现出来。

一般的，应力函数是由积分形式给出的，在小应变理论下，各向同性的粘弹性本构方程可以写成如下形式：()()002t t de d G t d I K t d d d σττττττ∆=-+-⎰⎰ (1) 其中σ＝Cauchy 应力()G t ＝为剪切松弛核函数()K t ＝为体积松弛核函数e ＝为应变偏量部分（剪切变形）∆＝为应变体积部分（体积变形）t ＝当前时间τ＝过去时间I ＝为单位张量。

该式是根据松弛条件本构方程(1)，通过将一点的应变分解为应变球张量（体积变形）和应变斜张量（剪切变形）两部分，推导而得的。

这里不再敖述，可参考相关文献等。

ANSYS 中描述粘弹性积分核函数()G t 和()K t 参数表示方式主要有两种，一种是广义Maxwell 单元（VISCO88 和 VISCO89）所采用的Maxwell 形式，一种是结构单元所采用的Prony 级数形式。

实际上，这两种表示方式是一致的，只是具体数学表达式有一点点不同。

1.2 Prony 级数形式用Prony 级数表示粘弹性属性的基本形式为：()1exp G n i G i it G t G G τ∞=⎛⎫=+- ⎪⎝⎭∑ (2) ()1exp K n i K i i t K t K K τ∞=⎛⎫=+- ⎪⎝⎭∑ (3)其中，G ∞和i G 是剪切模量，K ∞和i K 是体积模量，G i τ和Ki τ是各Prony 级数分量的松弛时间(Relative time)。

再定义下面相对模量(Relative modulus) 0G i i G G α= (4)0K i i K K α= (5)其中，0G ，0K 分别为粘弹性材质的瞬态模量，并定义式如下：()010Gn i i G G t G G ∞====+∑ (6)()010Kn i i K K t K K ∞====+∑ (7)在ANSYS 中，Prony 级数的阶数G n 和K n 可以不必相同，当然其中的松弛时间G i τ和K i τ也不必相同。

基于ABAQUS的粘弹性动力人工边界精确自动施加徐磊;叶志才;任青文;杜小凯【摘要】粘弹性动力人工边界在结构-地基动力相互作用问题的有限元分析中已得到了广泛的应用,但对于规模较大的复杂计算模型而言,粘弹性动力人工边界的施加存在着前处理工作量浩繁、施加精度不高等问题.为此,基于粘弹性动力人工边界理论及ABAQUS软件平台,提出了粘弹性动力人工边界精确自动施加方法,编制了相关程序并进行了验证.算例将本文方法和通常的近似施加方法进行了对比分析,结果表明,近似施加方法会引起较大计算误差,本文所提出方法的施加简便性、施加精度均优于通常的近似施加方法.【期刊名称】《三峡大学学报（自然科学版）》【年(卷),期】2010(032)001【总页数】4页(P20-23)【关键词】粘弹性动力人工边界;ABAQUS;自动施加【作者】徐磊;叶志才;任青文;杜小凯【作者单位】河海大学,水利水电工程学院,南京,210098;河海大学,土木工程学院,南京,210098;宿迁市水务局,江苏,宿迁,223800;河海大学,土木工程学院,南京,210098;中国水电工程顾问集团公司,北京,100120【正文语种】中文【中图分类】P315在对结构-地基动力相互作用问题的计算分析中,如何考虑无限地基的辐射阻尼效应是核心问题之一[1] .若人工截取有限区域地基来模拟无限地基,则会使得散射波在人工截取的边界上产生反射而导致模拟失真.目前,在截取的有限域边界上设置合理的动力人工边界条件是解决这一问题的有效方法[2] .近年来,相关学者已提出了多种动力人工边界条件,如粘性边界、旁轴近似边界、透射边界和粘弹性边界等[3] .其中,粘弹性动力人工边界[ 4] 以其精度较高、能模拟人工边界外半无限介质弹性恢复性能且有良好的频率稳定性而在结构-地基动力相互作用问题的有限元分析中得到了广泛的应用.对于实际工程中的结构-地基动力相互作用有限元分析,一般均需采用三维计算模型,随着工程仿真建模的日益精细化,有限元分析模型存在大量的边界结点,数量通常可以达到数千甚至更多.当采用粘弹性动力人工边界时,由于各边界结点与散射波源(结构)之间距离不等、各结点所控制的边界面积不等以及各点所在区域的地基材料也会有所差别,因而,对于各边界结点,需逐点计算其弹簧系数及阻尼系数并施加弹簧单元及阻尼器单元,有限元前处理工作量浩繁.为了便于处理,通常采用一些近似方法来进行粘弹性边界的施加,如各边界结点取散射波源至其所在人工边界面的平均距离作为其与散射波源实际距离的近似等[5] ,这些近似势必会对边界条件的精度造成一定程度的影响.为此,基于粘弹性动力人工边界理论大型商用有限元软件ABAQUS,开发相关程序,实现了对粘弹性动力人工边界的精确自动施加.算例分析验证了程序的合理性和有效性.此外,还将本文方法和通常的近似施加方法进行了对比分析,结果表明,近似施加方法会引起较大计算误差,本文所提出方法的施加简便性、施加精度均优于通常的近似施加方法.1 粘弹性动力人工边界粘弹性动力人工边界可以等效为在计算模型截断边界上连续分布的并联弹簧-阻尼器系统,其中,弹簧系数K b 及阻尼系数Cb 的计算公式如下:式中,G、ρ分别为人工边界处介质的剪切模量、密度;R 为散射波源至人工边界的距离;c 为人工边界处介质的波速,法向阻尼系数取为纵波速cp,切向阻尼系数取为剪切波速cs;α的取值取决于粘弹性动力人工边界的类型及设置方向,见表1.表1 α的取值类型方向α二维平面内切向1.5出平面切向0.5平面内法向2.0三维切向2.0法向4.0粘弹性动力人工边界是一种连续分布的人工边界条件.当采用有限元法将计算区域离散化后,连续的人工边界面也随之离散化,此时,在有限元计算模型中,为了实现粘弹性动力人工边界的施加,需要在边界结点的法向和切向分别设置并联的弹簧单元和阻尼器单元,相应的弹簧系数K 和阻尼系数C 除了与式(1)有关外,还与边界结点的等效控制面积相关,计算表达式如下式中,Gi 、ρi 、ci 及Ai 分别为与边界结点相关联单元的材料物理力学参数及面积. 对于实际工程问题,由于有限元计算分析模型的边界结点数量众多,且需逐点按式(2)计算其弹簧系数及阻尼系数并施加弹簧单元及阻尼器单元,有限元前处理工作十分困难,若采用近似处理又难免引入误差.为此,在上述基础上,基于ABAQ US 编制了用于粘弹性动力人工边界精确自动施加的fortran 程序.2 粘弹性动力人工边界精确自动施加近年来,作为国际上最先进的大型通用商业非线性有限元分析软件之一,A BAQ US 在结构-地基动力相互作用等问题中的应用日益广泛[ 6] .ABAQUS 中的弹簧阻尼器单元也为施加粘弹性动力人工边界提供了方便.此外,用户还可以通过二次开发接口[7] 或编程读写input 计算文件等方式实现更为复杂的计算功能.为此,基于粘弹性动力人工边界理论以及ABAQ US 的input 计算文件的相关格式,编制了用以实现粘弹性动力人工边界精确自动施加的fortran程序.该程序在获取有限元计算模型相关信息的基础上,直接生成用以施加粘弹性动力人工边界的命令文本,该文本完全符合input 计算文件的相关格式,将其放入相应input 文件的指定位置即可完成基于ABAQ US 的粘弹性动力人工边界精确自动施加.为了实现粘弹性动力人工边界精确自动施加,首先需要完成有限元计算模型的网格剖分、材料分区等,在此基础上,获取所需的有限元网格的结点坐标、各单元的结点编号及其材料号、需要施加粘弹性动力人工边界的边界结点编号、各边界结点所在边界面的法线方向以及散射波源位置等信息.在上述基础上,编制了相应的fortran 程序,实现了基于ABAQUS 的粘弹性动力人工边界条件精确自动施加,图1 给出了相应的程序流程图.图1 程序流程图下面对程序中较为关键的问题进行扼要的说明,首先是如何计算关联单元边界面(即位于截断边界上的单元面)的面积,一般而言,与某一边界结点相关联的单元有多个,对任一关联单元,可首先依据单元各结点坐标与边界结点坐标之间的关系,并结合边界结点所在面的法向来确定在单元边界面的所有结点.其次是如何计算单元边界面的面积,对于二维问题,边界面的面积计算实质上是计算两点间线段的长度,易于解决;对于三维问题,单元边界面形状一般为三角形或四边形,对于三角形,依其顶点坐标容易计算其面积,而对于四边形,当只知其4 个顶点坐标而不知其顶点顺序时,无法应用通常的四边形面积公式进行求解,此时,可先应用海伦公式计算由此4 个顶点中任意3 个顶点所确定的4 个不同三角形的面积,则四边形面积即为这4 个三角形面积之和的一半.需要说明的是,以上关于面积计算的说明中仅考虑不含中结点的单元,对于非此类单元,一般均可将其单元形态近似为不含中结点的单元来加以处理.3 算例验证及对比分析为了验证2 中所编程序的正确性和有效性,进行了如下算例分析.考虑半无限空间z ≥0,分析模型如图2 所示,人工截取计算区域取为|x|≤10 m 、|y|≤10 m、|z|≤10 m.模型介质的物理力学参数如下:弹性模量E =50 MPa;泊松比v=0.3;剪切波速cs =83.62 m/s;纵波速cp=156.48 m/s;密度ρ=2 750 kg/m3.为了便于对计算成果进行分析,在分析模型上设置了2 个观察点A、B,其坐标分别为(0,0,10),(5,0,10).在自由面上沿z 方向加一集中面源(作用范围为-5 m ≤x ≤5 m,-5 m ≤y ≤5 m,z=10 m)时变荷载f(t), f(t)的函数表达式如式(3)所示.式中,0 ≤t ≤2 s.图2 算例分析模型采用8 结点六面体单元进行网格剖分,有限元模型如图3 所示,网格尺寸Δx=Δy=Δz=1 m,散射波源取为荷载作用面的中心.在截断边界上,采用节2中编制的fortran 程序实现粘弹性动力人工边界的精确自动施加,记为C.B.;为了对比分析本文方法与通常近似方法之间的差别,还采用了后一方法对上述算例进行了粘弹性边界的施加,在施加过程中,各边界结点取散射波源至其所在人工边界面的平均距离作为其与散射波源实际距离的近似,记为A.B.;此外,还给出了远置粘弹性动力人工边界(本文采取在模型各个方向均扩大20 倍范围后设置人工边界的方法)的有限元解作为精确解,记为E.B..图3 有限元计算网格图动力时程分析的时间步长取为0.01 s,计算总时长为3 s.图4 ～5 给出了观测点A、B 在不同边界条件下的z 向位移反应时程曲线.从中可以看出,与精确解相比,本文所提出的粘弹性动力人工边界精确自动施加方法(C.B.)与近似施加方法(A.B.)均可以给出与精确解(E.B.)总体趋势一致的计算结果,但也可明显发现,C.B.方法所给出的计算结果与A.B.方法所给出的计算结果相比,与精确解更为一致,而A.B.方法所给出的计算结果与精确解相比有较大误差,这是由于其在设置粘弹性动力人工边界时所采用的简化近似处理而造成的.以上分析表明,本文所提出的粘弹性动力人工边界精确自动施加方法在计算精度上相比通常的近似施加方法有明显的优势,可以给出更为精确的计算结果,此外,由于本文所提出方法为自动施加,其在施加速度上较通常近似施加方法也有明显的优势,尤其是对于计算规模较大、边界结点较多的结构与地基动力相互作用分析问题.4 结语在结构-地基动力相互作用问题的有限元分析中,粘弹性动力人工边界得到了广泛的应用.但对于规模较大的复杂计算模型而言,粘弹性动力人工边界的施加存在着前处理工作量浩繁、施加精度不高等问题.为此, 本文基于粘弹性动力人工边界理论及ABAQ US 软件平台,提出了粘弹性动力人工边界精确自动施加方法,编制了相关程序并进行了验证.算例将本文方法和通常的近似施加方法进行了对比分析,结果表明,近似施加方法会引起较大计算误差,本文所提出方法的施加简便性、施加精度均优于通常的近似施加方法.参考文献:[1] 刘晶波, 谷音,杜义欣.一致粘弹性动力人工边界及粘弹性边界单元[ J] .岩土工程学报, 2006, 28(9):1070-1075.[2] 刘晶波,王振宇,杜修力等.波动问题中的三维时域粘弹性动力人工边界[ J] .工程力学,2005,22(6):46-51.[3] 廖振鹏.工程波动理论导论[ M] .2 版.北京:科学出版社, 2002.[4] Deeks A J, Randolph M F.Axisymmetric Time-domain Transmitting Boundaries[ J] .Journal of Engineering Mechanics, 1994,120(1):25-42. [5] 张燎军,张慧星, 王大胜等.黏弹性人工边界在ADINA中的应用[ J] .世界地震工程,2008,24(1):12-16.[6] 庄海洋,陈国兴,胡晓明.两层双柱岛式地铁车站结构水平向非线性地震反应分析[ J] .岩石力学与工程学报,2006, 25(S1):3074-3079.[7] 叶志才, 徐磊,王超.基于ABAQUS 的三维锚杆单元开发[ J] .三峡大学学报:自然科学版, 2008, 30(5):29-32.。