前处理一切顺利，在执行到phase 4时程序自动退出，计算结果也没有保存。第一反应是程序出现了问题，检查一遍未发现异常。然后尽可能关闭后台进程、及服务，还是无法通过计算。

最后，关闭了NOD3.2 防火墙后，Plaxis终于执行完毕，计算结果如预期一样。这确实一个很奇怪的现象，NOD32在无警告的情况下强行关闭了fns.exe计算进程，不知这是不是Plaxis的bug、还是防火墙的bug。

第四课整个计算过程大约10分钟，笔记本配置为：酷睿双核2.0，内存2.0G。

国际上早20世纪60年代初就投入大量的人力和物力开发具有强大功能的有限元分析程序。其中最为著名的是由美国国家宇航局（NASA）在1965年委托美国计算科学公司和贝尔航空系统公司开发的NASTRAN有限元分析系统。目前的通用有限元程序主要有德国的ASKA、英国的PAFEC、法国的SYSTUS、美国的ABQUS、ADINA、ANSYS、BERSAFE、BOSOR、COSMOS、ELAS、MARC和STARDYNE等公司的产品。
以下对一些常用的软件进行一些比较分析：

1. LSTC公司的LS-DYNA系列软件

ANSYS/LS-DYNA将显式有限元程序LS-DYNA和ANSYS程序强大的前后处理结合起来。用LS-DYNA的显式算法能快速求解瞬时大变形动力学、大变形和多重非线性准静态问题以及复杂的接触碰撞问题。程序可以用ANSYS建立模型，用LS-DYNA做显式求解，然后用标准的 ANSYS后处理来观看结果。也可以在ANSYS和ANSYS-LS-DYNA之间传递几何信息和结果信息以执行连续的隐式－显式/显式－隐式分析，如坠落实验、回弹、及其它需要此类分析的应用。该数值模拟软件实际上仍在诸多不足，特别是在爆炸冲击方面，功能相对较弱，其欧拉混合单元中目前最多只能容许三种物质，边界处理很粗糙，在拉格朗日——欧拉结合方面不如DYTRAN灵活。虽然提供了十余种岩土介质模型，但每种模型都有不足，缺少基本材料数据和依据，让用户难于选择和使用。

2. MSC.software公司的DYTRAN软件

当前另一个可以计算侵彻与爆炸的商业通用软件是MSC.Software Corporation ( MSC公司) 的MSC.DYTRAN程序。该程序在是在LS-DYNA3D的框架下，在程序中增加荷兰PISCES INTERNATIONAL公司开发的PICSES的高级流体动力学和流体——结构相互作用功能，还在PISCES的欧拉模式算法基础上，开发了物质流动算法和流固耦合算法。MSC.DYTRAN 是具有快速模拟瞬态高度非性结构、气体/液体流动、流体-结构相互作用，采用显式积分技术，支持材料模型和高度组合非线性分析及流体-结构的全耦合。特别对高速碰撞、结构大变形和瞬时内发生的流体结构相互作用事件的仿真。
MSC.DYTRAN的算法基本上可以概况为：MSC.DYTRAN采用基于Lagrange格式的有限单元方法（FEM）模拟结构的变形和应力，用基于纯Euler格式的有限体积方法（FVM）描述材料（包括气体和液体）流动，对通过流体与固体界面传递相互作用的流体—结构耦合分析，采用基于混合的Lagrange格式和纯Euler格式的有限单元与有限体积技术，完成全耦合的流体-结构相互作用模拟。MSC.DYTRAN用有限体积法跟踪物质的流动的流体功能，有效解决了大变形和极度大变形问题，如：爆炸分析、高速侵彻。
但MSC.DYTRAN本身是一个混合物，在继承了LS-DYNA3D与PISCES的优点同时，也继承了其不足。首先，材料模型不丰富，对于岩土类处理尤其差，虽然提供了用户材料模型接口，但由于程序本身的缺陷，难于将反映材料特性的模型加上去；其次，没有二维计算功能，轴对称问题也只能按三维问题处理，使计算量大幅度增加；在处理冲击问题的接触算法上远不如当前版的LS-DYNA3D全面。

