土体非线性对土－结构动力相互作用的影响——以Millikan图书馆为例

梁建文，杨晶丽，韩 冰，朱 俊

（1.天津大学建筑工程学院，天津 300350；2.南京航空航天大学土木与机场工程系，江苏南京 210016；3.北京工业大学建筑工程学院，北京 100124）

引言

地震激励下，土体与结构之间存在土－结构动力相互作用（SSI）效应。文献［1］采用间接边界元法（IBEM）研究了SV波入射下，层状半空间中刚性基础－剪力墙结构体系的SSI，研究了场地动力特性对SSI的影响。文中在文献［1］基础上，对近场土体进行有限元法（FEM）模拟，远场土体仍然采用间接边界元法模拟，采用等效线性化方法考虑土体非线性，进一步研究土体非线性对SSI的影响。

等效线性化方法采用等效剪切模量和阻尼比随土层动剪应变的变化近似考虑土体的非线性动力特性，具有形式简单，计算效率高的优点，是目前用来估计场地非线性的主要方法之一。Esteban 等［2］采用等效线性方法分析了土层的非线性动力特性对理想SDOF 结构模型的动力响应的影响，研究表明了当输入地震波的主频率接近考虑土体非线性后衰减的场地基频值时，土体非线性会增加上部的结构损伤，而且相比于结构高度、基础尺寸等参数，结构质量对SSI的影响更明显。Massumi等［3］采用等效线性模型模拟有限元土域，分析了SSI对结构动力响应的影响，结果表明土体刚度较小时，SSI会使结构的挠度和层间位移角显著增加；Quoc等［4］利用粘性土的等效线性骨架曲线考虑土体非线性，研究了非线性土－桩－结构动力相互作用的重要影响因素。Liu等［5］通过ABAQUS有限元软件研究土－结构相互作用的重要影响，分别结合等效线性模型和弹塑性模型考虑土体非线性，结果表明2 类非线性数值模型的模拟结果均与振动台试验结果吻合。韩冰等［6］采用有限元－间接边界元耦合法，通过等效线性化法考虑土层的非线性，研究了场地中存在透镜体时对土－结构动力相互作用的影响。

文中采用有限元－间接边界元（FEM－IBEM）耦合法，以著名的Millikan图书馆结构为例，计算El Centro波和Taft波从基岩面垂直入射时，不同地震波幅值，不同土体非线性参数下，结构基底剪力响应和顶部层间位移角响应，研究土体非线性对SSI的影响，得出一些有参考价值的结论。

1 模型与方法

采用FEM－IBEM 耦合法建立模型。采用有限元模拟近场土体，采用间接边界元法（IBEM）模拟远场土体，自动满足无穷远辐射条件。模拟近场土体的有限元（FEM）子域要足够大从而保证IBEM 子域几乎不受FEM子域散射的影响，使得远场土体非线性可视作与自由场土体非线性等效［7］，这样，SSI体系的等效线性迭代计算只需在有限元子域内进行，而将自由场非线性迭代后得到的土层特性参数作为IBEM 子域的计算输入参数，避免重复建立半无限空间的整体动力刚度矩阵。在此前提下，确定较小的FEM 子域尺寸以提高计算效率。文中不考虑基础与土体交界处可能发生的滑移等接触问题。

1.1 模型

土－结构动力相互作用（SSI）体系的物理模型如图1所示。

图1 SSI模型Fig.1 SSI model

剪力墙高为H，宽为W，刚性基础半径为r，弹性基岩上覆盖土层厚度为D。剪力墙与基础之间以及基础与土层之间都完全固接。地震波从基岩面垂直入射。

土层的特征参数包括土层质量密度ρL，阻尼比ζL和泊松比νL，以及拉梅常数λL和μL，土层剪切波速和压缩波速分别为βL和αL。弹性基岩的特征参数包括基岩质量密度ρR，阻尼比ζR和泊松比νR，以及拉梅常数λR和μR，基岩的剪切波速和压缩波速分别为βR和αR。

