做细胞力学仿真这几年FEBio一直是我主力工具之一。模型几何能画出来材料参数也能从文献里扒到但真正让很多新手卡住的往往是边界条件和载荷设置这一步——模型建得挺漂亮一提交计算就报错或者算出来结果明显违背物理直觉。这篇是这个系列的第六篇专门把边界条件和载荷设置这件事讲透。无论你是想模拟细胞压缩、基底拉伸还是微管吸吮实验这篇内容都能直接用上。先说个基本认知在FEBio里边界条件和载荷是两个不同的模块对应feb文件里BoundaryConditions和Loads这两个标签块。边界条件约束的是自由度载荷施加的是力或位移。理解这两者的区别是后面所有操作的基础。1. 先搞懂FEBio里边界条件与载荷的底层逻辑1.1 自由度、约束与载荷的关系FEBio底层是有限元求解器每个节点有6个自由度三个平动ux、uy、uz或者用sx、sy、sz表示和三个转动rx、ry、rz。对实体单元来说转动自由度在大多数情况下不参与计算所以你会看到很多模型里只约束平动自由度就够了。边界条件的本质是“告诉求解器哪些自由度被钉死了”。比如细胞底面黏附在基底上那这个面上节点的sy方向位移就是0你需要用fix把它的sy钉住。加载荷的本质是“告诉求解器哪些自由度被施加了外力或强制位移”。比如原子力显微镜探针压细胞探针接触区域会受到一个向下的力这就是节点力载荷。这里有个容易混淆的点强制位移prescribed displacement在FEBio里是加载荷模块里的不是边界条件模块。边界条件只管约束自由度载荷管施加力或位移。很多初学者把强制位移写成边界条件算出来结果完全不对。1.2 为什么FEBio要区分边界条件与载荷从有限元求解的角度看约束方程是在组装全局刚度矩阵之后、求解之前引入的它直接修改刚度矩阵的秩而载荷是作为右端项施加的它影响的是方程组的右端向量。这两者在求解流程中的位置不同混在一起会导致求解器无法正确识别。FEBio之所以这么设计是为了灵活处理非线性问题。在大变形、接触等非线性分析中约束和载荷的施加方式会随迭代过程变化。比如接触边界条件是在每个增量步内动态更新的如果你把它当成固定约束接触区域的力传递就会出现问题。我自己的习惯是先问自己“这个边界条件会不会随时间或变形状态改变”如果会大概率不是简单的fix能解决的要考虑用刚性约束、接触或者强制位移来处理。2. 边界条件设置从基础约束到接触边界2.1 最基础的fix约束与节点集选取fix是FEBio里最常用的边界条件用于固定节点的某些自由度。在feb文件里它的写法长这样BoundaryConditions fix bcfix namefixed_bottom node_setbottom_nodes dofssx, sy, sz/dofs /fix /BoundaryConditions其中node_set可以是你在PreView里预先定义好的节点集也可以直接在feb文件里用node_set标签临时定义。dofs里填的是需要约束的自由度sx、sy、sz分别对应三个平动方向。实际操作中我强烈建议所有边界条件都建立在命名清晰的节点集上。比如固定底面就定义bottom_nodes固定侧面就定义side_nodes。有个同事曾经图省事直接在几何上选了一堆节点结果后面想改约束范围找节点找了一个下午。命名规范这东西前期花五分钟后期省五小时。如果约束节点太多而且你是在PreView里手动选的注意检查有没有漏选。可以用PreView的“select by box”配合“highlight”功能检查节点集的覆盖范围。另外对于大模型的底面约束我更喜欢用面集来定义约束而不是节点集这样后续网格重划分时不容易出错。2.2 刚性约束rigid处理细胞-基底相互作用的利器细胞力学仿真里经常需要模拟刚性基底、刚性压头或者细胞内部的刚性结构。FEBio的rigid边界条件可以把一组节点绑定到刚体上刚体的运动可以由指定的参考节点控制。BoundaryConditions rigid_constraint bcrigid namerigid_indenter node_setindenter_surf rb_node100/rb_node /rigid_constraint /BoundaryConditions这里的rb_node是刚体参考节点的编号。刚体的平移和旋转自由度集中在这个参考节点上你可以对这个节点施加力或位移来控制整个刚体的运动。这个功能在做细胞压缩模拟时非常有用——压头是刚体细胞是变形体压头的下压量直接通过参考节点的位移载荷控制。用rigid而不是简单fix的好处是刚体可以有六个自由度的运动并且可以施加力和力矩而fix只是把自由度钉死不能施加载荷。如果你的模型里有个托盘或者压板需要边移动边压细胞rigid是首选。2.3 接触边界细胞与压头之间怎么定义接触是细胞力学仿真里最常用的非线性边界条件。接触不是直接约束自由度而是动态地在两个表面之间建立约束关系防止相互穿透。在FEBio里接触通过Contact模块定义常用的有contact类型下的sliding滑移接触和tied绑定接触。一个典型的接触定义长这样Contact contact typesliding namecell_indenter surfaceindenter_surf/surface surfacecell_top_surf/surface penalty0.01/penalty penalty_comp1/penalty_comp two_pass0/two_pass laugon1/laugon /contact /Contact这里有两个关键参数penalty是罚函数刚度laugon是是否启用拉格朗日乘子法进行增广迭代。罚函数法相对简单但容易穿透尤其是刚度不够时拉格朗日乘子法更精确但会增加计算量。对细胞这种比较软的材料我一般先用较小的penalty值试算如果发现穿透再逐步增大而不是一上来就设很大——罚函数太大容易让方程病态反而更难收敛。接触表面需要是面集不是节点集。在PreView里创建面集时注意面的法线方向——接触计算中主面和从面的法线方向会影响接触检测方向。如果法线反了接触会失效或者完全穿透。2.4 周期性边界条件单胞模型的必修课很多细胞力学模型为了提高计算效率只建一个代表性单胞然后施加周期性边界条件。FEBio从2.0版本开始支持周期性边界条件通过periodic约束实现。BoundaryConditions periodic bcperiodic nameperiodic_x surfacex_left surfacex_right dofssx/dofs /periodic /BoundaryConditions这个功能在做细胞微结构力学分析时特别有用。比如研究细胞骨架的等效力学性能建一个包含少量微丝、微管的单胞加周期性边界条件模拟无限大材料的响应。不过说实话周期性边界条件的设置比较繁琐需要把相对面上的节点一一对应。如果网格不是周期性生成的这一步会非常痛苦。我的建议是在做周期性子结构分析时网格划分阶段就要考虑周期性的节点匹配用PreView自带的周期网格工具或者用Python脚本提前生成匹配的网格。3. 载荷设置从节点力到位移加载3.1 载荷类型总览与选择逻辑FEBio的载荷模块支持多种类型节点力载荷nodal_load、面压力载荷surface_load、体载荷body_load、强制位移prescribed displacement、刚性体载荷rigid_load等。每种载荷类型对应不同的物理场景选错的典型后果是模型得不到收敛结果或者结果与预期相差一个数量级。我见过不少新手把面压力载荷用成节点力载荷结果力的总量完全不对。比如一个真实场景是给细胞顶面施加0.1 kPa的压力正确的做法是在顶面上定义面载荷而不是在节点上给力。原因很简单面载荷的总力等于压强乘以面积而节点力载荷是你手动指定每个节点上的力值两者对网格密度的依赖完全不同。具体来说细胞力学仿真中的载荷分两大类力载荷和位移载荷。力载荷模拟的是实际外力如流体剪切力、探针压力位移载荷模拟的是被控制的形变如基底拉伸量、压头下压距离。选择哪个要看你的实验条件对应哪个控制模式。原子力显微镜压痕实验常采用位移控制这时用强制位移载荷流体剪切实验常采用力控制这时用面压力或节点力载荷。3.2 节点力载荷nodal_load与单位换算节点力载荷是最直接的载荷形式。写法如下Loads nodal_load lc1 bcnodal nameforce_center node_setcenter_node dofssy/dofs value-1e-8/value /nodal_load /Loads这里的value单位是FEBio计算时的单位系统中的力单位。这是一个巨大的坑——FEBio本身不强制单位系统但你必须在整个模型中保持单位一致。如果你用微米μm做长度单位、毫牛mN做力单位那么刚度单位就是mN/μm应力单位是mN/μm²换算下来压强单位是kPa要特别注意。给个参照一个典型的细胞压缩实验细胞直径约10 μm探针施加的力在1~10 nN量级。如果你建的是轴对称二维模型力的单位换成nN/μm可以简化计算但要格外注意维度问题。我踩过最惨的一次坑是单位搞混算出来的应力大了三个数量级排查了两天最后发现是力的单位不对。我的建议是在模型文件夹里放一个README标明单位系统每次打开模型前先看一眼。3.3 面压力载荷surface_load与方向控制面压力载荷模拟的是垂直于表面分布的压强比如细胞浸润式流体压力、压缩板对细胞的接触压力如果用接触就不需要这个等。Loads surface_load bcpressure nametop_pressure surfacecell_top pressure lc10.001/pressure linear0/linear /surface_load /Loads这里pressure的值是压强值正负号决定方向正值指向表面法线方向负值指向法线反方向。linear设为1时压力随面积线性变化模拟的是非均匀压力场。面压力载荷最需要注意的还是法线方向。在PreView里创建面集时默认法线方向可能不是你想要的。检查方法很简单在PreView里显示面的法线方向如果指反了用“flip surface”工具反转。如果法线方向错了载荷会从向内压变为向外拉结果可能直接发散。3.4 强制位移载荷prescribe的应用原子力显微镜压痕模拟强制位移载荷在细胞力学仿真中使用频率极高。以原子力显微镜AFM压痕实验为例探针压入细胞表面一定深度这个深度就是控制变量对应的力是输出量。Loads prescribe bcprescribe nameindenter_displacement node_setindenter_nodes dofssy/dofs value lc1-0.5/value /prescribe /Loads这里的value是位移值单位与你的几何模型一致。如果几何模型单位是μm压入0.5表示下压0.5 μm。注意强制位移载荷施加的对象是节点集而且这个节点集通常就是刚性压头表面如果你用了rigid那就是刚体参考节点上的位移载荷。在非线性分析中强制位移载荷往往比力载荷更容易收敛。原因在于位移边界条件直接约束了位移场只要增量步足够小求解器不会因为载荷突变而发散。这也是为什么我在调试复杂模型时会先用位移控制做一遍确认模型稳定后再换成力控制。3.5 载荷曲线loadcurve控制载荷随时间变化FEBio的载荷曲线模块是控制所有时变载荷的核心。不管是力载荷还是位移载荷它们的值都可以通过载荷曲线随时间变化。LoadCurves loadcurve id1 typesmooth point0, 0/point point1, 0.2/point point2, 1/point /loadcurve /LoadCurvesloadcurve的type有四种step阶跃、linear线性、smooth平滑和linear_staircase线性阶梯。做细胞力学仿真我强烈推荐smooth类型因为它给出了一个连续可导的加载过程避免了阶跃载荷造成的计算不收敛。载荷曲线的x坐标单位是时间或者你求解器里设置的步数单位y坐标是载荷的缩放因子。举个例子你在nodal_load里设置了value-1e-8同时指定lc1那么实际施加的力是-1e-8 * 曲线y值。曲线在t0时为0t2时为1对应载荷从0逐渐增加到完整值。有个细节容易被忽略如果value没有加lc属性那么这个值是恒定不变的如果加了lc则需要乘以曲线因子。有些复杂的模型里我会在同一根曲线上做多阶段加载——先加载、保载、卸载这样就能模拟完整的力学循环。4. 细胞力学典型场景实操细胞压缩模型完整配置4.1 场景定义与模型简化以最常见的细胞压缩实验为例细胞被放置在刚性平板上上方一个刚性压头向下压缩细胞。实验目标是获得力-位移曲线拟合细胞的弹性模量。模型可以做轴对称简化如果用二维轴对称单元或全三维模型。这里用三维模型演示单元类型用六面体或四面体都可以。几何尺寸按典型细胞~10 μm直径建模。4.2 边界条件与载荷的完整配置示例这个模型里边界条件有三处。底部细胞与刚性平板接触区域约束所有平动自由度模拟黏附在平板表面。顶部细胞与压头接触区域不直接约束自由度而是通过接触条件与刚性压头交互压头本身作为刚体在参考节点上施加位移载荷。对称面如果有约束法向位移。对应的feb文件片段BoundaryConditions fix bcfix namebottom_fixed node_setcell_bottom dofssx, sy, sz/dofs /fix rigid_constraint bcrigid nameindenter_rigid node_setindenter_surface rb_node2000/rb_node /rigid_constraint /BoundaryConditions Loads prescribe bcprescribe nameindenter_down node_setrb_node dofssy/dofs value lc1-1.0/value /prescribe /Loads Contact contact typesliding namecell_indenter_contact surfaceindenter_surface/surface surfacecell_top_surface/surface penalty0.01/penalty laugon1/laugon /contact /Contact LoadCurves loadcurve id1 typesmooth point0, 0/point point1, 1/point /loadcurve /LoadCurves这里有个关键点刚体参考节点2000同时出现在rigid_constraint的rb_node里和prescribe的node_set里。这表示压头是刚体参考节点在y方向下压1 μm。注意对刚体参考节点施加位移载荷时dofs里也可以包含旋转自由度激活后可以模拟压头倾斜下压。4.3 求解器步长与载荷曲线的配合边界条件设好了载荷曲线也设了还要在Control里把时间步长设好。一个经典的设置在模拟中是这样的Control time1/time dtmin1e-6/dtmin dtmax0.01/dtmax max_inc1000/max_inc solvermax_ugs/solver /Controlsolver有两种常用选择max_ugs无对称矩阵求解器和max_ur有理数迭代法。细胞力学模型通常尺寸不大几千到几十万元素max_ugs足够。如果模型太大可以考虑max_ur但收敛性控制更复杂。时间步长的配置逻辑time是总模拟时间dtmin和dtmax是自动时间步长的下限和上限。载荷曲线把1个时间单位内的载荷从0变到1那么dtmax0.01意味着最多100步加载完成。如果模型复杂接触收敛慢求解器会自动缩小时间步但下限不能低于dtmin否则报错。初次试算时我习惯把dtmax设得更小比如0.005让加载更平缓换稳定性。确认模型没问题后再放大时间步长提高效率。4.4 结果验证与合理性检查算完之后第一步不是看应力云图而是检查力-位移曲线是否合理。在FEBio的后处理PostView里可以输出接触面的反力画一条力-位移曲线。对于线弹性小变形这应该是一条近似直线斜率和弹性模量成正比。如果曲线形状怪异比如出现负刚度那多半是接触穿透或者约束设置有误。另外检查应变能密度分布是否连续。如果有孤立单元的应变能密度异常高说明网格质量在某个区域有问题或者边界条件导致应力集中。对于细胞这种软材料应力集中往往出现在压头接触边缘和底边固定区域这是物理合理的但应力值不应该高出周围几个数量级。5. 常见问题与排查技巧实录5.1 不收敛、负体积、穿透、错误结果排查用FEBio做细胞力学仿真绕不开下面这些常见问题。这里直接给排查顺序和解决思路。现象典型原因解决办法第一步就发散单位系统不一致、约束刚度冲突检查所有材料参数和载荷的单位检查是否同一自由度被重复约束接触穿透罚函数刚度太小、法线方向反了增大penalty检查并翻转接触面法线负体积tet单元变形过大、网格畸变减小加载增量步细化网格或改用六面体单元力-位移曲线震荡载荷曲线突变、时间步太大改用smooth类型载荷曲线减小dtmax刚体漂移刚体自由度没有完全约束检查刚体参考节点的约束和载荷设置这些问题的排查逻辑我总结成一句话先检查单位再检查网格然后检查约束最后检查载荷。单位问题会直接导致数值爆掉网格畸变会导致局部应变异常约束冲突会导致刚度矩阵奇异载荷设置错误往往是数值发散的直接诱因。5.2 接触设置里的经典坑接触是细胞力学仿真里最占调试时间的部分。这里分享几个我踩坑踩出来的经验。第一接触面的网格密度要匹配。主面通常选较硬的材料表面网格可以粗一点、从面较软的材料表面网格要细一点。如果两个面上网格尺寸差异太大接触搜索会出现漏点导致接触力分布不连续。我做细胞压缩时细胞顶面网格尺寸约为底面网格的一半以下这样接触压力分布才光滑。第二penalty不是越大越好。罚函数刚度取材料刚度的1~10倍通常就够。太大会增加刚度矩阵的条件数导致迭代收敛慢甚至发散。我的调试习惯是先用penalty0.001试算如果有穿透每次×10直到穿透消失。第三two_pass选项在需要更精确的接触压力分布时打开。这是指主面和从面互换后再执行一轮接触检测。代价是计算量增加但对接触压力分布精度有要求的场景非常值得。5.3 单位系统的通病前面提了多次单位问题这里再单独强调。FEBio不限制单位系统这既是优点也是坑。现在做个速查表长度力时间应力/模量典型场景μmnNskPa细胞力学推荐mmNsMPa组织力学mNsPa结构性分析用推荐的单位组合细胞模量1~100 kPa对应数值是1~100压力载荷0.1 kPa对应0.1都很好写。最忌讳的是长度用μm、力用N这种混搭数值会跨好几个数量级收敛性极难保证。我的做法是在PreView里新建模型时就设好单位命名时也在文件名里带上单位后缀比如cell_compression_um_nN.feb。这个方法听起来土但真的能避免很多后面排查时的痛苦。6. 系列后续与个人补充这里再补充一个后面会具体展开的细节。FEBio的边界条件与载荷设置其实和材料本构模型的选择有很强的耦合关系。比如细胞模型用超弹性本构如Ogden、Mooney-Rivlin那么加载方式对收敛性和结果的影响会比线弹性模型大得多。超弹性材料在大变形下体积不变特性要求载荷增量不要太大这也是为什么很多细胞模型里强制位移载荷比力载荷更好收敛。另外FEBio 3.0以上版本提供了更丰富的接触算法包括绑定接触、粘弹性接触等做细胞-基底黏附、细胞-细胞聚集这类场景时可以关注一下。新版本的contact参数里还多了fric_coeff摩擦系数设置模拟细胞与基底间的滑动摩擦时非常实用。如果你在把边界条件和载荷设置应用到自己的细胞模型时遇到了难题比如压头形状变了不知道怎么约束、或者想加一个剪切流场却不知道怎么施加速度边界建议先回到这篇文章把每个边界条件的物理含义捋一遍再对照FEBio官方文档里的输入格式逐项核对。边界条件和载荷设置这件事说到底是“给模型一个合理的物理环境”的技术活。把自由度、约束、载荷、接触这几件事想清楚你的FEBio模型就已经成功了一半。下一步就可以放心地去探索更复杂的多物理场耦合问题了。