3. HKS公司的ABAQUS软件

ABAQUS是一套先进的通用有限元系统，也是功能最强的有限元软件之一，可以分析复杂的固体力学和结构力学系统。ABAQUS有两个主要分析模块：ABAQUS/Standard提供了通用的分析能力，如应力和变形、热交换、质量传递等；ABAQUS/Explicit应用对时间进行显示积分求解，为处理复杂接触问题提供了有力的工具，适合于分析短暂、瞬时的动态事件，但对爆炸与冲击过程的模拟相对不如DYTRAN和LS-DYNA3D 。

4 ADINA

ABAQUS 其解决问题的范围从相对简单的线性分析到许多复杂的非线性问题。 ABAQUS 包括一个可模拟任意几何形状的单元库。并拥有各种类型的材料模型库，可以模拟典型工程材料的性能，其中包括金属、橡胶、高分子材料、复合材料、钢筋混凝土、可压缩超弹性泡沫材料以及土壤和岩石等地质材料。作为通用的模拟工具， ABAQUS 除了能解决结构（应力/位移）问题，还可以模拟其他工程领域的许多问题，例如热传导、质量扩散、热电耦合分析、声学分析、岩土力学分析（流体渗透 / 应力耦合分析）及压电介质分析。
大量的复杂问题可以通过选项块的不同组合很容易的模拟出来。例如，对于复杂多构件问题的模拟是通过把定义每一构件的几何尺寸的选项块与相应的材料性质选项块结合起来。在模拟非线性问题，用户只需提供一些工程数据，像结构的几何形状、材料性质、边界条件及载荷工况，ABAQUS 能自动选择相应载荷增量和收敛限度。他不仅能够选择合适参数，而且能连续调节参数以保证在分析过程中有效地得到精确解。用户通过准确的定义参数就能很好的控制数值计算结果。

5 ANSYS和NASTRAN

因为和NASA的特殊关系，msc nastran在航空航天领域有着崇高的地位。ANSYS软件是融结构、流体、电场、磁场、声场分析于一体的大型通用有限元分析软件。由世界上最大的有限元分析软件公司之一的美国ANSYS开发，在铁道，建筑和压力容器方面应用较多。它能与多数CAD软件接口，实现数据的共享和交换，如Pro/Engineer, NASTRAN, Alogor, I－DEAS, AutoCAD等， 是现代产品设计中的高级CAD工具之一。软件主要包括三个部分：前处理模块，分析计算模块和后处理模块。分析计算模块包括结构分析（可进行线性分析、非线性分析和高度非线性分析）、流体动力学分析、电磁场分析、声场分析、压电分析以及多物理场的耦合分析，可模拟多种物理介质的相互作用，具有灵敏度分析及优化分析能力；

后处理模块可将计算结果以彩色等值线显示、梯度显示、矢量显示、粒子流迹显示、立体切片显示、透明及半透明显示（可看到结构内部）等图形方式显示出来，也可将计算结果以图表、曲线形式显示或输出。

软件提供了100种以上的单元类型，用来模拟工程中的各种结构和材料。该软件有多种不同版本，可以运行在从个人机到大型机的多种计算机设备上，如PC，SGI，HP，SUN，DEC，IBM，CRAY等。

国际上数值模拟软件发展趋势

1. 由二维扩展为三维
   现在，计算程序一般都由二维扩展到了三维，如LSDYNA2D和LSDYNA3D，AUTODYN2D和AUTODYNA3D，但也有完全在三维基础上开发的，如MSC.DYTRAN，就没有二维功能。
2. 从单纯的结构力学计算发展到求解许多物理场问题
   数值模拟分析方法最早是从结构化矩阵分析发展而来，逐步推广到板、壳和实体等连续体固体力学分析。近年来数值模拟方法已发展到流体力学、温度场、电传导、磁场、渗流和声场等问题的求解计算，最近又发展到求解几个交叉学科的问题。例如内爆炸时，空气冲击波使墙、板、柱产生变形，而墙、板、柱的变形又反过来影响到空气冲击波的传播……这就需要用固体力学和流体动力学的数值分析结果交叉迭代求解，即所谓“流—固耦合”的问题。
