第四章 FLAC3D数值模拟

第四章 FLAC3D数值模拟 4.1 FLAC3D数值分析软件介绍 4.2 模型建立与运行

4.2.1 建立模型 4.2.2 各工况的数值模拟 （1） （2） （3）

4.3 水平荷载下刚性单桩工作性状分析

4.4 水平荷载下带帽刚性桩工作性状分析

4.5 水平荷载下带帽刚性桩复合地基工作性状分析

4.5.1 桩长、桩径及长径比、桩帽大小、褥垫层厚度带帽刚性桩

复合地基应力场和位移场的影响

4.5.2 各计算参数的敏感性分析

注：本章与第三章要相对应，分三类（刚性单桩、带帽刚性桩、带帽

刚性桩复合地基），每一类又有多少工况，依据是什么？

第四章 FLAC3D数值模拟

4.1 FLAC3D数值分析软件介绍

自R.W.Clough 1965年首次将有限元引入土石坝的稳定性分析以来，数值模拟技术在岩土工程领域获得了巨大的进步，并且成功的解决列入许多重大的工程问题。近代个人电脑的出现以及其计算能力的飞速发展，使得分析人员在室内进行岩土工程数值模拟成为可能，也使得数值模拟技术逐渐成为岩土工程研究和设计的必不可少的方法之一。

数值模拟的优势在于有效的延伸和扩展了分析人员的认知范围，为分析人员洞悉岩土体内部的破坏机理提供了强有力的可视化工具。因此，岩土工程数值模拟软件必须做到专业性、可视化和完善的信息输出能力，才能更方便的帮助分析人员研究问题。FLAC3D等软件的出现是数值模拟工程发展的一个里程碑。

FLAC3D软件是由Itasca公司研发推出的一款数值分析软件，其界面简单明了，特点鲜明，使用特征和计算特征别具一格，因此在岩土工程中应用广泛，并享有盛誉。

FLAC3D是一个三维有限差分程序，它是二维有限差分程序FLAC2D的扩展，能够进行土质、岩石及其他材料的三维结构受力特性模拟和塑性流动分析。FLAC3D可对分析的单元进行线性或非线性本构模型的定义，当材料发生屈服流动后，网格能够相应的发生变形和移动（大变形模式）。其采用了显示拉格朗日算法和混合-离散分区技术，能够非常准确的模拟材料的塑性破坏和流动。相对与我们熟知的ansys得有限元程序，它不需要形成刚度矩阵，因此大大其对内存空间的需求，能轻易的求解大范围的三维问题。

FLAC3D的优点

1、对于模拟破坏和塑性流动采用的是“混合离散法”。该方法相比于有限元法中采用的“离散集成法”更加精确、合理。

2、即使模拟的系统是静态的，仍然采用动态运动方程，这使得FLAC3D在模拟屋里上的不稳定过程不存在数值上的障碍。

3、采用了一个“显式解”方案。因此显式解方案对非线性应力-应变关系的求解所花费的时间几乎与线性本构关系相同。而隐式解求解方案将花费较长的时间求解非线性问题。尤其是它不需储存较大的刚度矩阵，使得采用中等容量的内存就能求解多单元结构，模拟大变形问题几乎并不比小变形问题花费的计算时间多，因为没有任何刚度矩阵需要修改。

FLAC3D的缺点

毋庸置疑，FLAC3D软件是一个十分优秀的岩土工程数值模拟软件，它的实

用性和专业性得到了广泛的证实。但是不可否认，FLAC3D软件也存在诸多问题，主要如下：

1、FLAC3D求解工程问题受划分的网格的密度的影响很大。对于一般的弹性问题，它比一般的有限元程序运行得要慢，其求解的时间与网格单元的三分之四次方成正比。由此可见，FLAC3D对网格的多少十分敏感，对于同一个模型，不同密度的网格单元可以导致求解时间相差数倍。

2、FLAC3D在某些模式下的计算求解时间很长。由于很多物理过程如固结过程、长期动力影响、材料流变等都与时间密不可分，不得不考虑时间的影响。对于这一问题，FLAC3D采用的是真实时间，因而造成求解时间很长。

