应力场数值模拟方法

Question

如题所述

Answer 1

近30年来，人们采用现场测试、实验室试验、理论分析与模型试验等多种方法，使岩土力学研究取得很大进展^{［162～166］}。如今随着计算机技术的快速发展，岩土力学的研究进入了一个新的阶段，其中数值计算方法已成为解决岩土力学问题的重要手段之一。

6.1.1 概述

许多工程分析问题，如固体力学中的位移场和应力场分布分析、电磁学中的电磁场分析、振动特性分析、传热学中的温度场分析以及流体力学中的流场分布等，都可以通过在给定边界条件下对其控制方程进行求解得到，但是利用解析方法只能求出一些方程性质比较简单且几何边界相当规则的极少数问题。对于大多数实际工程技术问题，由于物体的几何形状比较复杂或者问题的某些特性是非线性的，因而一般无解析解。为了解决此类问题，一般采用两种处理方法：一种是进行简化处理，将方程和边界条件简化为能够处理的问题，从而得到在简化情况下的解，但这种方法应用非常有限，且假设过多将会导致错误的解；另一种是在广泛接收现代数学和力学理论的基础上，借助于计算机和计算软件来获得工程上要求的数值解，这就是目前应用非常广泛的数值模拟方法。

目前在工程技术领域内常用的数值分析方法包括：有限单元法、边界元法、离散单元法以及有限差分法。最初常用的是有限差分法，它可以处理一些相当复杂的问题。但对于几何形状复杂的边界条件，其解的精度受到影响。20世纪60年代出现并得到广泛应用的有限单元法，使经典力学解析方法难以解决的工程力学问题都可以用有限元方法求解。它将连续的求解域离散为一组有限个单元的组合体，解析地模拟或逼近求解区域。由于单元能按各种不同的联结方式组合在一起，且单元本身又可有不同的几何形状，所以能适应几何形状复杂的求解域。但有限单元法需要的存贮容量常非常巨大，甚至大得无法计算。由于相邻界面上只能位移协调，对于奇异性问题（应力出现间断）的处理比较麻烦，这是有限单元法的不足之处。70年代末期，出现了另一种重要的数值方法为边界元法。边界元方法是把求解区域的边界剖分为若干个单元，将求解简化为求单元结点上的函数值，通过求解一组线性代数方程实现求解积分方程。上述两种数值方法的主要区别在于，边界元法是“边界”方法，而有限元法是“区域”方法，它们都是针对连续介质，只能获得某一荷载或边界条件下的稳定解。对于具有明显塑性应变软化特性和剪切膨胀特性的岩体，无法对其大变形过程中所表现出来的几何非线性和物理非线性进行模拟，这就使得人们去寻求适合模拟节理岩体运动变形特性的有效数值方法。

1971年Cundall，P.A^［167］提出了一种不连续介质数值分析模型——离散单元法。该方法优点在于适用于模拟节理系统或离散颗粒组合体在准静态或动态条件下的变形过程。离散单元法的基本原理不同于基于最小总势能变分原理的有限单元法，也不同于基于Betti互等定理的边界单元法，而是建立在牛顿第二运动定律基础上。最初的离散元法是基于刚性体的假设，由于没有考虑岩块自身的变形，在模拟高应力状态或软弱、破碎岩体时，不能反映岩块自身变形的特征，使计算结果与实际情况产生较大出入。Maini，T.，Cundall，P.A.^{［168～169］}等人针对刚体单元没有考虑岩块自身变形的缺点，利用差分方法提出了考虑岩石自身变形的改进的离散单元法，编制了通用的离散元程序UDEC（Universal Discrete Element Code），将离散元推广到模拟岩体破碎和变形情况，推动了离散元的进一步发展。我国学者也相继开展这方面的研究，王泳嘉教授^［170］等将离散单元法应用于采矿工程方面的研究。

6.1.2 FLAC数值模拟方法

（1）概述