3. 从单一坐标体系发展多种坐标体系
   数值模拟软件在开始阶段一般采用单一坐标，或采用拉格朗日坐标或采用欧拉坐标，由于这两种坐标自身的缺陷，计算分析问题的范围都有很大的限制。为克服这种缺陷，采用了三种方法，一是两个程序简单组合，如CTH—EPIC，爆炸与侵彻由不同的程序分开计算；二是在同一程序中采用多种坐标体系，如DYNA3D中早期采用的是拉格朗日坐标，而LSDYNA3D的最新版除原有类型外，新加了欧拉方法以及拉格朗日与欧拉耦合方法，而最近几年才发展的DYTRAN则是拉格朗日型的LSDYNA3D(1988版)与欧拉型的PISCES的整合体；三是采用新的计算方法，如SPH等，SPH法不用网格，没有网格畸变问题，所以能在拉格朗日格式下处理大变形问题，同时，SPH法允许存在材料界面，可以简单而精确地实现复杂的本构行为，也适用于材料在高加载速率下的断裂等问题的研究。
4. 由求解线性工程问题进展到分析非线性问题
   诸如岩石、土壤、混凝土等，仅靠线性计算理论就不足以解决遇到的问题，只有采用非线性数值算法才能解决，它涉及到很多专门的数学问题和运算技巧，很难为一般工程技术人员所掌握。为此，国外一些公司开发了诸如LSDYNA3D、ABAQUS和AUTODYN等专长于求解非线性问题的有限元分析软件，并广泛应用于工程实践。这些软件的共同特点是具有高效的非线性求解器以及丰富和实用的非线性材料库。
5. 增强可视化的前置建模和后置数据处理功能
   目前几乎所有的商业化数值模拟程序系统都有功能很强的前置建模和后置数据处理模块。在强调“可视化”的今天，很多程序都建立了对用户非常友好的GUI（图形用户界面—Graphics User Interface），使用户能以可视图形方式直观快速地进行网格自动划分，生成有限元分析所需数据，并按要求将大量的计算结果整理成变形图、等值分布图，便于极值搜索和所需数据的列表输出。
6. 与CAD软件的无缝集成
   即在用CAD软件完成结构设计后，自动生成有限元网格并进行计算，如果分析的结果不符合设计要求则重新进行构造和计算，直到满意为止。今天，工程师可以在集成的CAD和数值模拟软件环境中快捷地解决一个在以前无法应付的复杂工程分析问题。所以当今所有的商业化有限元系统商都开发了和著名的CAD软件（例如AutoCAD、Pro/ENGINEER、Unigraphics、SolidEdge、SolidWorks、IDEAS等）的接口。
7. 工作平台多样化
   早期的数值分析软件基本上都是在大中型计算机上开发和运行的，后来又发展到以工程工作站（EWS，Engineering Work Station）上，它们的共同特点都是采用UNIX操作系统。Microsoft Windows操作系统和32位的Intel Pentium处理器的推出，为PC机用于有限元分析提供了必需的软件和硬件支撑平台。因此当前国际上著名的有限元程序研究和发展机构都纷纷将他们的软件移值到Windows平台上。
8. 软件开发强强联合
   随着数值模拟软件商业化的进展，各数值模拟软件公司为扩大市场，追求共同的利润，出现了强强联合的局面。典型的如ANSYS与LSDYNA3D联合，MSC.software软件公司对ABAQUS、LS DYNA3D及PISCES等的购买。

后记：其实目前仿真软件数量远远不止上面所列举内容，每个专业领域都有各自的仿真程序，如岩土工程领域的FLAC3D（有限差分），PLAXIS等等。无论那种手段都有其优缺点，熟练掌握好2个工具足矣。所谓的“通用”有限元软件，在解决实际问题中效果未必跟得上“专业软件”，但是作为工作平台仍然可以。另外，以上大部分内容来源于互联网，未能找到原文作者，所以在文章最后无法提供参考出处。

