英语原文共 19 页
适用于地下爆炸的地下结构的全耦合数值分析方法
摘要：
土中爆炸的物理过程和埋地结构的后续反应是极其复杂的。将所有这些过程组合成一个单一的分析模型涉及到一些数值上的困难，但是这样的模型将使底层物理过程的再现更加真实。本文提出了一种全耦合数值分析方法，采用SPH(光滑粒子流体力学)方法对近场介质进行建模，以适应大变形，而采用常规有限元方法对中场和远场土壤介质及结构响应进行建模。采用作者建立的稳健三相土模型对土体进行建模。数值模型与经验预测结果吻合较好。
关键词：埋在地下的结构;地下爆炸;SPH-FEM耦合的方法;应力波;结构响应;结构内冲击
1、介绍
地下钢筋混凝土结构用于防止常规武器影响的基本设施。通常这种结构呈盒状，部分或全部埋置。控制地下结构响应的物理过程非常复杂，包括炸药、土壤和地下结构之间的动力相互作用。主要现象包括火山口的形成或爆炸所造成的伪装术;激波和弹塑性波在土壤中的传播;和土与结构的相互作用。土体和钢筋混凝土的非线性特性和大变形使得整个物理过程在材料和几何非线性方面都具有高度的非线性。因此，为了充分描述整个过程，有必要采用数值方法。
通常采用两种数值方法来分析地下结构在爆炸荷载作用下的响应，即非耦合法和“耦合法”。在“非耦合法”中，主要物理过程分为几个连续的阶段;一个阶段的输出是下一个阶段的输入。在这方面，所审议的问题可分为三个阶段:(1)电荷爆炸和形成弹坑或伪装术;(2)冲击波传播;(3)结构的响应。耦合方法可以分为两类，即部分耦合法和完全耦合法。在部分耦合法中，将上述三个阶段分解为两个阶段，合并前两个阶段或后两个阶段。完全耦合方法将所有三个阶段组合在一个模型中。
关于爆炸荷载作用下地下结构的数值分析，已有许多研究报道[1-6]。这些研究大多基于非耦合方法或部分耦合方法。在这些方法中，一个基本的问题在于定义结构负载的充分性。在非耦合方法中，首先计算自由场中的应力或速度的历史。然后将这些时间历程作为分析结构响应的边界条件应用到结构上。因此，不能以现实的方式来考虑土与结构之间的相互作用。Hamdan和Dowling[7]指出，如果界面发生共振，非耦合方法可能导致不安全的结构。Henrych[8]也认为耦合效应可能是显著的，特别是当结构处于高密度介质(如水、土壤)中时。此外，界面效应，如滑移、分离和再结合，也是影响结构响应的重要因素。为了考虑这些影响，出现了几种耦合分析技术。Nelson[9]采用土岛法分析了埋地结构的墙体，模拟了埋地结构前的一小部分土，并将应力输入土岛的自由表面。Stevens和Krauthammer等[1,10]采用了一种将有限差分技术(FDT)与有限元法(FEM)相结合的混合方法，即采用适合于分析连续非线性介质中波浪传播的FDT对土壤进行建模，而采用有限元法对结构进行建模。这些耦合方法考虑了土体与结构之间的动力相互作用和耦合效应，但爆炸荷载仍需根据应力或压力时间历程进行定义。虽然在相对简单和对称的情况下定义爆炸荷载是合适的，但在结构形状不规则或地表效应显著的情况下(如浅埋结构)，爆炸荷载的定义就有问题了。在这些情况下，需要一个包括爆炸源在内的完全耦合方法。
现有文献中很少有研究将爆炸源纳入其中。除了计算成本高的问题外，模拟爆炸与土壤介质之间的动态相互作用的困难也是一个主要障碍。由于电荷近场的应力条件变化很大，因此很难对该区域的土壤行为进行建模，特别是考虑到土壤介质的多相性。
本研究旨在建立一种完全耦合的数值方法来分析地下结构在爆炸荷载作用下的响应。该方法将爆炸源、应力波在土体中的传播以及土体与结构之间的相互作用统一为一个模型。采用最先进的流体动力学数值技术和材料模型。方面的数值技术,光滑粒子流体动力学(SPH)方法与拉格朗日fi -合并夜间元素方法(FEM),即使用SPH模型近场响应和使用有限元模型的中间和远场地面运动和结构响应。在材料建模方面，采用作者开发的用于冲击荷载的鲁棒三相土模型[11,12]对土体进行建模。对于埋地结构，采用Riedel - Thoma-Hiermaier (RHT)混凝土模型[13]对混凝土进行建模，钢材采用弹塑性硬化模型。采用JWL (Jones - Wilkins-Lee)状态方程模拟了炸药的爆炸过程。
最后给出了一个数值算例，验证了该方法的有效性。数值结果与一些经验和工程观测结果进行了对比验证。
2、SPH-FEM耦合分析的基本思路
将爆炸源纳入结构响应的数值分析超出了一般结构分析规范的能力，因为这些规范通常不考虑能量守恒的考虑。该方法适用于模拟爆炸和冲击波在土体中的传播、土体-结构相互作用和结构响应等复杂过程。爆炸产物迅速膨胀。电荷附近的土壤发生了很大的变形。如果靠近电荷，结构也会发生大变形。另一方面，结构的响应依赖于接口和边界以及与之相关的效果。
基于材料颗粒与网格之间的相对运动，描述连续介质的方法主要有两种;一个是欧拉描述，另一个是拉格朗日描述。在欧拉描述中，网格固定在空间中，不同的物质粒子通过网格移动。 在拉格朗日描述中，网格与材料粒子重合。欧拉描述适用于网格可能高度扭曲的情况;但是用欧拉方法对材料的滑移、接触面等边界条件进行建模是非常困难的。拉格朗日描述更适用于变形不大但界面和自由边界影响较大的情况。为了利用这两种描述，提出了任意的拉格朗日欧拉描述(ALE)，以便分析人员可以选择网格是遵循材料(即拉格朗日)还是固定(即欧拉)。这种方法涉及到一种复杂的重分区技术，重分区过程通常需要有经验的用户[14]的干预。此外，它并不能完全避免网格畸变严重和时间步长随之急剧减小的问题，从而降低了计算效率。
欧拉法和拉格朗日法结合的困难主要是由于网格问题。一旦网格生成，表示物理区域的元素或网格就不能轻易更改。对于拉格朗日网格，大变形会导致严重的网格畸变，从而大大降低精度和时间步长。在欧拉网格中可以出现混合细胞，从而两种或两种以上的材料聚集在一起。混合网格会模糊材料间的界面和边界，因此欧拉网格很难与拉格朗日网格耦合。
为了克服网格化带来的困难，提出了一些无网格化方法。最重要的无网格方法之一是SPH(光滑粒子流体动力学)方法。SPH方法的优点是不需要跟踪材料界面(在这个意义上，SPH可以看作是一种特殊的拉格朗日方法)，从而避免了欧拉法的上述困难。无论溶液中有多少湍流，计算都可以继续进行，可以处理较大的变形。由于SPH法中没有真实的网格，避免了拉格朗日法中网格畸变严重的困难。
目前，SPH技术已被纳入几个水文编码[15,16]。虽然大多数问题可以用SPH方法建模，但也存在一定的局限性。正如Attaway等人[17]所提到的，用光滑粒子建模薄壁结构是低效的，因为需要很多小粒子，时间步长会变得非常小。在这方面，SPH和拉格朗日方法的耦合似乎是一个有效的解决方案，这样一来，结构等问题的一部分 可以用固体或壳单元来建模，而其他部分可以用光滑粒子来模拟。这种耦合方法对于所考虑 的问题类型是非常有效的。
3、计算框架
3.1、守恒方程
在三维空间中，质量、动量和能量守恒方程表示为
密度
动能
能量
密度rho;,V是体积,下标0表示初始值,m是质量,是应力,是应变,u是空间速度,e是能量,f是质量力,标记rsquo; lsquo; rsquo;是时间的一阶导数,i和j范围从1到3。
应变用变形表示，
其中X为空间坐标，W = Wi(Xj)为位置向量。
一般来说，应力和应变可以分为两部分，流体静力分量和偏差分量。前者对应体积变形，后者对应剪切变形。
其中sij为偏应力，eij为偏应变。
应力应变本构关系也可分为强度模型和状态方程两部分，分别描述剪切变形和体积变形。边界条件为比位移或牵引力，
其中x为点的当前坐标，x为参考坐标，t为时间，n为外法线，g为具体位移函数，或为位移或牵引边界条件作用的表面。
3.2.光滑粒子流体动力学(SPH)方法
SPH是一种无网格拉格朗日技术，最早于1977年应用于天体物理学。该方法的主要优点是绕过了用数值网格计算空间导数的要求。这避免了拉格朗日分析中涉及大变形冲击和爆炸载荷事件时经常出现的网格缠结和变形等严重问题。虽然名称中包含了“水动力”一词，但实际上材料的强度是可以包含进去的。
在SPH方法中，材料是由固定质量的粒子来跟随其运动来表示的。与基于网格的方法(如拉格朗日法或欧拉法)不同，SPH采用一种基于随机分布插值点的核近似来计算导数，而不需要假设哪些点是相邻的。
粒子携带质量m、速度矢量v、位置矢量x等物理量，构成守恒方程的计算框架。在这种方法中，每个粒子I都与距离它一定距离(通常假设为2h)内的所有其他粒子J相互作用。距离h称为平滑长度。互动是加权函数W (x- ,h)叫做平滑(或内核)函数。利用这个原理，一个连续函数的值，或者它的导数，可以在任意粒子I上根据周围粒子J的已知值来估计，使用下面的核估计:
其中f是三维位置向量x的函数，dxrsquo;是体积。
图1(a)说明了核估计的概念。核近似的数学推导的全部细节可以在[20]中找到。式(9)的一个常用的对称公式是
其中nabla;W是xJ,m是质量,q是密度。函数f可以是计算中的任何变量，如密度、应力或应变等。注意，在SPH方程的推导过程中，没有假定插值点之间存在连接性或空间关系，从而避免了网格纠缠。另一个重点是SPH节点可以使用与FEM元素相同的本构模型。
3.3. SPH - FEM耦合方法
精确的SPH模拟需要大量粒子贯穿整个SPH区域。因此，如果要求高精度或需要一些特殊的几何形状，如薄壁等，那么运行时间就会成为一个问题。SPH与拉格朗日有限元的结合是解决这一问题的一个潜在的好方法。利用有限元方法可以对低变形区的材料进行建模。SPH区域的粒子大小也可以分级，从而减少了计算量。图1(b)给出了SPH粒子如何嵌入传统拉格朗日有限元网格的基本概念。SPH粒子与有限元元件耦合的方式有两种。当它们附着在有限元单元上时，SPH粒子与有限元单元结合在一起，那么来自其他SPH粒子以及有限元元素的力作用在粒子上，形成运动方程。如果SPH粒子和FEM元素没有附着在一起，它们会沿着FEM元素的表面滑动，在这种情况下，必须使用一种特殊的滑动界面算法[15]。在本研究中，由于SPH颗粒代表近场土壤介质，因此SPH颗粒与有限元单元结合;SPH网格和FEM网格之间的界面不是一个材料界面。
4. 材料模型
该问题涉及的材料有土体、结构中的混凝土和钢筋以及高能量荷。用于描述这些材料的模型如下。
4.1.三相土壤模型
针对土的静动力响应建模，提出了多种不同的方法，包括弹性模型、内慢性模型、含速率无关模型和速率相关模型的塑性模型、粘塑性模型、临界状态模型等。但考虑到土体爆炸及爆炸后冲击波的传播，土体应力的变化范围比一般土体动力学中通常遇到的变化范围要大得多。电荷附近的压力可以达到几个GPa (giga pascal)，并且随着距离电荷的增加而迅速衰减。由于土是由固体矿物颗粒、水和空气组成的多相混合物，不同相的变形机理和贡献随应力条件的突变而变化。因此，为了对土中爆炸事件进行建模，需要一个鲁棒土模型来满足整个荷载范围的要求;为此，有必要对变形机理进行现实的反思。不幸的是，现有的土壤模型似乎都不符合上述要求。为了填补这一空白，作者最近开发了一个三相土模型来模拟冲击波在土壤中的传播[11,12]。该模型源于[8]中引入的概念模型，将土壤视为不同大小和形状的固体颗粒的集合，形成骨架，其空隙被水和空气填充。当外界荷载作用于土体时，固体颗粒、水、空气以及固体颗粒形成的骨架在不同的规律下发生变形。图2给出了三相土模型的基本思想及其数学表示，其中A、B、C元素分别对应固体颗粒、水和空气的变形，D、E元素描述了固体颗粒之间键连接的摩擦力和阻力。固体粒子之间的键由一系列细丝表示。模型公式大致可分为两个主要部分;状态方程和强度模型。固相、气相和水相的体积比分别为a1、a2、a3。以下是三相土壤模型的概述。
4.2.状态方程(EOS)
为了满足连续性要求，多相系统的总体积变化必须等于各相的体积变化之和，即
其中V为土壤单元体积，V0为单元初始总体积，Vw为水体积，Vg和Vs分别为空气体积和土壤颗粒体积。将孔隙体积记为Vp, Vp = Vg Vw，因此V = Vs Vp
压力载荷引起各相的变形，以及固体颗粒之间的摩擦和固体颗粒之间键的变形。摩擦力和由键产生的力都作用在固相上。满足平衡就得到
p是总静水压力,ps压力施加在固相,pa的压力由固体颗粒之间的摩擦,pb是由水和气体的压力,或“孔隙压力”,pc是由固体颗粒之间的压力,和pe是土骨架的压力等于pa和pc的总和。
式(12)描述了静压作用下的体积变形，其中可以由状态或应力应变关系的独立方程得到。
水的状态方程为
其中pw、pw0分别为水的电流和初始压力;cw0为水的初始声速，qw0为水的初始密度，qw为电流密度;kw是一个常数。
对于固体颗粒，Lyakhov(在Henrych[8]中)建议使用类似的状态方程，下标w用s代替，
当压力波在土壤中传播时，气泡受到突然的压缩，因此，可以用多元气体的状态方程来模拟孔隙中的空气，
式中pg0为空气的初始压力;qw0为初始压力下的空气密度，qg为压力pg下的空气密度，kg为等熵指数。
在土体骨架中，固体颗粒之间的摩擦力pa依赖于颗粒之间的法向应力。一般认为，法向应力与土体骨架的变形成正比。因此,
<!--资料编号：[4325]
以上是毕业论文外文翻译，课题毕业论文、任务书、文献综述、开题报告、程序设计、图纸设计等资料可联系客服协助查找。
