FLAC3D 实战:Python 记录峰值质点速度 PPV 与保存策略(同理可实现时程分析包络云图)

摘要： 本文详细讲解 FLAC3D 官方 PPV（峰值质点速度）算例，展示如何通过 Python 回调函数在动力计算中实时追踪每个网格点的最大速度矢量，并将 Python 自定义数据持久化保存到 FLAC3D 存档文件中。文末附改良版包络云图计算脚本，可直接用于时程分析，适合从事动力分析、爆破振动、地震响应的岩土工程师阅读。回调频率并非越高越好，过密会拖慢计算。（点赞 + 转发 + 关注本公众号，后台回复「PPV包络」即可获取完整改良版，直接用于时程分析包络云图）

字数： 约 3200 字 | 预计阅读时间： 8 分钟

🎯 问题背景

在岩土动力工程中，峰值质点速度（Peak Particle Velocity, PPV） 是评价振动效应最核心的指标之一。无论是爆破振动安全评估、地震响应分析，还是动力机器基础设计，PPV 都是判断结构是否安全的关键依据。

FLAC3D 本身并不直接输出 PPV 云图——动力分析完成后，我们通常只能看到某个时刻的速度场，而整个时程过程中的最大值需要额外记录。

这个算例的特别之处在于：

• ✅ Python 回调实时计算——每个力学步自动更新 PPV，无需后处理
• ✅ numpy 向量化运算——利用 gridpointarray 模块高效批量计算
• ✅ Python 数据持久化——将 Python 变量存入 FLAC3D save 文件，重启不丢失
• ✅ save / restore / new 全生命周期管理——保存、恢复、新建模型时数据自动同步

📐 问题描述

模型概况

算例建立了一个 100 m × 100 m × 100 m 的立方体弹性模型，中心区域挖去一个 10 m × 10 m × 20 m 的长方体孔洞，用来模拟地下开挖空间。在模型一角施加集中应力作为震源，应力波在模型内传播，记录所有网格点在整个时程中的速度峰值。

模型几何特征：

• 整体尺寸：100 m × 100 m × 100 m
• 网格划分：100 × 100 × 100（共 100 万个单元）
• 中心开挖：x ∈ [40,50]，y ∈ [40,50]，z ∈ [40,60]
• 震源位置：模型角部 (8081, 8081, 80~81) 单单元应力加载
• 计算时长：0.08 s

材料参数

边界条件与震源

• 安静边界（quiet boundary）：模型六个外表面均设置粘性安静边界，吸收反射波
• 震源加载：在角部一个单元内施加 10 MPa 等向初始应力作为扰动源
• 求解终止：总动力时间 0.08 s

🔧 FLAC3D 建模要点

1️⃣ 网格生成

zone create brick size 100 100 100
zone delete range position-x 40 50 position-y 40 50 position-z 40 60

• 立方体 100×100×100 均匀网格
• 删除中心长方体区域模拟开挖
• 单元尺寸 1 m，满足动力波长分析要求

2️⃣ 本构与动力配置

model config dynamic
zone cmodel assign elastic
zone property density 2950 young 12e9 poisson 0.25
zone dynamic damping rayleigh 0.05 150

• 线弹性本构，聚焦波动传播规律
• Rayleigh 阻尼 5%，中心频率 150 Hz

3️⃣ 安静边界设置

zone face apply quiet-normal range position-x 0 position-x 100 ...
    position-y 0 position-y 100 position-z 0 union
zone face apply quiet-strike range position-x 0 position-x 100 ...
    position-y 0 position-y 100 position-z 0 union
zone face apply quiet-dip range position-x 0 position-x 100 ...
    position-y 0 position-y 100 position-z 0 union

• 所有外表面施加法向、走向、倾向三向安静边界
• 使用 union 关键字合并 6 个面为一个 range

4️⃣ 震源加载

zone initialize stress xx 10e6 yy 10e6 zz 10e6 ...
    range position-x 80 81 position-y 80 81 position-z 80 81

• 单个单元内瞬时施加 10 MPa 等向压应力
• 应力骤变产生球面波向四周传播

📊 核心技术：Python 回调记录 PPV

原理

PPV 定义：每个网格点在整个动力时程中经历的速度矢量模的最大值。

PPV_i = max( |v_i(t)| )   对所有时间步 t

实现方式：通过 it.set_callback 注册一个 Python 函数，在每个力学计算步自动调用，用 np.maximum 逐元素取最大。

核心代码

import itasca as it
import numpy as np
from itasca import gridpointarray as gpa
ppv =Nonedefstore_ppv(*args):global ppv
if ppv isNone:        ppv = np.zeros(it.gridpoint.count())    ppv = np.maximum(np.linalg.norm(gpa.vel(), axis=1), ppv)it.set_callback("store_ppv",0.1)

关键解读：

• gpa.vel() 返回所有网格点速度的 (N, 3) numpy 数组
• np.linalg.norm(..., axis=1) 批量计算每个点的速度矢量模
• np.maximum 逐元素与历史最大值比较并更新
• 回调频率应兼顾精度与效率

🔄 进阶：Python 数据存入 FLAC3D Save 文件

痛点

默认情况下，Python 变量不会随 FLAC3D 保存文件一起存储。执行 model save 后再 model restore，Python 里的 ppv 数组就没了——这对于大型动力分析来说非常致命。

解决方案

利用 FLAC3D 的 gridpoint extra 变量作为桥梁，在保存时写入、恢复时读回。

defppv_save(*args):"""保存时触发：将 ppv 写入网格点 extra 变量"""global ppv
print("saving PPV")if ppv.any():        gpa.set_extra(1, ppv)defppv_restore(*args):"""恢复时触发：从 extra 变量读回 ppv"""global ppv
print("restoring PPV")try:        ppv = gpa.extra(1)print("max ppv during restore", ppv.max())except:passdefppv_new(*args):"""新建模型时触发：重置 ppv"""global ppv
print("resetting PPV")    ppv =None# 注册三个事件回调it.set_callback("ppv_save","save")it.set_callback("ppv_restore","restore")it.set_callback("ppv_new","new")

三个回调事件的时机：

💡 另一种方案：也可以用 it.fish 函数将 Python 值转为 FISH 变量保存，适合标量或简单数组。但对于大规模网格数据（如百万级 PPV 数组），用 extra 变量更直接高效。

📈 后处理分析：PPV 衰减规律

计算完成后，用 numpy + matplotlib 分析 PPV 随震源距离的衰减关系。

分类统计

将网格点分为三类：

load_point =(80.5,80.5,80.5)distance = np.linalg.norm(gpa.pos()- load_point, axis=1)gpos = gpa.pos()gx, gy, gz = gpos.T
# 地表点（顶面）top_mask = gz ==100.0# 外表面点exterior_mask =reduce(np.logical_or,(gx==0, gx==100, gy==0, gy==100, gz==0, gz==100))# 边界网格点（关联单元数 ≠ 8）boundary_gridpoint_mask = np.array([len(v)!=8for v in gpa.zones()])# 开挖面 = 非外表面 + 边界点excavation_surface_mask = np.logical_and(    np.logical_not(exterior_mask),    boundary_gridpoint_mask)# 内部点interior_mask = np.logical_not(boundary_gridpoint_mask)

衰减曲线