剪力墙的特征参数包括沿纵向单位长度的质量Μb，泊松比νb和阻尼比ζb，以及剪力墙基频fb和等效剪切波速βb= 4fbH。基础沿纵向单位长度的质量为M0。

地震作用下基础的位移响应包括水平平动位移dx，竖向平动位移dz和关于圆心O的转角θ。剪力墙沿水平方向和竖向的绝对位移分别为dbx和dbz，结构顶部相对于底部（或基础）的相对位移，即剪力墙自身发生的弹性位移分别为drelx和drelz。

FEM－IBEM 耦合法模型如图2 所示。其中，FEM 子域由近场土体组成，子域宽度2*m1，深度m2；IBEM子域由远场土体和弹性基岩组成。将土层进一步划分为n个子层，子层厚度dj（j=1，2，…，n）。IBEM 子域与FEM子域的界面称为边界S1，土－结构界面称为边界S2。FEM子域需要全部划分有限单元，并将位于边界S1和S2上的单元称为耦合单元，其余为内部单元，而IBEM 子域和结构只需在各自边界S1和S2处离散边界单元即可。沿共同边界S1，需要使IBEM 子域的边界单元数与FEM 子域的耦合单元数相同，并且单元节点重合。沿共同边界S2，FEM子域的耦合单元数与结构在S2上的边界单元数相同，并且单元节点重叠。同时，边界S1，S2上各对应单元响应之间应满足以下条件：

图2 FEM－IBEM耦合法计算模型示意图Fig.2 FEM－IBEM model

式中，上标“－”和“+”分别代表边界单元和耦合单元。单元位移u和应力t分别在各自子域独立求解得到。另外，在边界S1和S2上引入虚拟均布荷载P1，P2和P3，可看作是分别施加在IBEM 子域，FEM 子域和结构上的外荷载。而通过式（1）～式（3）建立的联系最终可求解出荷载P1，P2和P3，这也是采用耦合法求解整个体系地震响应的核心步骤。

1.2 IBEM子域的求解方法

通常，在间接边界元法中，远场任意位置的响应可认为是该点自由场响应与散射场响应的叠加。首先，自由场响应u′（x，z）和t′（x，z）可通过直接刚度法计算得到：

式中，S为场地整体刚度矩阵，由每一子层的动力刚度矩阵Ej（j=1，2，…，n）和基岩刚度矩阵ER按一定顺序集整而成；外荷载仅在基岩面由入射波引起，因此，荷载向量那么，通过方程（4）可解出各土层界面上的位移当第j层土的上下界面位移已知后，即确定了（z=0，z=dl）处的边界条件，就可以在该土层的局部坐标系下求解土层中任意位置上的自由场位移响应u′(x,z)和应力响应t′(x,z)了。

求解散射场响应时，需要在边界S1上离散N个边界单元并引入格林影响函数为位移格林函数，其物理意义为在第j个边界单元上依次施加单位水平均布荷载，单位竖向均布荷载时第i层土中任意点的位移响应。其求解过程为：先固定荷载作用层上下界面，根据该土层波数域中的运动方程得到上下界面节点位移，应力的特解和齐解；再取消固定，将节点反力作为外荷载施加于荷载作用层的上下界面，由式（4）求得场地中各土层的界面位移，那么任意i土层内任意位置的位移响应ui(k′,z)也相应可解（k′为散射场波数），注意荷载作用层的响应还应增加上一步的方程齐解和特解部分；最后对ui(k′,z)进行傅里叶变换得到第i层土在空间域中的位移响应即为其位移格林函数。同样，场地的应力格林函数可由本构关系得到。

将格林函数与边界上虚拟均布荷载P1相乘即为散射场响应，其中单元荷载向量可认为是FEM子域施加在IBEM子域上的荷载(j= 1,2,…,N)，上标“H”、“V”分别代表水平方向和竖向的荷载。