3、FLAC3D的前处理功能相对不够完善。对于复杂的三维模型的建立仍然非常的困难。尽管FLAC3D软件提供了12中初始单元模型，通过对这些单元进行组合可以方便快捷的建立规则的三维工程地质模型；同时也可以通过内置的FISH语言，通过编写命令来调整、完善特殊的计算模型，使之更加符合工程实际。但是由于FLAC3D运用的是命令流的形式进行人机交互，加上FISH语言需要一定的编程能力，对操作人员来说是一项重大的挑战。既是是有经验的分析人员，复杂模型的建立也是一项费时费力的苦差。这也是严重制约FLAC3D进一步推广的主要原因。

FLAC3D的应用范围

尽管最早开发FLAC3D软件是运用于岩土工程和采矿工程的力学分析，但由于该软件具有强大的解决复杂力学问题的功能，因此FLAC3D的应用范围已经拓展到土木建筑、地质、交通、水利等工程领域，成为这些专业领域进行工程分析和设计中不可缺少的工具。其研究的范围主要有：

1、岩、土体的渐近破坏和崩塌现象的研究。 2、岩体中断层结构的影响和加固系统的模拟研究。 3、岩、土体材料固结过程的模拟研究。 4、岩、土体材料流变现象的研究。

5、高放射性废料的地下储存效果的研究分析。

6、岩、土体材料的变形局部化剪切带的演化模拟研究。

7、岩、土体的动力稳定性分析、土与结构的相互作用分析及液化现象的研究等。

FLAC3D的求解流程

采用FLAC3D进行数值模拟时，有三个基本部分必须指定：有限差分网格、本构关系和材料特征、边界和初始条件。

网格是用来定义分析模型的几何形状，本构关系和与之对应的材料特性用来

表征模型在外力作用下的力学响应特性，边界和初始条件用来定义模型的初始状态。

在处理完这些条件后，便可通过运行软件求解得模型的初始状态；接着在初始状态的基础上进行开挖等其他模拟条件的加载，进而解得模型在模拟条件变化后作出的响应。

FLAC3D是拉格朗日连续介质法（Fast Lagrangian Anaysis of Continua）简称，拉格朗日连续介质法属于有限差分法，因此FLAC3D为有限差分软件，并非有限元软件。除此之外，FLAC3D采用的是混合离散法和动态松弛法，这也与有限元软件有着显著地不同。

有限差分法

运用数值计算的方法求解偏微分方程时，将任意处的导数用有限差分近似代替，从而把复杂的偏微分方程转化成简单的代数方程进行求解的方法就叫做有限差分法。其求解偏微分方程的步骤如下：

1、区域离散化。把待求的偏微分方程的求解区域细分成由有限个格点组成的网格，如图4-1。

2、近似替代。采用有限差分公式替代每一个格点处的导数。

3、逼近求解。对替换成的代数方程进行求解，其求得的精度取决于步骤1中区间的离散程度。

求解过程

1、导数的有限差分近似

FLAC3D采用的是混合离散法，因此其计算都是在四面体上进行的。现在以一个四面体说明计算时导数的有限差分近似过程。如图4-2所示，四面体的节点编号为1～4号，第n面表示与节点n相对的面，假设其内任意点的速率分量为vi，则由高斯公式可得：

?vVi,jdV??vniSjdS （4-1）

式中，V为四面体的体积，S为四面体的外表面，nj为外表面单位法向向量分量。

对于常应变单元，vi为线性分布，nj在每个面上位常量，由式（4-1）可得：

vi,j??13V4?vni?1li(l)jS(l) （4-2）

式中l表示节点l的变量，(l)表示面l的变量。 2 运动方程

FLAC3D软件中，体系中力和质量都集中在四面体的节点上，以节点作为计算对象，然后通过运动方程在时域内进行求解。节点的运动方程可表示为：

?vi?tl?Fi(t)mll （4-3）

式中，Fil(t)为早t时刻l节点的在i方向的不平衡力分量，可以通过虚功原理推导出；ml为i节点的集中质量，在分析静态问题时，采用虚拟质量以保证数值稳定，在分析动态问题时这采用实际的集中质量。

将（4-3）式左端用中心差分来近似，则可得到：

v(t?lt?t2)?v(t?lt?t2)?Fi(t)mll?t (4-4)

3 应变、应力及节点不平衡力

FLAC3D由速率来求某一时步的单元应变增量，如式（4-5）：

?eij?12(vi,j?vj,i)?t （4-5）

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

4 阻尼力