参考：[1] LS-DYNA使用指南.

1. 程序提供结构单元(element)、网格(zone)等数值手段构建数值单元集合体以表征复杂性状模拟对象。FLAC中内置了完善的线性/非线性应力/应变屈服及流动法则，完美解决了大变形这一复杂问题的数值实现，并且，拉尔朗日算法和混合离散技术的应用确保塑性屈服/流动过程算法的完备性和结果的精确性；
2. 快速拉格朗日法不需要组集模型的整体刚度矩阵，避免了大规模计算机内存消耗，这一特点大大提高了数值计算的效率；
3. 与FEM等隐式方法不同地，计算过程中引入了节点质量缩放及阻尼技术使得计算效率进一步提升，而这一技术的应用对材料破坏模式并无影响；
4. 尽管FLAC设计开发之初，其目标应用对象是岩土工程、采矿工程师，但FLAC的强大功能已被广泛应用于力学领域复杂问题的分析，缘于FLAC内置了完善、丰富的材料库，为用户模拟高度非线性、力学响应不可逆材料提供了强大工具；
5. FLAC内置了强大程序编译器FISH，用户可以使用FISH编译器定制自定义函数满足特定功能，甚至是自定义材料本构模型。高级用户通常期望与程序内核实现交互访问的意图，FISH编译器无疑为此提供了快捷途径。

FLAC3D - 基于FLAC连续介质理论的三维分析、试验、和设计软件，被岩土工程、土木工程和采矿工程师广泛用于模拟FEM理论所不能企及的复杂问题，如：复杂工况问题、大变形问题、复杂非线性行为和不稳定系统问题（甚至是大范围屈服/破坏，或坍塌）。

UDEC - (通用离散单元程序，Universal Distinct Element Code) 是为涉及节理岩体、散体材料准静态或动态问题分析量身定制的一款离散单元法程序。目前，UDEC主要用于采矿、核废料处理、能源回收、岩质基础坝体稳定、地震/微震和深埋地下结构等一系列超常规问题。

3DEC - (三维离散单元法程序，3-Dimensional Distinct Element Code) 是基于ITASCA-UDEC理论并扩展至三维情形的离散单元法程序，专为模拟节理岩体、散体材料准静态、动力模拟而开发。

PFC2D & PFC3D - (颗粒流分析程序，Partical Flow Code) 特别用于模拟任意性状、大小的二维圆盘或三维球体集合体的运行及其相互作用的强大颗粒分析程序。除了模拟大体积流动和混合材料，程序更适合于描述在固体材料中细观/宏观裂纹扩展、破坏累积并断裂、破坏冲击和微震响应等高水平课题的深化研究。

岩土工程专业软件FLAC、FLAC3D、UDEC、3DEC、PFC2D、和PFC3D，这些程序的共同特点如下：

1. 针对岩土体问题开发、但不限于岩土体问题
2. 可以解决大变形、甚至几何形态破坏问题
3. 可以追踪记录破坏过程
4. 多种岩土本构
5. 地质结构面模拟
6. 真时间历程动力模拟
7. 地下水模拟
8. 复杂施工过程模拟(分布开挖、回填、加固)
9. 多种加固措施
10. 内置外接程序语言满足用户特定要求

除了用于高级数值分析(建模、连续非线性分析、非连续非线性分析、非连续刚性块体分析)之外，广泛应用于解决采矿工程、爆破工程、水电边坡工程、工程可视化、检测设计。

参考：ITASCA.CN

1. 检查模型的运行时间 FLAC3D运行的求解时间与N ^4/3成比例，其中N是单元的数目。这个关系式仅适用于求解弹性问题达到平衡, 对于塑性问题，运行时间会有些变化，但变化不大。但如果发生持续的塑性流动则时间要比这多得多。.对某一特定问题，检查你的计算机的计算速度非常重要，运行手册中的测量计算速度的例子（主要使用了FISH中的CLOCK），得出你的计算机的运行速度，然后对单元数目进行插值，就可以估计出该模型的计算速度。