间接边界元法的详细计算方法见文献［1］和［8］。

因此，根据间接边界元法，IBEM子域在共同边界S1上的边界单元的平均位移响应和平均应力响应如下所示：

1.3 FEM子域的求解方法

近场有限单元法子域内的土体需要按四节点等参单元划分网格。FEM 子域中沿S1边界划分的耦合单元数为N，沿S2边界划分的耦合单元数为M。根据有限单元法理论中平面应变问题的求解思路，集整所有土单元的单元刚度矩阵，建立FEM子域的整体动力刚度矩阵KF。KF可分区表示为：

因此地震作用下，由PF引起的FEM子域的动力响应为：

式中：A为将单元均布荷载转化为单元等效节点荷载的转换矩阵；uu为内部节点的位移向量；us为耦合边界节点的位移向量。

1.4 结构响应的求解方法

结构顶部的位移响应为［9］：

结构的基底剪力Fbx为：

1.5 耦合方程

由式（5），式（8）和式（9）给出的沿S1、S2的单元位移和应力响应，式（1）～式（3）可具体化为：

由式（16），可将P3写作用P2替代的关系式如下：

然后将式（17）代入式（14）和式（15），有

联立式（18）和式（19）并写成矩阵方程形式，即为FEM－IBEM耦合法的总耦合方程，如下：

由耦合方程式（20）最终可解得间接边界元法子域的虚拟荷载向量P1和有限单元法子域的虚拟荷载向量P2，再由式（17）可得到虚拟荷载向量P3。最后，通过式（11）～式（13）可求解结构的地震响应。

1.6 等效线性化迭代计算

等效线性化的计算流程如下。根据土层的初始剪切模量G0和阻尼比ζ进行首次运算，将土单元经过前一次完整动力计算后的整个地震持时内的最大剪应变γmax，乘以等效系数0.65得到单元等效剪应变γeff；再分别根据等效剪切模量比、阻尼比与等效应变的关系曲线G/G0－γ，ζ－γ得到土单元的新剪切模量值Gi和新阻尼比值ζi，并将Gi和ζi作为开始本次动力迭代计算的输入参数。若本次迭代计算后得到的单元等效剪应变与前一次的相差|γi-γi-1|/γi-1在某允许值σ范围内，即认为该单元收敛［10］，否则，继续求取新的Gi+1和ζi+1进入下一次动力计算。文中取σ=0.05，并以某次迭代有95%以上的单元收敛作为收敛条件。将最后一次动力迭代计算后得到的系统响应作为完成非线性计算的最终输出结果。

2 方法验证

图3为文中耦合方法和文献［1］的间接边界元法的计算结果对比。基础宽度r=10 m，土层厚度D=30 m，40 m，土层和基岩参数为ρR=ρL=1 900 kg/m3，ζL=0.05，ζR=0.02，νL=νR= 1 3，βL=200 m/s，βR=400 m/s，1 000 m/s。上部结构H=20 m，W=20 m，βb=200 m/s，ε=βLH βbr=2，ζb=0，νb= 1 3，MbM0=2，M0MS=1。以基岩面垂直入射SV 波时的基础位移和结构顶部的相对位移为例进行了对比验证。可以看出，文中计算结果与文献［1］的结果相吻合，证明了文中方法的有效性。

图3 垂直入射SV波时基础位移响应和顶部相对位移响应与文献［1］的结果对比Fig.3 Comparison of the foundation displacement and the relative displacement at the top of the structure for vertical incidence of SV-wave with the results in the paper［1］

图3 （续）Fig.3 （Continued）

3 算例分析

本节以著名的Millikan图书馆为例进行算例分析，具体参数见文献［11－12］。