数值模拟技术通过计算机程序在工程中得到广泛的应用。一直到20世纪80年代初期，国际上较大型的面向工程的通用程序有：ANSYS、NASTRAN、FLAC、UNDEC、ASKS以及ADINA等程序。它们功能越来越完善，不仅包含多种条件下的有限元分析程序，而且带有功能强大的前、后处理程序。

连续介质快速拉格朗日差分法（Fast Lagrangian Analysis of Continua，简写FLAC）是近年来逐步成熟完善起来的一种新型数值分析方法。把拉格朗日法移植到固体力学中，即将所研究的区域划分为网格，节点相当于流体质点，然后按照时步用拉格朗日方法来研究网格节点的运动，这就是固体力学变形研究中的拉格朗日数值研究方法。

FLAC与基本离散元法相似，但它克服了离散元法的缺陷，吸取了有限元法适用于各种材料模型及边界条件的非规则区域连续问题解的优点。FLAC所采用的动态松弛法求解，不需要形成耗机时量较大的整体刚度矩阵，占用计算机内存少，利于在微机的工程问题。同时，FLAC还应用了节点位移连续的条件，可以对连续介质进行大变形分析。

（2）数学模型

显式有限差分法的基本方程主要包括：平衡方程、几何方程、物理方程和边界条件。在FLAC^3D2.0中采用的拉格朗日描述方程，一般规定介质中一点由向量分量x_i，u_i，v_i，dv_i/dt（i=1，2，3）来表征，其分别代表位置、位移、速度和加速度分量。

其基本原理和基本公式简单叙述如下：

空间导数的有限差分近似

三维FLAC方法中采用了混合离散方法，区域被划分为常应变六面体单元的集合体；而在计算过程中，又将每个六面体分为常应变四面体，变量均在四面体上进行计算，六面体单元的应力、应变取值为其四面体的体积加权平均。

如图6.1所示，所研究区域任一四面体，节点编号为1～4，规定与节点n相对的面为第n面，设定其内任一点的速度分量为v_i，则由高斯散度定理得

煤岩动力灾害力电耦合

式中：V——四面体体积，m³；S——四面体外表面，m²；n_j——外表面单位法向向量分量。

图6.1 四面体

对于常应变单元，n_j在每个面上为常量，因此通过上式积分可得

煤岩动力灾害力电耦合

式中上标f表示f面的变量值，对于为线性分布的速率分量，速度分量的平均值为

煤岩动力灾害力电耦合

式中上标l表示节点l的变量值。将（6.3）式代入（6.2）式可得

煤岩动力灾害力电耦合

经过变换可得节点速率计算公式：

煤岩动力灾害力电耦合

1）平衡方程（运动方程）

显式有限差分法采用的平衡方程就是人们熟知的牛顿第二运动定律，即

煤岩动力灾害力电耦合

式中：F_i——节点合力在i方向分力，N；m_i——节点质量，kg；a_i——节点加速度在i方向分量，m/s²。

作用于各个节点的合力：外力（集中力、均布力、重力等）和内力（单元变形引起的应力在单元节点上的分量）。节点质量是根据节点相邻单元的面积（体积）和密度，按照面积（体积）加权求出。

FLAC^3D以节点为计算对象，将力和质量均集中在节点上，然后通过运动方程在时域内进行求解。节点运动方程可以表示为如下形式：

煤岩动力灾害力电耦合

式中：（t）———t时刻l节点在i方向的不平衡力分量，可以由虚功原理导出；m^l———l节点的集中质量，在分析静态问题时，采用虚拟质量；而在分析动态问题时，则采用实际的集中质量。

将（6.7）式左端用中心差分来近似，则可得

煤岩动力灾害力电耦合

2）变形协调方程——几何方程

作为连续介质力学，变形体之间必须满足变形协调方程（几何方程），否则变形体就会出现分离或嵌入。变形协调方程反映了位移与应变间的关系，对于某一时步的单元应变增量可由下式确定：