2.运行时间的影响 如果出现下面的情况，FLAC3D将需要较长的时间才能达到收敛： (a) 单元材料的刚度差异较大，或者单元、结构元和界面的刚度差异很大；或者 (b) 单元尺寸的差异较大。当这些刚度差异较大时，FLAC3D变的不是非常有效。应当在进行详细的分析之前调查刚度的对比影响。例如，非常刚性的载荷面可以用一系列的具有恒定速度的网格点来代替。(记住：FIX命令固定的是速度而不是位移) 包含地下水将增加力学体积模量。

3. 考虑单元密度 FLAC3D 使用常应变单元。如果应力/应变的梯度很高，你将需要更多的单元来表示这种变化的分布。可以对同一问题使用不同的单元密度来检查这种影响。FLAC3D之所采用常应变单元是因为模拟塑性流动时，用较多的低阶单元比用较少的高阶能够达到更好的精度。尽量保持单元尺寸均匀，特别是那些感兴趣的区域。避免纵横比大于5:1的长而窄的单元，避免单元尺寸的跳跃（也就是说，使用平滑过渡的网格），在GENERATE命令中用ratio关键字平滑地从细网格区域到粗网格区域过渡。

4. 平衡状态的自动监测 缺省情况下，使用SOLVE命令会自动监测力学的力的平衡。当模型中所有网格点的最大不平衡力值与模型中所有节点的平均施加力的比值低于1 ×10^-5时，就认为达到了平衡状态。注意网格点施加的力由内力(例如由于重力载荷)和外力(例如由于施加应力边界条件)共同作用。在大多数情况下，这种不平衡力与施加力之比的定义对静态平衡提供了精确的的界限，并且由于这个比值是无量纲的，该极限值适合于不同单位系统下建立的模型。 通过SET ratio命令可以用户定义不同的ratio值，如果默认的比例极限不能为静态平衡提供足够准确的界限，那么应当检查更改的比例极限。默认的比例极限也可用于监测热计算和流体流动计算的稳态解。对热计算来说，使用不平衡的热流量值和施加的热流量值来评价；对流体流动计算来说，使用的是不平衡的流体流动量和施加的流体流动量来评价。

5. 考虑选择阻尼 静态分析中默认的力学阻尼是局部阻尼(local damping)，当大多数节点的速度分量周期性地过穿零值时，用局部阻尼来移去动能是最有效的办法。这是因为质量调整过程取决于速度符号改变。为了不引入错误的阻尼力而达到静态平衡解，局部阻尼是一种非常有效的算法。

如果问题中有重要部分的网格在解答的最后状态具有非零的速度分量, 那么默认的damping不适合去达到平衡状态.一种不同的damping形式,即 combined damping, 在网格中有重要的刚体运动发生时, 能比local damping提供更好的收敛以达到稳定状态 . This may occur, 例如, 在蠕变模型中或在计算轴向受荷的桩的极限承载力时. 用命令 SET mechanical damp combined 去触发combined damping. Combined damping 在removing kinetic energy不那么有效, 因此要注意minimize dynamic excitation of the system (see Example 3.14). 用命令SET mechanical damp local可以转回默认的damping .

6. Check Model Response 分析前做简单的tests确认模型如所想建立的那样. For example, if a loading condition and geometry are symmetrical, make sure that the response is symmetrical. 在模型上做出变动后, 执行几步 (say, 5 or 10)去验证初始反应和位置是否正确. 分析应力或位移值的发展情况, 并与FLAC3D的输出结果比较. 如果对模型施加什么样的条件, 就会得到什么样的response. 例如, 如果模型各方面显示结果比较理想, 但某角的网格具有较大的 velocities, 一定要先找出原因再继续进行运算. 此种情况有可能是对某边界节点没有固定.

7. 变量初始化 在运算期间常需对节点位移初始化以便实现多次不同开挖阶段的模拟. 因为运算过程中不需要位移值因此可以这样做 .velocities的初始化则有点困难. 如果节点velocities固定为一恒定值, 那么将一直持续为该值除非给它设置其它值. 因此, 不要将网格的velocities初始化为零 — 这将会影响模拟结果. 但是, 有时设velocities为零很重要(for example, to remove all kinetic energy).

