柔性容器简化建模方法，有效减少模型箱效应

1.柔性容器模拟

在振动台试验中，只能用有限尺寸的容器来装模型土。容器边界上的波动反射以及体系振动形态的变化将会给试验结果带来一定的误差，即所谓“模型箱效应”。为了减少模型箱效应，试验设计了一个柔性容器，详见第2.3节。

对土-结构相互作用振动台模型试验进行建模时，应反映柔性容器的性状。建模时容器侧壁采用三维壳单元SHELL63划分。容器底部与振动台台面用螺栓连接，且由于试验中采取适当措施使容器底板粗糙，可以忽略土与容器底部之间的相对滑移，建模时将土体底部考虑为固定端。

圆筒外侧的钢筋加固，在建模时比较难实现。若按实际情况将钢筋一圈圈地加上去，建模就会相当复杂。考虑到钢筋加固的目的是提供径向刚度，且允许土体作层状水平剪切变形，建模时将其考虑为：同一高度处沿圆筒周边的节点在沿地震波输入方向上(振动台X方向)具有相同的位移，并利用ANSYS软件中的自由度耦合功能来实现[128]。采用这样的建模方法，建模比较简单，而且能够满足精度的要求。本书第3.3节验证了这种建模方法的合理性。

2.网格划分

在对相互作用体系进行网格划分时，要考虑波动对网格划分的影响，考虑模型的节点与试验的测点相对应，并要采用合适的网格单元尺寸，具体原则如下：

(1)满足波动对网格划分的要求。如果单元尺寸过大，则波动的高频部分难以通过。研究表明[84]，对于一般沿竖向传播的剪切波，单元高度可取为：

式中：vs——剪切波速(m/s2)；

f max——截取的最大波动频率(Hz)。

单元的平面尺寸比高度尺寸的限制要求不甚严格，它取决于土层的情况，一般取h max的3～5倍。

(2)满足计算仿真与试验结果对照分析对网格划分的要求。使节点与试验的测点布置相对应，便于计算结果与试验结果进行对照分析。

(3)满足求解精度对有限元网格划分的要求。单元网格越细，自由度数越多，计算精度越高，但计算所需时间和计算量也会很大。因而，应采用合适的网格单元尺寸。

3.土体的动力本构模型与材料非线性模拟

土体本构模型的选取，不仅要考虑其能否比较真实地反映在给定环境下土的物理-力学特征，而且还要考虑是否可以获得该模型中所包含的参数的可靠数值。另外，还要考虑计算机的容量、速度以及计算费用与计算结果有效性之间的关系。土体动力本构模型主要有：双线性模型、等效线性模型、Iwan模型、Martin-Finn-Seed模型和内时模型等。在描述土性非线性方面，等效线性模型具有概念明确、应用方便等优点，在实际中应用最广泛。

采用等效线性模型时土层的划分一般应遵循下列原则：性质接近的场地土划为一层，但每层层厚不宜过大，否则由于土层深度变化引起的应力、位移变化难以确切反映。通常土层的卓越周期越长需要分割的数目越多。I.Midriss等根据集中质点系模型按不同分割的计算结果与剪切梁法得到的解析解进行比较，给出土层基本自振周期与分割数N的关系曲线(图3-1)，可作为划分土层时的参考。

图3-1　土层剪切振动的基本自振周期与土层分割数目的关系

本书在BC20试验的计算模型中，按土体的竖向划分，每层单元划分为一个土层，共划分了11个土层，符合图3-1对土层划分的要求。

采用等效线性模型的计算过程如下：计算时首先假定各层土一对动剪切模量Gd 1和阻尼比D 1，据此算出相应的有效动剪应变γd 1，由此γd 1值在土的动剪切模量和初始动剪切模量G 0之比与有效动剪应变γd关系曲线G d/G 0～γd、阻尼比与有效动剪应变关系曲线D～γd上找出对应的动剪切模量G d2和阻尼比D 2，重复以上步骤直至前后两轮的动剪切模量和阻尼比相差在允许范围内为止。由于剪应变随时间而变化，若按其最大值查找对应的G、β值，是不合适的，因而计算时取0.65倍的最大动剪应变作为有效动剪应变。