图书馆结构参数H=44 m，W=25 m，r=12.5 m，fb=2.33 Hz，νb= 1 3，ζb=0.01，MbM0=7.5，M0MS=0.2。图书馆场地参数如表1 所示。沿弹性基岩面垂直入射的地震波分别选用幅值为0.1 g，0.2 g的El Centro 波和Taft 波，波型为SV 波。地震波的加速度时程曲线如图4（a）和（b）所示。对上述4 个工况分别进行线性和非线性分析。其中，非线性分析时各工况下又分开采用了2类非线性参数，一次计算过程中所有土层均使用同一组土动力参数。2组土层动剪切模量、阻尼比与剪应变的关系曲线如图5所示，取自文献［13］。

图4 基岩输入地震波的加速度时程曲线Fig.4 Acceleration time-history curves of the seismic waves

图5 两组等效线性G/G0－γ和ζ－γ关系曲线Fig.5 Two types of the relation curves of G/G0－γ and ζ－γ

表1 Millikan图书馆场地参数［11］Table 1 Parameters of the Millikan library building site［11］

为了方便描述考虑土体非线性后非线性SSI响应与线性SSI响应之间的差异，定义2个参数如下：

式中：χF表示考虑土体非线性后的结构基底剪力响应较线性响应的变化幅度；χD表示考虑土体非线性后的结构顶部层间位移角响应较线性响应的变化幅度（Db=drelx/H）。

图6 和图7 分别给出了Taft 波和El Centro 波入射下，结构基底剪力Fbx非线性响应与线性响应的比较。由图可见，在考虑土体非线性后，Fbx时程响应峰值均有不同程度的减少。当基岩输入幅值为0.1 g 的Taft 波时，采用两类非线性参数（Ⅰ和Ⅱ）进行等效线性计算后的基底剪力分别比线性情况时降低了|χF|Ⅰ=8.4%，|χF|Ⅱ=20.54%；当基岩输入幅值为0.2 g 的Taft 波 时，|χF|Ⅰ=24.18%，|χF|Ⅱ=43.08%。当基岩输入0.1 g 幅值的El Centro 波时，Fbx非线性响应的降低幅度分别为|χF|Ⅰ=22.13%，|χF|Ⅱ=31.75%；当基岩输入0.2 g 幅值的El Centro 波时，|χF|Ⅰ=44.53%，|χF|Ⅱ=54.28%。显然，采用非线性更强的Ⅱ类参数计算后的基底剪力下降幅度|χF|Ⅱ均大于Ⅰ类情况时的下降幅度|χF|Ⅰ；而且地震波幅值为0.2 g时的基底剪力下降幅度均比0.1 g时更明显。

图6 基岩输入Taft波情况结构基底剪力Fbx时程曲线Fig.6 Time-history curves of the base shear for Taft wave input from bedrock

图7 基岩输入El Centro波情况结构基底剪力Fbx时程曲线Fig.7 Time-history curves of the base shear for El Centro wave input from bedrock

同样，图8 和图9 分别给出了Taft波和El Centro 波入射下，结构顶部层间位移角Db的非线性响应与线性响应的比较。由图可见，在考虑土体非线性后，层间位移角Db的时程响应峰值均有不同程度的减少。当基岩输入幅值为0.1 g 的Taft 波时，层间位移角Db的非线性响应降低幅度分别为|χD|Ⅰ=6.65%，|χD|Ⅱ=15.71%；当基岩输入幅值为0.2 g 的Taft 波时，|χD|Ⅰ=17.22%，|χD|Ⅱ=34.89%。当基岩输入0.1 g 幅值的El Centro 波时，Db非线性响应的降低幅度分别为|χD|Ⅰ=14.32%，|χD|Ⅱ=24.95%；当基岩输入0.2 g 幅值的El Centro 波时，|χD|Ⅰ=40.93%，|χD|Ⅱ=53.31%。显然，对于结构的层间位移角响应，非线性响应规律与基底剪力类似，土层的非线性越强，地震波幅值越大，非线性响应越明显。

