ALE 方法在沉箱码头地震液化分析的应用

刘 顺，唐小微，栾一晓

(大连理工大学海岸和近海工程国家重点实验室，大连 116023)

重力式沉箱码头因其耐久性、方便建造而成为港口常见的滨水构筑物之一，在神户港(Japan)、阿尔戈斯托港(Greece)、斯里兰卡港 (SriLanka)等重要国际港口都有应用，我国青岛港、厦门港、花莲港(台湾)等也有使用重力式码头。港口往往滨海或滨河建造，地基软弱，为满足承载力和稳定性要求，需采取地基处理措施以支撑上部结构物和工作载荷。图1 为重力式码头的典型断面形式，重力式码头的基础和墙后通常使用砂土、碎石土等进行处理。

图1 重力式码头墙的典型断面Fig. 1 Typical section of gravity wharf

饱和砂土地基在地震发生时可能发生液化，进而引发结构破坏，因此饱和砂土液化及减防灾措施研究始终是岩土工程领域研究的一个重点问题[1-6]。在过去几十年，也有很多重力式码头墙遭遇地震液化破坏的记录。1995 年，日本兵库县地震，神户港码头的大量重力式沉箱结构遭到破坏，整个港口陷入瘫痪。此外，中国台湾Chi-Chi 地震、土耳其Kocaeli 地震都观测到重力式码头液化破坏的现象。破坏原因往往是基底置换砂和墙后回填土液化引起的码头墙沉降、向海侧位移、倾斜等。

Abu 等[7]研究了地震作用下重力式码头土和结构相互作用，并分析了置换区和回填区密砂土实度对重力式岸壁结构的影响，Alyami 等[8]探究了渗透系数对重力式码头结构孔隙水压力的影响以及回填土和地基土的相对密实度对岸壁残余变形的影响。王丽艳等[9]详细分析了液化地基中沉箱码头残余变形随各影响因素的变化规律，Tong 和Schaefer[10]研究了振冲加密优化设计对重力式码头墙液化灾害的减弱作用。

上述研究均是基于有限元方法进行的数值模拟，作为一种有效的工具，有限元方法在工程领域有重要的应用。但在一些变形比较大的情况下，可能由于网格畸变导致计算中断、精度丧失等。地震调查报告显示，兵库县地震中神户港区地面最大加速度达到0.6g，部分沉箱顶部水平残余位移达5 m，残余沉降达1 m～2 m，向海侧倾角3°～5°。在强震作用下，基于传统有限元方法的沉箱码头动力分析可能由于较大的变形引起分析中断或者结果失真，不能合理反映强震作用下沉箱码头的最终变形，难以对沉箱码头的抗震优化设计起到有效的参考作用。基于算子分裂技术和自主开发的水土耦合动力液化数值分析平台，笔者发展了一种ALE 有限元方法，并通过饱和土柱一维固结和离心机模型试验模拟进行了验证[11]，详细实现过程可参考文献[11]，这里将其用于强震作用下的沉箱码头地震液化分析。

1 饱和砂土控制方程

基于水土二相混合理论，以土骨架位移u和孔隙水压力p作为变量，Uzuoka[12]推导了u-p形式的饱和砂土的场方程式：

2 ALE 方法简介

基于拉格朗日描述的有限元方法可能由于网格畸变导致计算中断、精度丧失，而基于欧拉描述的有限元方法处理运动的物质边界和相互作用问题是十分困难的。结合两种方法的优点，Donea 等[15]发展了一种任意的拉格朗日欧拉方法(Arbitrary Lagrangian Eulerian Method) 用于解决流体力学问题，随后该方法被众多学者用于解决固体力学中的非线性问题。

图2 物质域、网格域、空间域之间的映射关系Fig. 2 Mapping relationship among material domain, grid domain and spatial domain

根据ALE 方法的实现途径，可以将其分为两类：1) 完全耦合的(Fully Coupled)；2) 解耦合的(Decoupled)。完全耦合的ALE 方法将式(8)代入基于计算域的控制方程进行线性化，将物质运动和网格运动整合于离散化的有限元方程中，其基本形式如式(9)[16]所示：

式(10)中，耦合的ALE 方法需事先指定物质点和网格点的运动关系，这种关系通常比较简单，如超限映射法[17]等，对于结构形式复杂或变形较大的情况，往往不能作为有效的解决方案。在实践中，解耦合的ALE 方法更为普遍，该方法基于算子分离技术将一个完整的ALE 分析步分解为[18]：1) 忽略网格和物质相对运动的拉格朗日分析步，即更新拉格朗日方法；2) 将场变量在物质域和网格域之间传递的欧拉分析步。本文发展的ALE 有限元方法是在课题组开发的更新拉格朗日格式的水土耦合地震液化分析平台基础上，基于算子分离技术的解耦合ALE 有限元程序[11]。