对于静态问题，FLAC3D在式（4-3）的不平衡力中加入了非粘性阻尼，以使系统的振动逐渐衰减直至达到平衡状态（不平衡力趋向于零）。此时式（4-3）则变为：

?vi?tl?Fi(t)?fi(t)mlll （4-6）

阻尼力fil(t)为：

fi(t)???Fi(t)sig(vi)lll （4-7）

式中?为阻尼系数（默认为0.8），而

sign(y)?1

y>0 y<0 (4-8) y=0

sign(y)??1 sign(y)?0

5 计算循环

由以上步骤可以看出FLAC3D的计算循环，将其用图形表示出来如下：

模拟工况简介

运动方程 对每个节点 1、由应力及外力利用虚功原理求节点不平衡力 2、由不平衡力求节点速度 本构方程 对每个单元 1、由节点速度求应变增量 2、由应变增量求应力增量及总 图4-2 FLAC3D计算循环图

对于单桩在水平荷载作用下的数值模拟，主要研究桩体尺寸及荷载对其受力响应的影响。即主要研究桩长、桩径、荷载作用等工况的不同对水平荷载作用下的单桩桩身位移和弯矩的影响。

模型为某混凝土桩，桩长分两组，分别为20m和10m，桩径分两组，分别为0.6m和0.8m，弹性模量为25Gpa，桩顶设置载荷板。桩周土体为均质沙土，土体用长宽高分别为8m、8m、28m、及8m、8m、18m的长方体来模拟，弹性模量为20MPa，泊松比为0.3，桩顶（即泥面处）受水平荷载的作用，荷载大小分4组，分别为9 kN、22.5 kN、45 kN、90kN。桩与土的接触位置设置接触面，土体的约束条件为：土体对底层约束竖向位移（z方向），土体侧面约束该面的法相位移。由此建立起几何模型如图4.3和4.4所示。

Itasca Consulting Group, Inc.Minneapolis, MN USAFLAC3D 3.00Step 22973 Model Perspective11:28:12 Wed Feb 15 2012Center: X: -2.665e-015 Y: -1.776e-015 Z: -1.395e+001Dist: 9.008e+001Rotation: X: 50.000 Y: 0.000 Z: 30.000Mag.: 1Ang.: 22.500Block Groupclaypile 4.3 单桩计算模型 FLAC3D 3.00Step 22973 Model Perspective11:26:04 Wed Feb 15 2012Center: X: -2.665e-015 Y: -1.776e-015 Z: -1.395e+001Dist: 6.905e+001Rotation: X: 50.000 Y: 0.000 Z: 30.000Mag.: 1Ang.: 22.500Plane Normal: X: 0.000e+000 Y: 1.000e+000 Z: 0.000e+000

Plane Origin: X: 0.000e+000 Y: 0.000e+000 Z: 0.000e+000Block Group Plane: onclaypileSketch Magfac = 0.000e+000 LinestyleItasca Consulting Group, Inc.Minneapolis, MN USA 4.4单桩计算模型中心切片图 材料参数的选取

表4.1为FLAC3D软件在建模与求解过程中所需要的材料参数：

表4.1 桩与桩周土体材料参数

混凝土桩 桩周土体

干密度 2500kg/m

3

弹性模量 25GPa 100MPa

泊松比 0.20 0.30

粘聚力 30kPa

内摩擦角 25°

1530kg/m3

数值模拟内容

本文运用FLAC3D软件对单桩在不同工况条件下进行数值模拟，主要对不同工况条件下单桩的水平位移、桩身弯矩等情况进行分析。

3 数值模拟结果及分析

3.1 水平力—位移梯度(H??Y0?H)关系曲线分析

单桩水平临界荷载和单桩水平极限承载力都能根据水平力—位移梯度关系曲线确定[8]。根据曲线的特征，取曲线的第一拐点所对应的荷载为水平临界荷载。取曲线的第二拐点所对应的荷载为水平极限荷载。

图3为三种不同工况的水平—位移梯度关系曲线。由图可看出带帽桩曲线的第一拐点和第二拐点都较普通桩曲线对应的荷载要大，说明带帽桩的水平极限承载能力高于普通桩，且达到水平临界荷载时的水平力要比普通桩大。带帽桩曲线较普通桩曲线平缓，即带帽桩对水平荷载增加导致位移增大的敏感性比普通桩弱，因此在受到水平荷载作用时也更稳定。