图8 基岩输入Taft波情况层间位移角时程曲线Fig.8 Time-history curves of the structure drift for Taft wave input from bedrock

图9 基岩输入El Centro波情况层间位移角时程曲线Fig.9 Time-history curves of the structure drift for El Centro wave input from bedrock

图10和图11分别给出了Taft波和El Centro波入射下，土体非线性（I和II）和线性情况下结构顶部绝对加速度Ab的傅里叶谱比较。由图可见，考虑土体非线性后，系统频率明显降低。当基岩输入幅值为0.1 g的Taft波时，SSI系统频率由2.03 Hz降到了1.54 Hz；基岩输入幅值为0.2 g的Taft波时，SSI系统频率由2.03 Hz降到了1.17 Hz。当基岩输入0.1 g幅值的El Centro波时，考虑土体非线性后，SSI系统频率由1.95 Hz降到了1.46 Hz；输入0.2 g幅值的El Centro波时，SSI系统频率由1.95 Hz降到了1.17 Hz。显然，地震波幅值越大，系统频率降低幅度越大。

图10 基岩输入Taft波时结构顶部绝对加速度时程曲线及其傅里叶谱Fig.10 Time-history curves of the absolute acceleration at the top of the structure and their Fourier spectra for Taft wave

图10 （续）Fig.10 （Continued）

图11 基岩输入El Centro波时结构顶部绝对加速度时程曲线及其傅里叶谱Fig.11 Time-history curves of the absolute acceleration at the top of the structure and their Fourier spectra for El Centro wave

4 结论

文中建立了层状半空间中SSI体系的二维FEM－IBEM 耦合法模型，采用等效线性化方法同时考虑近场和远场土体的非线性，以著名的Millikan图书馆为例，计算了Taft波和El Centro波垂直入射情况下，分别采用2 类非线性参数时的结构基底剪力和层间位移角。最后通过与线性SSI 响应比较，研究了土体非线性对SSI的影响。主要结论如下：

（1）土体非线性对SSI 有明显影响：考虑土体非线性后，结构基底剪力响应和层间位移角均比线性情况时有明显降低；结构系统频率均明显低于线性情况。文中算例中，基岩输入Taft波时基底剪力最小和最大降幅分别为8.4%和43%；层间位移角最小和最大降幅分别为6.7%和35%。基岩输入El Centro 波时基底剪力的最小和最大降幅分别为22%和54%；层间位移角最小和最大降幅分别为14.3%和53.3%。

（2）土体的非线性越强，对非线性SSI 的影响越明显。文中算例中，采用非线性更强的Ⅱ类土体参数时的非线性SSI 响应衰减幅度均比采用Ⅰ类参数时更明显。基岩输入幅值为0.1 g 的Taft 波时，采用两类非线性参数（Ⅰ和Ⅱ）的基底剪力的降幅分别为8.4%和20.5%，层间位移角的降幅分别为6.7%和15.7%；基岩输入幅值为0.2 g的Taft波时，基底剪力的降幅分别为24%和43%，层间位移角的降幅分别为17%和35%；基岩输入幅值为0.1g 的El Centro 波时，采用两类非线性参数（Ⅰ和Ⅱ）的基底剪力的降幅分别为22%和32%，层间位移角的降幅分别为14.3%和25%；基岩输入幅值为0.2 g的El Centro波时，基底剪力的降幅分别为44.5%和54%，层间位移角的降幅分别为41%和53.3%。

（3）地震波入射幅值对非线性SSI 有明显影响，且地震波幅值越大，非线性SSI 越明显。文中算例中，基岩输入幅值为0.2 g 的Taft 波或El Centro 波时，结构基底剪力，层间位移角以及系统频率的非线性衰减幅度均比幅值为0.1 g时更明显。基岩输入Taft波（0.1 g和0.2 g）时系统频率的降幅分别为24%和42%；基岩输入El Centro波（0.1 g和0.2 g）时系统频率的降幅分别为25%和40%。