在假定每层初始剪切模量时，考虑了有效围压对初始剪切模量的影响。文献[129]给出下列关系：土的初始剪切模量与有效固结压力之比的平方根成正比，即：

式中：Gi——第i层土的初始剪切模量(Pa)；

Gi+1——第i+1层土的初始剪切模量(Pa)；

——第i层土的有效固结压力(Pa)；

3——第i+1层土的有效固结压力(Pa)。

本书利用ANSYS的参数化设计语言将上述土的等效线性模型及其计算过程并入ANSYS程序中，实现了材料非线性模拟。

试验中采用的重塑土进行了共振柱、循环三轴联合试验，得到了动剪应变γd在10-6～10-2范围内的Gd/G 0～γd、D～γd曲线[129]，如图2-5～图2-7所示。

4.土体与结构接触界面上的状态非线性模拟

混凝土材料和土这样两种材性相差很远的介质的界面，当应力水平超过一定限制时，位移的连续性就会受到破坏，会发生相对滑移和分离。在一定的荷载条件下，界面又会在张开后重新闭合。以前的研究者主要采用Goodman单元、薄层单元以及薄层土单元这些界面单元来模拟土体与结构(或基础)接触界面上的状态非线性。本书利用ANSYS程序的接触单元实现接触分析。交界面处的土表面作为接触面、结构(或基础)表面刚度相对土体要大，将其作为目标面，在接触面上形成接触单元、目标面上形成目标单元，然后通过相同的实常数将对应的接触单元和目标单元定义为一个接触对，并假定接触面上存在库仑摩擦。通过选择合理的参数，可实现土与结构界面上的黏结、滑移、脱离、再闭合的状态模拟。

在ANSYS程序中，考虑一个由N个物体组成的接触系统，根据虚功原理得到平衡方程如下：

等式的左边项表示内应力σ对虚应变δε所作虚功；等式右边第一、二项分别表示外力f B(面力)和f S(体力)对虚位移δu所作虚功；等式右边最后一项为接触力f C对虚位移δu所作虚功。

在式(3-3)的平衡方程中，必须满足相应的接触条件如下：

(1)法向条件

接触对法方向必须满足的条件是：

式中，g表示接触对之间的最短距离；λ表示接触对之间牵引力的法向分量。gλ=0表示如果g＞0则必须有λ=0成立，反之λ＞0则必须有g=0成立。无论接触对之间是否存在摩擦，此条件都必须满足。

(2)切向条件

本书在求解接触问题时，假设接触面上存在着库仑摩擦，摩擦系数为μ。令无量纲变量τ为其中T为接触对之间牵引力的切向分量；μλ表示接触对之间的库仑摩擦力。令u ·为接触对之间的切向相对速度。由库仑摩擦定理得到接触面上必须满足的条件是：(www.xing528.com)

令w为g和λ的函数，且方程w(g，λ)=0的解满足式(3-4)给出的条件；令v为τ和u ·的函数，且方程v(，τ)=0的解满足式(3-5)给出的条件。因此，接触条件可以表示为：

以上条件可以通过罚函数法或拉格朗日乘子法引入到由虚功原理得到的平衡方程中。与罚函数方法相比，拉格朗日方法不易引起病态条件，对接触刚度的灵敏度较小。可将变量λ和τ取为拉格朗日乘子，δλ和δτ分别为变量λ和τ的变分。δλ乘式(3-6)与δτ乘式(3-7)之和在物体I和J的接触面SIJ上积分，可得到约束方程为

对连续力学方程式(3-3)和式(3-8)进行离散化，得到以有限元网格节点的位移为基本未知量的方程组，对其叠代求解就得到了接触问题的答案。以上平衡方程以及接触条件都仅仅考虑了静力接触状态。在动力分析中，分布体力还应包括惯性力，而且在任何时刻，运动接触条件还必须满足接触体之间的位移、速度和加速度的协调性。

对于不考虑摩擦的接触问题，刚度矩阵是对称的，而考虑摩擦后会导致刚度矩阵不对称，叠代求解时采用不对称的求解器比使用对称的求解器花费机时更多，在ANSYS程序可以采用将不对称矩阵转化为对称矩阵的对称化算法。本文在求解接触问题时，假定接触面上存在库仑摩擦，采用了此对称化算法。