3 有限元模型

3.1 模型简介

本文考虑了邵琪[19]在饱和砂土地震液化网格自适应的研究中采用的两种重力式码头地基处理形式。如图3 所示，沉箱位于海底黏土之上，沉箱底部用砂土置换，沉箱后部以砂土回填，两种方案不同之处在于沉箱地基的置换范围和形式。两种方案中，模型两侧边界设置宽高比较大(＞100)的单元和等位移边界(EDB)模拟自由场边界[11]，两侧边界和底边界为不透水边界，模型顶部(粗虚线)可自由排水。各区域的尺寸见文献[19]。

图3 沉箱模型及有限元网格Fig. 3 Caisson model and Finite element meshes

3.2 输入地震动

文中选择Kobe (Japan，1995，PA=0.834g)，Kocaeli (Turkey，1999，PA=0.349g)，Northridge(USA，1994，PA=0.568g) 三种地震波作为输入荷载，如图4 所示。图中，左侧为地震波，右侧为对应傅里叶谱，其中Kocaeli 波的主要频率范围低于1 Hz，Kobe 和Northridge 波的主要频率范围在1 Hz～10 Hz。考虑竖向地震动的影响，取水平向输入荷载的1/2 作为竖向地震动输入。

图4 地震波加速度及傅里叶谱Fig. 4 Acceleration and Fourier spectra of earthquake waves

图5 为近几十年来，记录的一些地震的峰值加速度分布，可以看出，强震动峰值加速度主要集中在0.5g左右。因此，本文对于强震动作用下沉箱码头的地震液化分析，所采用的三种地震波时程将按照其峰值强度与0.5g的比例进行调整到0.5g。

图5 近几十年地震最大加速度[22]Fig. 5 Peak acceleration of earthquakes in recent decades[22]

表1 模型参数Table 1 Model parameters

4 数值结果分析

4.1 网格形状

图6 和图7 显示峰值加速度为0.5g时，两种方案在三种地震荷载作用下，最终的网格变形，左侧为ALE 方法的网格，右侧为UL 方法的网格。

图6 方案1 网格变形图(Peak=0.5 g)Fig. 6 Mesh deformation of Case 1 (Peak=0.5 g)

图6 显示，方案1 中使用UL 方法的沉箱趾部前端在Kobe 波和Kocaeli 波作用下发生大范围的网格畸变，沉箱墙体后方的回填区在Northridge 波作用下存在局部的网格畸变；ALE 方法中网格在三种地震荷载作用下均处于健康的状态。

图7 显示，方案2 中在三种地震波作用下，UL 方法在沉箱趾部前端局部(Kobe)和较大范围(Northridge)以及墙体后部(Kocaeli)产生网格畸变，导致该方法失效；ALE 方法仅在沉箱趾部发生局部网格变形，整体网格处于健康状态。需要说明的是，在Kocaeli 波作用下ALE 方法失效是由于方案2 中沉箱发生较大的倾斜，置换区产生较大隆起变形，沉箱趾部侵入置换区，导致局部网格由凸变凹，同时也反映出此种处理方案下，沉箱在强震作用下的变形特点。

图7 方案2 网格变形图(Peak=0.5 g)Fig. 7 Mesh deformation of Case 2 (Peak=0.5 g)

综合图6 和图7，ALE 方法能够在强震作用下保持沉箱结构网格的整体质量，从而保持计算的持续性和准确性，反映沉箱结构的变形特点；UL 方法在强震作用下可能由于结构的大变形引起网格畸变导致数值方法失效。

4.2 残余位移

图8 为三种地震荷载作用下，方案1 沉箱顶点(图3(a)A点)在不同峰值加速度下的水平位移。可以看出，在峰值加速度PA=0.1g时，加载初期变形较小，两种方法的沉箱顶点位移基本一致，随着加载进行，UL 方法的顶点残余位移稍大于ALE 方法。强震作用时(PA＞0.3g)，两种方法对应的沉箱顶点位移随着加载进行差异明显，并且峰值加速度越大差异越大。图8(a)、图8(b)显示，在强震作用下，由于沉箱变形太大引起网格畸变导致UL 方法失效而不能完整反映沉箱位移变化，而ALE 方法则完整反映了沉箱顶点位移变化，并提供了最终的沉箱残余位移。图8(c)中，在Northridge 波作用下，UL 方法和ALE 方法均能反映沉箱顶点位移变化过程，但二者的差异随峰值加速度增大而变大。需要指出的是，强震时ALE 方法对应的沉箱顶点最终位移高于UL 方法中断时对应的位移。如表2 所示，在Kobe 和Kocaeli 波作用下，二者的位移差值最大达16%。