FLAC3D 3.00Step 22973 Model Perspective09:37:08 Wed Feb 15 2012Center: X: 1.356e+000 Y: 7.057e-001 Z: -1.255e+001Dist: 8.204e+001Rotation: X: 47.500 Y: 0.000 Z: 62.500Mag.: 1Ang.: 22.500ZYZFLAC3D 3.00Step 22973 Model Perspective09:30:20 Wed Feb 15 2012YXCenter: X: -2.665e-015 Y: -1.776e-015 Z: -1.395e+001Dist: 8.570e+001Plane Origin: X: 0.000e+000 Y: 0.000e+000 Z: 0.000e+000Rotation: X: 0.000 Y: 0.000 Z: 0.000Mag.: 1Ang.: 22.500Plane Normal: X: 0.000e+000 Y: 1.000e+000 Z: 0.000e+000Block Group Plane: onclaypileAxes LinestyleItasca Consulting Group, Inc.Minneapolis, MN USAFLAC3D 3.00Step 22973 Model Perspective09:35:15 Wed Feb 15 2012Center: X: 8.746e-001 Y: -5.049e-001 Z: -1.013e+001Dist: 3.157e+001Rotation: X: 100.000 Y: 30.000 Z: 30.000Mag.: 1Ang.: 22.500Interface LocationsItasca Consulting Group, Inc.Minneapolis, MN USAXContour of Y-Displacement Magfac = 0.000e+000-2.6531e-003 to -2.5000e-003-2.5000e-003 to -2.2500e-003-2.2500e-003 to -2.0000e-003-2.0000e-003 to -1.7500e-003-1.7500e-003 to -1.5000e-003-1.5000e-003 to -1.2500e-003-1.2500e-003 to -1.0000e-003-1.0000e-003 to -7.5000e-004-7.5000e-004 to -5.0000e-004-5.0000e-004 to -2.5000e-004-2.5000e-004 to 0.0000e+000 0.0000e+000 to 1.5897e-004 Interval = 2.5e-004Axes

FLAC3D 3.00Step 22973 Model Perspective09:50:12 Wed Feb 15 2012Center: X: -2.665e-015 Y: -1.776e-015 Z: -9.935e-001Dist: 6.938e+000Rotation: X: 0.000 Y: 0.000 Z: 0.000Mag.: 1Ang.: 22.500Contour of SYY Magfac = 0.000e+000 Gradient Calculation-9.0779e+004 to -7.5000e+004-7.5000e+004 to -5.0000e+004-5.0000e+004 to -2.5000e+004-2.5000e+004 to 0.0000e+000 0.0000e+000 to 2.5000e+004 2.5000e+004 to 5.0000e+004 5.0000e+004 to 7.5000e+004 7.5000e+004 to 1.0000e+005 1.0000e+005 to 1.2500e+005 1.2500e+005 to 1.2975e+005 Interval = 2.5e+004Itasca Consulting Group, Inc.Minneapolis, MN USAFLAC3D 3.00Step 22973 Model Perspective09:42:11 Wed Feb 15 2012Center: X: 4.000e+000 Y: 4.000e+000 Z: -1.400e+001Dist: 8.142e+001Plane Origin: X: 0.000e+000 Y: 0.000e+000 Z: 0.000e+000Rotation: X: 0.000 Y: 0.000 Z: 90.000Mag.: 1Ang.: 22.500Plane Normal: X: 1.000e+000 Y: 0.000e+000 Z: 0.000e+000Contour of Y-Displacement Plane: on Magfac = 1.000e+002 Exaggerated Grid Distortion -2.6099e-003 to -2.5000e-003-2.5000e-003 to -2.2500e-003-2.2500e-003 to -2.0000e-003-2.0000e-003 to -1.7500e-003-1.7500e-003 to -1.5000e-003-1.5000e-003 to -1.2500e-003-1.2500e-003 to -1.0000e-003-1.0000e-003 to -7.5000e-004-7.5000e-004 to -5.0000e-004-5.0000e-004 to -2.5000e-004-2.5000e-004 to 0.0000e+000 0.0000e+000 to 1.1941e-004Itasca Consulting Group, Inc.Minneapolis, MN USA

本文来源：https://www.bwwdw.com/article/62l2.html