ABAQUS模型：基于CEL算法的砂土静力触探贯入分析。 使用ABAQUS有限元软件，利用CEL流固耦合算法计算了砂土的探针贯入模型，还原了清华大学的《静力触探贯入机理的数值模拟研究》硕士论文第3张模型，模型分为两个分析步，首先地应力平衡，然后是探针贯入模拟，分析了砂土的变形受力情况。 包括模型教程和文件。

这个模型核心在于控制欧拉体与拉格朗日体的交互。砂土层用欧拉网格避免畸变，探针作为刚体采用传统拉格朗日网格。地应力平衡阶段有个小技巧：先给砂土施加重力场，通过*Geostatic分析步让土体自然沉降。这里容易翻车的是材料参数设置，特别是砂土的Drucker-Prager模型参数：

mdb.models['Model-1'].materials['Sand'].DruckerPrager(
    table=((30.0, 0.7, 25.0), ), 
    deviatoricEccentricity=0.1,
    meridionalEccentricity=0.5
)

参数表里的0.7是材料硬化参数，实测发现这个值调低会显著影响贯入阻力曲线形态。建议先做参数敏感性分析，别像我第一次跑出个探针像捅豆腐的鬼畜结果。

贯入阶段的关键在于接触定义。CEL特有的*Coupling interaction里有个隐藏关卡——欧拉体边界处理。代码里这个设置直接影响砂土是否从模型侧面漏出：

mdb.models['Model-1'].EulerianSection(
    name='Section-sand',
    controls=EulerianBC(
        outflowFaces=(LEFT, RIGHT, FRONT, BACK), 
        inflowFaces=()
    )
)

记得把上表面设为OPEN允许土体隆起，下表面设为FIXED模拟硬底层。有个骚操作是在初始步设置0.99体积分数的砂土层，顶部留1%空隙给材料流动腾空间，实测能有效避免计算初期发散。

结果后处理时别光盯着位移云图。提取锥尖阻力qc有个隐藏路径：在历史输出里勾选探针反力RF，用Python脚本处理时域数据：

from odbAccess import openOdb
odb = openOdb('penetration.odb')
rf = odb.steps['Step-2'].historyRegions['Node PART-1-1.1'].historyOutputs['RF3'].data
qc = [abs(v[1])/1000 for v in rf]  # 单位换算kN→MN

注意探针锥角处的应力集中现象。对比清华论文的数据时发现，当网格尺寸缩小到0.2D（D为探针直径）后，qc值的误差能控制在15%以内。不过计算时间也从2小时暴涨到8小时，建议用多核并行计算。

模型文件已打包，包含完整的材料参数和边界条件设置。需要的小伙伴可以评论区留言，记得在Job提交前检查欧拉域尺寸是否足够——别问我怎么知道探针贯入时砂土从屏幕边缘喷出来的酸爽。