表2 方案1 沉箱顶点位移Table 2 Displacement of caisson apex of Case 1

图8 沉箱顶点水平位移时程(方案1)Fig. 8 Horizontal displacements of caisson apex of Case 1

图9 为三种地震荷载下，方案1 沉箱顶点在不同加速度峰值下ALE 方法的残余水平位移(图9(a))和竖向位移(图9(b))。从图中可以看出，沉箱的残余位移与地震荷载峰值加速度近似呈线性关系，与文献[9]中残余位移与输入荷载相关性一致。Abu 等[7]指出低频荷载对应的沉箱码头向海侧位移相对更大。可以发现，图9 中Kobe 波和Northridge 波对应的残余位移增长率相似，且数值相近，而Kocaeli 波作用下，残余位移随峰值加速度的增长率高于其他两种地震波，且强震时相同峰值加速度对应的残余位移远高于其他两中地震波对应的残余位移。考虑到输入荷载的频谱特性，沉箱码头的地震响应受低频地震波影响更大，这与文献[7]的结论一致。

图9 沉箱顶点残余变形(方案1)Fig. 9 Residual displacements of caisson apex of Case 1

图10 所示为三种地震荷载作用下，方案2 中沉箱顶点(图3(b)B点)在不同峰值加速度下的水平位移。在小震(PA=0.1g)和强震(PA＞0.3g)作用下，方案2 沉箱顶点位移均匀与方案1 的沉箱顶点位移发展趋势相似。不同之处在于，方案2在强震(PA＞0.3g)作用下，UL 方法失效时间比方案1 有所提前，并且ALE 方法在一些情况下(Kobe波PA=0.5g和Kocaeli 波PA=0.4g、0.5g)也发生失效。不同于UL 方法导致的大范围网格畸变，ALE 方法失效的共同特点是由于沉箱向海侧倾斜，侵入置换砂地基导致沉箱趾部某些网格失效(图7)，并未影响网格的整体质量。此种情况下，除更换整体网格或计算方法外并无有效解决途径。此外ALE 方法失效时，均已超过了输入地震波的主震范围，震动基本结束，且残余位移水平已经远超沉箱码头正常工作所允许的范围，此时的残余位移对于沉箱码头抗液化措施有效性评价仍具有参考意义。表3 为在峰值加速度为0.4g和0.5g的Kobe 和Northridge 波加载下，UL 方法和ALE 方法在计算中止时对应的沉箱顶点位移，可以看出在方案2 情况下，ALE 方法的位移远大于UL 方法，最大差值达56%。

图10 沉箱顶点水平位移(方案2)Fig. 10 Horizontal displacements of caisson apex of Case 2

表3 方案2 沉箱顶点位移Table 3 Displacement of caisson apex of Case 2

图11 为三种地震荷载作用下，方案2 中沉箱顶点在不同峰值加速度下使用ALE 方法的残余水平位移(图11(a))和竖向位移(图11(b))。与方案1中结果相似，沉箱残余位移与地震波峰值加速度近似呈线性关系，并且Kocaeli 波的残余位移大于Kobe 波和Northridge 波，反映出沉箱地震响应的低频敏感性。文献[7]中显示地基处理方案对将沉箱码头系统的模态频率有较大影响，具体表现为随着置换区地基区域的变大，模态频率逐渐变大；沉箱码头墙底部置换区加固对模态频率的影响高于沉箱码头墙后回填区。由图9 和图11 可以看出，在相同峰值加速度下，方案2 中Kobe 波和Northridge 波对应的残余位移大于方案1，这是由于地基处理方案的不同，改变了沉箱码头的自然频率，使得沉箱码头对Kobe 波和Northridge 波的地震响应变大，也减小了Kocaeli 波的地震响应与Kobe 波Northridge 波地震响应的差距。

图11 沉箱顶点残余变形(方案2)Fig. 11 Residual displacements of caisson apex of Case 2

4.3 相对误差分析

在地震液化的数值分析中，土体和结构由于液化引起的变形是关注的重点，因此，本文将单元应变作为误差评估对象来评估UL 方法和ALE方法在地震液化计算中对计算结果的影响。基于应变的平均相对误差定义形式可参考文献[19]。

图12 显示，方案1 中，在小震(PA=0.1g)作用时，除初期两种方法的相对误差相差不大以外，UL 方法平均相对误差随着变形增大而持续增大，ALE 方法对应的平均相对误差处于很低的水平。强震(PA＞0.3g)作用下，计算初期两种方法的平均相对误差亦差别不大。随着变形的增大，UL 方法对应的平均相对误差持续增大，在计算失效时，相对误差出现陡升，对应图6 中大范围的网格畸变，而ALE 方法则保持平均相对误差始终处于较低的水平，直至计算结束。