煤岩动力灾害力电耦合

求出应变增量后，即可由本构方程得到应力增量，各时步的应力增量叠加即可得到总应力，在大变形时，还需根据本时步单元的转角对本时步前的总应力进行旋转修正，然后即可由虚功原理求出下一时步的节点不平衡力，进入下一时步的计算。

3）物理方程——本构关系

物理方程反映应力与应变之间的关系，在程序中通常被称为材料模式或材料模型。在FLAC^3D2.0中提供了10种基本材料模型，它们是：①Null；②Elastic，isotropic；③Elastic，transversely isotropic；④Druck-Prager plasticity；⑤Mohr-Coulomb plasticity；⑥Ubiquitous joint plasticity；⑦Strain-hardening／softening Mohr-Coulomb plasticity；⑧bilinear strain-hardening/softening ubiquitous-joint plasticity；⑨Modified Cam-clay plasticity 和⑩elastic，orthotropic。

本文进行应力场数值模拟时采用的是Mohr-Coulomb应变硬化软化破坏准则，在FLAC^3D2.0中，Mohr-Coulomb 模型的破坏准则以主应力σ₁，σ₂，σ₃来描述，相应的应变为三个主应变ε₁，ε₂，ε₃。根据Hooke定律，应力、应变增量具有如下表达形式：

煤岩动力灾害力电耦合

式中α₁，α₂为材料常数，可以由体积模量K和剪切模量G确定：

煤岩动力灾害力电耦合

不失一般性，令σ₁≥σ₂≥σ₃，摩尔—库仑准则为

其中：

煤岩动力灾害力电耦合

式中C，φ分别为煤岩的粘聚力和内摩擦角。

FLAC^3D2.0的Mohr-Coulomb 破坏准则如图6.2所示。

图6.2 FLAC^3D的Mohr-Coulomb 破坏准则

本著作中就是选用上述的Strain-hardening／softening Mohr-Coulomb plasticity模型，对单轴压缩煤岩以及矿山地下煤岩独巷掘进时围岩的变形破坏过程进行模拟。

4）阻尼力

对于静态问题，FLAC^3D2.0在式（6.7）的不平衡力中加入了非黏性阻尼，以使系统的振动逐渐衰减直至达到平衡状态（即不平衡力接近零），此时节点运动方程变为：

煤岩动力灾害力电耦合

式中阻尼力（t）由下式确定：

煤岩动力灾害力电耦合

上式中α为阻尼系数，其默认值为0.8；而：

煤岩动力灾害力电耦合

5）初始条件与边界条件

边界条件包括面积力、集中载荷等应力边界条件和位移边界条件。此外也可加载体力和初始应力。在编写程序代码时，一般所有的应力和节点速度初始化为零，然后指定初始化应力。集中载荷则加载在面节点上，位移边界条件则以运动方程形式施加到相应的边界节点上。

边界条件分为应力边界条件和位移边界条件，应力边界条件为：

煤岩动力灾害力电耦合

式中：F_i———作用于节点i上的力；——作用于边界上的应力；n_j———边界上的法线沿j方向的矢量大小；Δs———边界的长度。

若是位移边界条件，应将边界条件以运动方程的形式施加到相应的边界节点上。

FLAC^3D2.0^［171］与FLAC^2D3.3也是由美国Itasca Consulting Group Inc开发的三维显式有限差分法程序，它可以模拟岩土或其他材料的三维力学行为。FLAC^3D2.0的计算循环过程如图6.3所示。

图6.3 FLAC^3D2.0的计算循环

6.1.3 FLAC数值模拟方法在采矿工程中的应用^{［172～179］}