5.阻尼模型的选取

结构的阻尼主要与材料的特性有关，在结构与地基的相互作用问题中，地基的阻尼往往大于结构本身的阻尼，应分别输入各自材料的阻尼，按直接集成法，组成阻尼矩阵[130]。而ANSYS程序中采用的瑞利阻尼，即：

式(3-9)只能算出针对整个体系的质量阻尼系数α和刚度阻尼系数β，不能考虑地基与结构阻尼的不同。而在大多数结构问题中，α阻尼(质量阻尼)可忽略，在这种情况下：

ANSYS程序按上述只考虑β阻尼的思路提供了材料阻尼的输入方法，可针对不同材料输入相应的阻尼。总的阻尼矩阵按下式集成：

式中：

βj——每种材料的刚度阻尼系数；

[K j]——每种材料的刚度矩阵。

按上述输入材料阻尼后集成阻尼矩阵的方法，解决了结构和地基土阻尼不同的问题。本书计算中结构的阻尼比取5%，土体的阻尼比利用前述试验得到的D～γd曲线进行迭代。

6.重力的考虑

抗震设计及研究中，通常分别计算静力荷载及动力荷载作用下结构的内力，然后将它们线性叠加作为总的内力。本书采用的计算模型在土与结构(或基础)间加入接触单元模拟界面上的黏结、滑移、脱离、再闭合的状态变化。若在动力计算时不考虑重力的影响，重力引起的初始应力对接触状态的影响得不到考虑，也就不能真实反映土与结构的受力和变形情况。为此，本书将重力作为一种动力荷载并入动力计算。在计算用的地震波尚未施加前，将重力作为竖向的加速度场施加到体系上进行瞬态分析，作用一段时间后体系的反应趋于稳定，此时的反应为静平衡位置时的反应；然后，再将水平方向的地震波和重力同时作用到体系上，此时得出的反应是总反应，将总反应减掉静平衡位置的反应就可得到通常意义上的动力计算结果。

图3-2是模型比为1∶20的BC20试验模型体系、输入El Centro波激励时，考虑重力与不考虑重力情况的基底中点接触压力时程。从图中可看出，不考虑重力作用时的基底接触压力明显小于考虑重力作用时的基底接触压力，而且不考虑重力时基础底面发生了土与基础接触面的脱离现象，考虑重力作用时则没有发生该现象。可见，在考虑土与基础之间的状态非线性时，若不考虑重力的作用，将造成较大的误差，不能反映真实情况。

图3-2　基底中点接触压力时程(BC20试验模型)

7.结构中钢筋的处理

图3-3　1∶20柱截面配筋图

钢筋混凝土结构有限元模型一般有整体式、组合式和分离式三种[131]。其中分离式模型相对比较复杂，主要用于分析结构构件内微观受力机理，本书没有采用这种模型。这里以基础固定的钢筋混凝土单柱、受X方向El Centro地震波激励的情况，来说明建模时结构中钢筋的不同处理方法对计算结果和计算时间的影响，柱截面尺寸及其配筋如图3-3所示。分三种情况对其进行对比讨论：①忽略钢筋作用。将柱子材料简单地当作各向同性均一的素混凝土材料处理。②整体式模型。将钢筋弥散于整个单元中，并把单元视为连续均匀材料，钢筋对整个结构的贡献，通过刚度EI等效原则提高材料的弹性模量来实现，混凝土的初始弹模是22.058×103 MPa，按刚度EI等效原则调整后的弹模是38.584×103 MPa。③组合式模型。采用ANSYS软件中的钢筋混凝土组合单元SOLID65划分钢筋混凝土单柱。以上三种情况下的计算时间分别为0.5小时、0.5小时和11小时，结构的基频分别为28.12、37.10和35.70，通过进一步比较加速度时程图，发现忽略钢筋作用时的结果与组合式模型的结果相差较远，而整体式模型的结果与组合式模型的结果较接近。可见，采用刚度EI等效方法调整弹模的整体式模型，不仅可以得到令人满意的结果，而且建模简单、计算时间大为节省。

8.对称性的利用

为了降低自由度、减少计算时间，建模中可利用对称性原理。

