简介DCB.rar是一份面向地质力学、土木工程及结构分析研究者的DDA非连续变形分析源程序包适合具备一定编程基础并希望深入理解块体系统静动力响应机理的用户。包内程序基于离散介质建模思想通过接触模型处理块体破裂、滑移与碰撞可用于地震、爆炸等动载下的变形预测及地质灾害评估。压缩包共32个文件以C语言源文件、Visual Studio工程文件vcproj/vcxproj/sln/filters及调试辅助文件sdf、ncb等为主整体约7.78MB结构紧凑便于直接打开工程并对照源码学习。目前已有171人学习适合用于二次开发、算法复现或课程实验。通过阅读和运行这套DDA程序使用者可以掌握非连续介质建模的核心步骤了解接触判断、时间积分与动力响应的实现细节为实际岩土工程中的稳定性分析提供可扩展的代码基础。1. 从“DCB.rar”到非连续变形这份 DDA 源程序在解什么问题手上捏着一份名为 DCB.rar 的源程序包压缩包名里同时写着 DDA 源程序、deformation、动力分析。这通常不是一份普通有限元代码而是用非连续变形分析Discontinuous Deformation Analysis模拟双悬臂梁试件从预制切口起裂、界面张开到整体变形全过程的工程程序。DCB 试件是岩石和混凝土 I 型断裂测试的标准几何两臂端部被拉开裂纹沿界面扩展在连续介质框架里裂纹一旦起裂就需要网格重构或扩展有限元这类特殊处理DDA 走的是另一条路把试件切成离散块体让裂纹只能沿块体界面出现用接触弹簧的开闭迭代判断界面何时退出工作。动力分析在这里并不是加载地震波而是指块体在惯性、阻尼和接触力共同作用下运动DDA 用隐式时步求解每个时步落在一个平衡态上宏观上能看到块体脱离、转动和重新接触。这篇文章写给手上有 DDA 源程序却不懂怎么配置输入文件、或者参数改一次就发散需要把 DCB 计算结果和实验曲线对上的人。2. DDA 为什么适合做 DCB 动力分析块体离散与开闭迭代的本质2.1 非连续变形分析与有限元、离散元的边界在哪里DDA 和有限元最大的区别在于放弃了对连续位移场的强制逼近。计算区域被离散成凸多边形块体块体内部允许弹性变形位移场用完整多项式逼近常见是一阶或二阶块体之间没有共享节点接触面上用弹簧传递法向和切向力。全局方程仍然写成[K]{D} {F}的形式但刚度矩阵[K]不是单纯由材料积分得到而是由块体弹性刚度、接触弹簧、惯性力、加载约束等多个子矩阵叠加装配而成。和离散元相比DDA 不是显式中心差分的刚体动力学块体自身的变形参与刚度矩阵求解方式是每个时步解一个隐式方程组因此时步可以取得相对大低速断裂问题不需要像 DEM 那样跑到微秒量级。这一点决定了 DCB 模拟的工具选型。双悬臂梁的两个臂在开裂过程中不是刚体运动而是伴随明显的弯曲变形和裂纹尖端应力集中块体内部应变场必须参与计算。离散元也能做断裂但通常把块体视为刚体弯曲效应只能靠大量小尺寸块体拼出来计算开销和接触数量都会失控。DDA 保留了块体变形自由度同时把非连续变形界面张开、滑移、脱离作为一等公民这正是它在岩石裂纹扩展问题里比连续类方法更自然的原因。2.2 最小势能原理和时步推进DDA 的“运动”是怎么走出来的DDA 的求解核心是每个时步构造一个最小势能泛函并对其求极值。泛函中包含块体弹性应变能、惯性力所做的功、外荷载功、接触弹簧势能以及位移约束的罚函数。对未知的位移增量求偏导并令其为零就得到线性方程组。一个时步结束更新块体顶点坐标和接触状态再进入下一时步。所以 DDA 的动力分析本质上是在时间域上一步一步静力推进每一步终点满足平衡宏观运动通过坐标更新累积出来速度和加速度用相邻时步的位移差分反算。因此 DDA 对“时步”的语义和显式程序完全不同。显式程序的时步由稳定性条件决定DDA 的时步要同时兼顾接触探测分辨率和动力响应精度后文会给出工程区间。源程序里每步解方程之前都要做一次“开闭迭代”上一时步的接触状态作为初值装配接触弹簧后求解再检查所有接触是否处于一致的状态。闭着的接触如果法向拉力超过抗拉强度就在本时步内把接触置为开开着的接触如果出现重叠就重新置为闭然后重新组装、重新求解如此循环直到没有接触状态翻转为止。DCB 受拉时界面接触从前端逐个打开每打开一个荷载就转移给下一个接触宏观力-位移曲线上的渐进降载或突然跌落都来自这个开闭切换过程。2.3 DCB 试件建模的第一个决定裂纹路径怎么预埋DDA 里的裂纹只能沿块体界面扩展这是所有用 DDA 做断裂模拟的人都必须接受的第一性约束。因此 DCB 模型的第一步不是找网格生成器而是先画“潜在裂纹路径”。常见做法是把两臂之间的对称面从预制切口端开始切出一排细长块体或一条薄层界面切口尖端附近块体加密远离尖端逐渐稀疏。块体尺寸的过渡不能太陡因为接触弹簧刚度按接触边长折算相邻块体尺寸差一个量级就会让刚度矩阵出现悬殊项解出的位移会在界面处振荡。块体尺寸的起步值一般取 DCB 臂厚 h 的 1/10 到 1/20。块体太粗裂纹尖端的应力梯度捕获不足峰值荷载偏高块体太细时步必须跟着缩小计算量上升接触数量多到难以定位错误。真实 DCB 实验里切口是直的DDA 模型可以沿切口延长线预埋一条平直界面如果关心斜向断裂路径就需要把整个可能破裂区域都切成随机或规则块体让 DDA 的开闭迭代自己选择最薄弱路径但这会显著增加接触搜索成本。这个权衡在建模阶段就要定下来不要指望算法自动穿块破裂。3. DCB 模型的 DDA 源程序组织输入文件、核心子程序和数据流3.1 拿到 DCB.rar 后先做的三件事确认编译链、读目录、跑自带算例解压之后不要急着改参数。DDA 源程序的发行版本极多最常见的是一套 Fortran 77 老代码不同渠道流传的包在文件命名上差异很大但数据流基本一致一个主程序控制时步循环若干子程序分别负责输入读取、刚度装配、接触判断、位移更新和结果输出。常见结构是main.for或dda.for负责总控readin.for读取ddainput.datcontact.for做开闭判断output.for写结果文件还有一批编译脚本和示例输入文件。拿到包后先做三件事。第一确认编译器能把这套老代码编译通过用 gfortran 编译时建议打开越界检查选项老代码里的数组越界很常见开-fcheckall能提前暴露问题。第二找出输入文件的准确命名和数据格式不同版本对DDAINPUT.DAT的字段顺序有出入必须以源程序里的 READ 语句为准不能照抄网上的行号说明。第三完整跑一遍自带算例确认输出文件里的STEP分隔行和块坐标格式符合后处理脚本预期。这三步做完再替换成自己的 DCB 模型否则改了参数发散分不清是模型问题还是程序问题。3.2 ddainput.dat 关键参数检查表材料、时步、接触刚度与加载DCB 模型涉及的输入参数不算多但互相耦合。以下是一份经验性的参数检查表以岩石 DCB 工况为例实际取值必须以你手上源程序的量纲为准参数作用常规量级改错的典型现象块体数、节理数决定模型规模和接触探测范围与网格划分一致读取报错或算到一半崩溃弹性模量 E、泊松比 ν块内变形和接触刚度基底E 10~60 GPaν 0.2~0.3整条荷载曲线等比例偏移密度 ρ惯性力项和弹性波速2400~2700 kg/m³应力波到达时间错乱抗拉强度 σt界面开裂阈值按劈裂试验取值开裂提前或严重滞后粘聚力 c、摩擦角 φ压剪接触的破坏判据按直剪试验取值裂纹路径偏移时间步 Δt每步推进的时间增量10⁻⁴~10⁻³ s大穿透或计算量暴涨接触弹簧刚度 K界面罚刚度E 的 10~50 倍量级块体穿透或解振荡重力加速度 g体积力开关0 或 9.8整体位移漂移这张表刻意不写行号。DDA 老代码的输入文件格式在不同版本里差得很远有的版本第一行是控制开关有的版本第一行是块体数目照抄任何一套行号都可能在下一个版本失效。正确做法是打开源程序里读文件的子程序对照 READ 语句里的格式描述符逐个确认字段顺序然后用正文的参数表核对功能。K 值尤其要单独验证因为它不参与破坏判据只负责阻止穿透改错不会被程序报错只会在结果里留下莫名其妙的变形。3.3 用 Python 后处理从 OUT 文件提取 DCB 张开位移DDA 的输出文件通常按时步分段每个 STEP 段后跟着所有块体的顶点坐标。老版本输出格式相似都能用按行扫描的方式解析。用一段简单的 Python 提取每个时步的块体顶点坐标import re def parse_dda_out(path): steps [] cur None with open(path, r, errorsreplace) as fp: for line in fp: if line.startswith(STEP): cur {blocks: []} steps.append(cur) continue if cur is None: continue if line.strip().startswith(BLOCK): cur[blocks].append([]) continue nums re.findall(r-?\d\.\d(?:E[-]?\d)?, line) if len(nums) 2 and len(nums) 6 and cur[blocks]: cur[blocks][-1].append([float(nums[0]), float(nums[1])]) return steps逻辑说明程序按行扫描输出文件遇到STEP开头的新块就创建一条记录遇到BLOCK行说明接下来的数据属于一个新的块体常规坐标行的浮点数量在 2 到 6 个之间用正则抽取前两位作为 x、y 坐标。解析完成后指定左右臂端点的块体编号就能提取出每个时步的张开位移取左臂端点块和右臂端点块在切口位置的x坐标值两者差再减去初始间隙就是当前时步的张开距离。加载点位移同理用加载块体的顶点位移代替。这个后处理步骤要先跑通之后调参才有客观依据。4. 动力分析的参数校准与常见陷阱时步、刚度、阻尼的相互牵制4.1 时步上界弹性波穿越最小块体的时间动力分析时步的第一个约束来自波动传播。如果时间步大于弹性波穿过最小块体所需的时间力的传递在空间上就会延迟接触探测会漏掉本应发生的碰撞。时步上限按Δt_max L_min / Cp估算L_min是最小块体边长Cp是纵波波速。对岩石材料E 取 30 GPa、密度取 2600 kg/m³ 时Cp 约 3400 m/s如果最小块体边长是 2 mm时步上限约 0.0006 s。DDA 虽然是隐式求解接触探测仍然依赖几何位置更新时步超过这个上限会导致块体在一个步内越过彼此接触搜索直接失效。实际调试时建议从上限的一半开始即 0.0002 到 0.0003 s再逐步放大观察曲线变化。时步缩小的代价在 DCB 模型里非常明显因为裂纹尖端附近预埋了细长块体最小边长往往只有 1 到 2 mm总时步数会达到几万步。这时要把输出间隔调大只记录每 50 或 100 步的坐标否则后处理文件会膨胀到难以打开的规模。时步修改后还要同步检查动力响应的频率成分如果加载速率较快结构性振动周期远大于时步曲线不会出现异常如果加载速率极慢总时步数乘以时步可能超出实际物理时间太多模型还没开裂就已经跑了几十万步。4.2 接触弹簧刚度的收敛区间穿透与振荡之间怎么找平衡接触弹簧刚度是 DDA 里最需要校准的参数。它不决定界面何时破坏只决定破坏前界面能承受多大相对位移所以本质是一个罚参数。刚度太小块体之间出现明显穿透接触力滞后裂纹路径会沿着穿透区串出去刚度太大罚刚度与块体弹性刚度量级悬殊矩阵条件数恶化开闭迭代在一个时步内反复振荡能量曲线出现尖峰。经验做法是把 K 取在 E 的 10 到 50 倍量级具体数值要根据网格尺寸折算成单位面积刚度再通过一次静力试算验证。验证方法很直接单独跑一个不含破坏的 DCB 模型施加一个固定荷载然后检查接触处的穿透量。穿透量应小于最小块体边长的 5%如果超过就加倍 K 再算直到穿透量进入容差范围如果穿透量已经很小但求解开始振荡就把 K 减半再观察振荡是否消失。要注意这个标定必须在和正式算例相同的网格密度下进行因为穿透量与接触长度和块体尺寸直接相关换网格就要重新标。标定完成后 K 值一般不再改动后续调开裂行为应该只动抗拉强度和粘聚力。4.3 阻尼和加载速率在动力分析里做出准静态的 DCB 试验DCB 室内试验是准静态加载而 DDA 天然是动力学演化这两者之间的落差靠阻尼和加载速率来弥合。常见做法有两种一是用位移控制加载把加载速率放慢使整个过程中动能峰值占系统总能量的比例控制在百分之几以内这样惯性效应可以忽略算出来的就是准静态响应二是在程序中加入全局粘性阻尼每时步更新速度时乘以一个接近 1 的系数工程上常用 0.9 到 0.99作用相当于给整个系统一个背景阻尼。两者也可以结合先用小阻尼配合慢加载跑出曲线再逐步去掉阻尼检查峰值荷载是否改变。阻尼过小的表现是力-位移曲线伴随明显的高频震荡两臂像弹簧一样反复弹跳阻尼过大的表现是峰值荷载被压低破坏过程变成缓慢滑移而不是突发开裂。识别标准是看荷载峰值阻尼只应该消耗动能和震荡不应该改变准静态破坏的荷载阈值。如果你发现加大阻尼后峰值荷载明显下降说明加载速率仍然太快应该继续放慢加载而不是靠阻尼硬压。真实动力标定场景则要反向操作阻尼必须对应材料实测的阻尼比不能用准静态那套取值。4.4 三种失败图像和对应排查飞块、锯齿、能量不守恒调试 DDA 算例时最常遇到三种典型失败现象和原因基本固定。第一种是块体飞出去表现为某个块体的坐标在相邻两个时步间突然跳到远离主体的位置通常由时步太大和接触刚度不足共同引起块体在一次穿透后没有被接触弹簧弹回直接脱离了接触搜索范围。处理办法是把时步减半、接触刚度加倍重新运行。第二种是力-位移曲线在后半段持续锯齿状震荡且不收敛开闭迭代在同一批接触上反复翻转常见原因是加载速率过快或抗拉强度与接触刚度搭配不当先放慢加载再看是否需要把 K 下调半个量级。第三种是系统总能量随步数持续上升没有收敛趋势这通常是罚函数接触的伪能量在累积检查输出间隔和阻尼系数确认能量统计口径里包含了所有接触做功项。5. 从开裂荷载反推 I 型断裂韧度从 OUT 的 P-δ 曲线到 K_ICDCB 模型跑通后最终要输出的是一条荷载-张开位移曲线曲线峰值对应开裂临界状态峰值荷载结合试件几何就能换算 I 型断裂韧度。DDA 输出文件里通常只保存块体坐标和系统能量不直接保存加载点反力这时可以用外力功对位移做数值微分来还原荷载。外力功在每个时步有累计值张开位移已经从块体端点坐标提取出来两者对时间同步后P dW/dδ即可得到加载点荷载import numpy as np def load_from_work(steps): d np.array([s[opening] for s in steps]) w np.array([s[external_work] for s in steps]) p np.diff(w) / np.diff(d) d_center 0.5 * (d[:-1] d[1:]) return d_center, p逻辑说明np.diff计算相邻时步的外力功增量和张开位移增量两者相除得到该区间的平均荷载d_center取相邻位移的中点保证荷载与位移在同一位置对应。开裂前张开位移几乎为零外力功也没有增量数值微分结果是无意义的噪声开裂后才有连续曲线分析时从峰值所在的区段截取即可。因为数值微分放大了输出文件的离散误差通常还要对p做滑动平均处理窗口大小取总步数的 1% 到 2%足以压掉锯齿又不至于削平峰值。拿到峰值荷载后代入 DCB 的断裂韧度公式计算 K_IC。实际报告中更要紧的是验证这条公式用的峰值荷载是否可靠。两个检查必做一是收敛性检查把裂纹路径附近的块体尺寸从 h/10 加密到 h/20 再到 h/40峰值荷载的变化应逐渐收窄两次加密之间变化小于 5% 才认为网格无关二是能量守恒检查在无阻尼条件下跑同一个算例末时步的应变能加动能加摩擦耗能应等于外力功总和偏差超过 5% 就要回头查时步和接触刚度。DDA 给出的峰值荷载通常比实验值偏高原因是离散块体的裂纹尖端应力集中无法模拟真实材料的微损伤区域如果偏高太多优先检查预埋界面是否偏离了实际裂纹路径确认切口尖端附近网格密度足够。把不同网格尺寸对应的峰值荷载画成双对数曲线斜率趋平的位置对应的 K_IC 才是可以写进报告的值。本文还有配套的精品资源点击获取