采矿过程中围岩活动规律及巷道围岩稳定性问题涉及岩体力学特性、围岩压力、支护围岩相互作用关系及巷道与工作面时空关系等一系列复杂力学问题。随着我国经济建设的高速发展，岩土工程稳定性分析问题日益突出，除采矿工程外，在水利、交通（铁道和公路）、高层建筑的地基等行业也都存在着大量的岩土力学数值计算分析问题。能否用计算机数值模拟分析采矿岩层控制问题和岩土工程问题已成为一个大学岩层控制技术和岩土力学学科水平高低的标志之一。

与ANSYS、ADINA相比，FLAC 和UDEC的最大特点是计算分析岩土工程中的物理不稳定问题，因而特别适用于岩土工程中几何和物理高度非线性问题的稳定性分析，如采场的采动影响规律，软岩巷道的大变形问题，采动后的地表沉陷，露天矿的边坡稳定，水坝的稳定性等问题。

从力学计算方法上讲其主要特点

1）可以直接计算非线性本构关系；

2）物理上的不稳定问题不会引起数值计算的不稳定；

3）开放式程序设计（FISH），用户可以根据需要自己设计程序；

4）既可以分析连续体问题（FLAC），也可以分析非连续体问题（UDEC）；

5）可以模拟分析很大的工程问题；

6）高度非线性问题不增加计算时间。

在采矿工程中，许多学者利用FLAC软件对采矿过程中围岩活动规律及巷道围岩稳定性问题涉及到岩体力学特性、围岩压力、支护围岩相互作用关系及巷道与工作面的时空关系等一系列复杂的力学问题进行了一系列的研究，取得了显著的效果。梅松华等以施工期监测结果为基础，在正交设计原理的基础上，选定反演参数与水平，采用二维显式差分法FLAC进行弹塑性位移反分析。朱建明等在分析FLAC有限差分程序的基础上，提出了变弹性模量方法模拟时间因素对巷道围岩稳定性影响的衰减曲线，为揭示巷道围岩变形机理和有效指导围岩支护提供了有效的分析方法。来兴平等探讨了岩石力学非线性计算软件FLAC^2D3.3在地下巷道离层破坏数值计算中的应用。康红普对回采巷道锚杆支护影响因素进行了FLAC分析，认为FLAC^2D3.3在分析几何非线性和大变形问题方面性能优越。

在煤岩动力灾害预测中，这些方法的优点

1）可以提前知道煤与瓦斯突出、冲击矿压等煤岩动力灾害防治的重点区域；

2）可以得到大范围内的空间信息；

3）可以提前预测预报煤岩动力灾害的危险性；

4）可以确定在采掘过程中，应力的分布状况和集中程度。

在煤岩动力灾害预测中，这些方法也具有以下缺点

1）对实际问题均进行了简化处理；

2）对于煤岩体的力学特性，如弹性模量、泊松比等力学参数，也进行了简化，没有考虑其局部非均质性和各向异性；

3）只能作为一种近似方法使用。

Answer 2

通过Mechanical APDL Product Launcher 19.0，对轴承座铸件进行浇注结束后，凝固冷却这段时间的温度场、应力场数值模拟，可以发现轴承座在凝固冷却成型阶段温度最高、应力最大、最容易失效的部位还是轴承座的底座。所以综上所述，应当降低底座的厚度，提高其散热性能，避免因温度过高而引起的缩孔缩松，同时增大底座的长度宽度，目的是为了使得螺栓孔的间距变大，从而减少应力集中。必要情况下，可以采取自然时效、热时效、振动时效等方式消除应力，但是自然时效周期长不建议使用，热时效成本太大并且能耗不低，所以可以采用振动时效设备，对铸造过程中的等效应力进行一定量的清除。本次研究旨在通过结果分析，来探明轴承座铸件在凝固冷却成型过程中温度场和应力场的一个变化的情况，并根据这种变化，对铸件本身的结构加以一定的改进，以此来提升轴承座铸件的表面质量、结构精度、工作效率。因此，可以依据本次研究的结果，在以后轴承座的铸造生产过程中，加以运用，以此来优化轴承座凝固冷却的效果，具有一定的意义。