在本试验中，土体、基础、上部结构组成的体系，几何形状关于X轴对称；地震波沿振动台X方向单向输入，可认为动荷载关于X轴正对称，因而整个结构沿X轴正对称。可沿Y=0平面将体系一分为二，取体系的二分之一为研究对象，在Y=0处加上正对称边界，约束边界上节点沿Y轴方向的位移、沿X轴和Z轴方向的转动。

试验的整个体系，几何形状关于Y轴也是对称的；地震波沿振动台X方向单向输入，可认为动荷载关于Y轴反对称，因而整个体系沿Y轴是反对称的。在不考虑大变形且进行线性计算时反对称是成立的，主要原因在于：考虑大变形时，结构的刚度矩阵要随结构的变形不断的变化，而发生大的变形以后，结构的刚度矩阵不再满足对称性，因而反对称不成立。考虑材料非线性时，由于在反对称荷载作用下，结构一侧受拉，一侧受压，且往往材料拉压不同性，造成结构刚度矩阵不满足对称性，因而反对称不成立。考虑接触面的状态非线性时，由于受拉和受压侧的状态比一致，即脱离、滑移等情况不一致，不满足几何形状的对称性，因而反对称不成立。

可见在不考虑大变形、材料非线性和接触面的状态非线性的情况下进行有限元计算，可以取体系的四分之一为研究对象，在Y=0处加上正对称边界，约束边界上节点沿Y轴方向的位移、沿X轴和Z轴方向的转动；在X=0处加上反对称边界，约束边界上节点沿X轴和Z轴方向的位移、沿Y轴方向的转动。

利用对称性，可以大大减少单元数目，降低自由度数，在求解时间的节约上是非常可观的。对于地震动输入相同的同一体系，采用相同网格划分。以整体、体系的二分之一和体系的四分之一为计算对象时，自由度数目分别为17228、8084和3950，所需计算时间则分别为35小时、14小时和6小时。经过计算验证，考虑与不考虑对称性，结果完全一致。

本书在进行计算时，由于采用了等效线性模型，因此计算中仅考虑线性或考虑材料非线性但不考虑接触面的状态非线性的情况采用了四分之一的模型，而在计算时同时考虑材料非线性和接触面状态非线性的情况则采用二分之一模型。

9.自由度不协调的处理

在计算模型采用的单元中MASS21、SHELL63单元每个节点具有6个自由度：3个平动自由度和3个旋转自由度，而采用的SOLID45单元每个节点只具有3个平动自由度，这样在不同单元交界处就存在自由度不协调的问题。例如在柔性容器(采用SHELL63单元)与土体(采用SOLID45单元)的交界面处，同一节点既属于SHELL63单元，每个节点具有6个自由度，该节点又属于SOLID45单元，每个节点只具有3个自由度。问题同样也存在于MASS21和SOLID45单元之间。

对于MASS21和SOLID45、SHELL63和SOLID45之间的自由度不协调问题，有两种解决方案：一是把SOLID45单元改为SOLID73单元。SOLID73单元与SOLID45单元的区别在于每个节点具有6个自由度，因而不存在自由度的不协调问题，缺点在于采用SOLID73单元，模型的自由度数目增大很多。以BC20试验模型为例，采用SOLID45单元时模型的自由度数为7793，改用SOLID73单元后模型的自由度数增加为15928，显然采用SOLID73单元从计算时间上看是不经济的。第二种办法是采用放松自由度的做法。即在MASS21和SOLID45或SHELL63和SOLID45两种单元的交界面处，将公共节点的3个平动自由度耦合起来，即令两个单元的X、Y、Z三个方向的位移分别相等，同时将公共节点的3个转动自由度放松，不加任何约束。由于转动自由度只对应结构动能的极小部分，这样处理对结果的影响很小。为了检验效果，计算了一个采用SOLID73单元的模型和一个采用放松自由度方法的模型，两种方法的计算结果非常接近，图3-4为采用SOLID73单元和放松自由度方法的A2、S6点的结果比较。采用放松自由度的方法可以合理地解决MASS21和SOLID45、SHELL63和SOLID45单元之间自由度不协调的问题。

图3-4　采用SOLID73单元和放松自由度方法的A2、S6点的结果比较(BC20试验模型)