8. Minimizing Transient Effects on Static Analysis 对一系列的静态分析, 在不同阶段逐步得到解答很重要—i.e., make the solution more “static” by minimizing the effects of transient waves when problem conditions are changed suddenly. 有两种方法使FLAC3D 解答更静态化. 1. 当有突变产生时 (e.g., by nulling zones to simulate excavation), 将强度参数设置为较大值并step至平衡. 再将参数设为真实值,再次 step以确保不平衡力很小. 这样, 在过渡时不会激发破坏产生. 2. 当移除材料时, 用FISH 函数或 table history 逐步减少荷载 (e.g., see examples in Section 3.3.1.2).

9. Changing Material Models FLAC3D 对采用的模型材料的数量没有限定. 这种原则允许用户对对最大尺寸网格的每一网格采用不同的材料 (如果需要的话).

10. Running Problems with In-situ Field Stresses and Gravity 许多问题中都需要对模型施加原位应力和重力. 深矿井的开挖与回填就是其中一例, 其中岩体受原位应力场作用 (i.e., gravity stresses for the limited mesh size can be ignored), 但其中the emplaced backfill pillars会发展gravitational stresses而导致在荷载作用下发生破坏. 在这些模型(as in any simulation in which gravity is applied)中需要注意的重要一点是网格至少有三点在空间上必须固定—要不然, 整个网格会在重力作用下 translate. 如果你发现在重力加速度方向整个网格发生translating的话, 可能就是忘了在空间上 fix the grid(e.g.,see Example 3.16).

11. 确定破坏荷载 为确定破坏荷载, 最好采用“应变控制”边界条件而不是“应力控制” — i.e., apply a constant velocity and measure the reaction forces rather than apply forces and measure displacements. 当所施加荷载接近破坏荷载时, 系统很难控制, 不管对真实系统还是模型系统都是如此.

12. 安全系数的确定 FLAC3D中对任何参数都可以通过给定条件下的参数值与破坏时的参数值的比值来确定.例如, Fw =water level to cause collapse/actual water level FL =applied load to cause failure/design load Fφ =tan (actual friction angle)/tan (friction angle at failure) 注意往往是大值除以小值(假定系统在实际条件下不会破坏). 破坏的定义由用户建立. Dawson and Roth (1999) 及Dawson et al.(1999)对在强度折减基础上得到的安全系数的有限元解答作了比较.见section 1 in the Examples volume.

13. 使用体积模量和剪切模量 FLAC3D的弹性参数应该使用体积模量K和剪切模量G，这比使用杨氏模量E和泊松比v好。（K，G）适用于所有不违反热－动力学定律的弹性材料，而(E，v)不适用于某些材料。一个极端是体积变化但不发生剪切的材料；另一个极端是发生剪切但体积不变的材料。前者对应于限定的K和零G，后者对应于零K和限定的G。如果使用（E,v）就不能表示这两种材料。如果我们排除这两种特殊情形(也就是v=0.5和v=-1)，那么这两对数据集有如下的关系：

3K(1-2v)=E 2G(1+v)=E

只要不达到上述的两种极限情况，此关系式就有效。我们不必把他们与那些可行或不可行的物理试验相联起来，这些等式仅仅是两种定义比例系数的结果。假定我们有一种材料，变形的阻力逐渐减少，而体积变化的阻力保持恒定，在这种情况下v达到了0.5，还必须满足关系式3K(1-2v)=E。有两种可能性（争论在几何层面上而不是物理层面上）：或者E保持限定（非零），K趋于任意的大值，或者K保持限定，E趋于零。第一种可能性是因为对所有已知的材料有一个限制的压缩性，比如水2GPa，其泊松比是0.5，第二种可能性E发生彻底的变化，尽管我们假定材料弹性阻力的主要模式不改变。因此我们得出结论：用（E,v）表示材料的行为是不合适的。

原文作者：刘波博士