图13 显示，方案2 中，小震作用下，UL 方法和ALE 方法的相对误差发展趋势与方案1 相似。强震(PA＞0.3g)作用下，UL 方法计算失效，对应了图7 中部分网格畸变，ALE 方法在Kobe 波和Kocaeli 波作用下尽管比UL 方法计算持时较长，也未能完成整个计算过程。如第4.2 节所述，ALE 方法失效但并未影响整体网格质量，从图13可以看出UL 方法失效时，由于较多网格畸变，平均相对误差出现陡升，而ALE 方法虽然失效，平均相对误差始终处于很低的水平，且未出现由于大范围网格失效导致的误差陡增。

综合图12 和图13，ALE 方法用于沉箱码头的地震液化分析，在保证网格质量的基础上，提供了相比UL 方法精度更高的结果。

图12 平均相对误差曲线 (方案1)Fig. 12 Average relative error curves of Case 1

图13 平均相对误差曲线 (方案2)Fig. 13 Average relative error curves of Case 2

4.4 超孔隙水压力

为表征饱和砂土在地震作用下的孔压水平，这里定义超孔隙水压力比(excess pore water pressure ratio)为EPWPR=Δu/σ′y0，其中， Δu为单元超孔隙水压力， σ′y0为单元当前竖向有效应力。

图14 为方案1 在不同峰值的Kocaeli 波(0.3g～0.5g)作用下，ALE 方法中E1(图3(a))的超孔压比曲线和云图。从图14(a) 可以看出在不同峰值作用下，E1 的超孔压水平都较低，未达到液化水平。整体来看，孔压发展在该地震波激励下波动较为剧烈，0.5g对应的孔压曲线加载后期与0.3g和0.4g相差较大，达到较高的水平。从图14 超孔压云图可以看出，随着峰值加速度的增大，沉箱底部虽未达到液化，但高超孔压区逐渐扩展。回填区液化主要发生在距离墙后一定距离的区域，紧邻码头墙后方区域未发生液化，这种现象与实际地震中墙后液化现象一致[23]。随着峰值加速度增大，回填区和置换区的高孔压区逐渐产生连通，这将加大沉箱的残余变形，因此，设计中可考虑设置阻断置换区和回填区孔压转移的措施，如振冲加密置换回填区土体[10]或设置水平止水层。

图15 为方案2 在不同峰值的Northridge 波(0.3g～0.5g)作用下，ALE 方法中E1(图3(b))的超孔压比曲线和云图。方案2 中E1 孔压随加载进行逐渐累积，且峰值加速度增大，相应超孔隙水压力越大。图15 云图显示在方案2 的沉箱下部，有较大范围的液化区，随着峰值加速度的增大，液化区范围逐渐增大，并且有与图14 类似置换区和回填区高孔压区贯通扩展的现象。方案1 沉箱底部未出现明显液化区，而方案2 沉箱底部存在较大范围的液化区，说明方案2 的形式在沉箱底部易产生孔压积聚。

图15 超孔隙水压力曲线和云图(方案2)Fig. 15 Curves and contours of EPWPR (Case 2)

5 结论

基于饱和砂土控制方程和算子分裂技术，发展了一种解耦和的ALE 有限元方法，并将其用于饱和砂土地基沉箱码头地震液化分析。使用UL 方法和ALE 方法对比了两种沉箱地基处理方案在三种不同地震波激励下的地震响应，主要结论如下：

(1) 小震作用下，UL 方法和ALE 方法均能完成分析，并且UL 的计算结果略大于ALE 计算结果；强震作用下，UL 方法失效引起计算中止时的沉箱顶点位移与ALE 方法完成计算时相差较大。因此，文中发展的ALE 方法可用于对于UL方法不能够提供完整的沉箱码头地震响应的情况。此外，相对误差分析结果显示，无论UL 方法是否能完成计算，ALE 方法的平均相对误差均低于UL 方法，表明ALE 方法能够提供精度可靠的数值解。

(2) 沉箱码头的残余位移与地震加速度峰值近似呈线性关系，且对低频的荷载较为敏感，相应的地震响应也较大。地基处理方案不同将改变结构的自然频率，进而影响沉箱码头对地震荷载的敏感频率范围。因此，在沉箱码头的设计中，应考虑地基的处理形式对自然频率的影响，进而影响到沉箱码头的地震响应。

(3) 置换区和回填区在地震荷载作用下，可能存在高孔压连通区，并且地震峰值加速度越大，连通区范围越大。连通区将增大沉箱码头地震响应的残余位移，因此，在设计中，可考虑设置阻断置换区和回填区孔压转移的措施。