基于三场混合有限元的近不可压缩超弹性材料拓扑优化代码_非线性_拓扑优化_MATLAB_UG_通信

基于三场混合有限元的近不可压缩超弹性材料拓扑优化代码

一、摘要

本文提出了一种基于三场混合有限元的近不可压缩超弹性材料拓扑优化方法。采用Mooney-Rivlin模型来描述超弹性材料的本构关系。采用三场混合有限元方法解决了近不可压缩材料导致的体积闭锁问题。在有限变形假设下建立了以柔度最小为目标函数，体积分数为约束的拓扑优化模型。采用能量插值方案来规避低密度单元引起的数值收敛问题。利用伴随方法获取灵敏度信息。数值算例研究了结构的载荷大小对拓扑优化结构的影响，通过与基于位移场有限元的拓扑优化结果进行对比验证了所提方法的有效性。最后，提供了完整的 299 行 MATLAB 代码及详细说明。

吉林大学左文杰教授团队撰写了“299-line topology optimization code of nearly incompressible hyperelastic materials using three-field mixed finite element”论文，并发表于Computational Mechanics期刊，其中通信作者为吉林大学左文杰教授和白建涛副教授，第一作者为吉林大学博士生郭会强。

二、研究背景

不可压缩材料在整个变形过程中保持恒定的体积，这是连续介质力学和计算力学中常见的理想假设。常见的工程结构如垫圈、衬套、非充气轮胎、人造血管等，往往表现出近似不可压缩或不可压缩的特性，如图1所示。因此，在设计此类工程结构时，考虑其不可压缩特性是非常重要的。

(a)垫圈和衬套

(b)免充气轮胎

(c)免人造血管

图 1不可压缩材料的应用

采用基于位移场的标准有限元方法分析不可压缩材料时存在体积闭锁问题。在此背景下，Sigmund和Bruggi等采用混合位移/静水压力这种两场混合变分公式实现了对可压缩材料和不可压缩材料的拓扑优化设计。Castañar等人提出了一种基于拓扑导数概念的拓扑优化算法，该方法不仅能够处理体积闭锁问题，而且可以计算得到更高精度的单元应力。Jang等人利用基于位移的非协调有限元，实现了不可压缩材料的最小柔度拓扑优化设计。此外，混合尺度边界有限元法也可用于求解不可压缩材料的拓扑优化问题。

据作者所知，有限变形假设下近似不可压缩材料的拓扑优化研究仍然有限。三场混合有限元公式（位移/静水压力/体积改变率）也可以用来模拟近似不可压缩特性，此时通过平衡方程可以直接求解更加准确的体积改变率场。对于近似不可压缩超弹性材料，目前也没有公开的拓扑优化开源代码。因此，本文的目的是在有限变形假设下，提出一种近似不可压缩超弹性材料的拓扑优化方法。同时，本文还提供了完整的299行MATLAB代码和详细的说明，供初级研究人员学习和使用。

三、拓扑优化列式

基于三场混合有限元公式，在几何非线性假设下以结构柔度为目标、以体积为约束的优化模型如下所示：

其中，是结构的终末柔度；是单元总数；和分别是最后一个增量步对应的的外力和全局位移向量；，和分别为位移场、静水压力场和体积改变率场对应的全局残余力向量；和分别是优化设计变量和物理密度；是体积分数，它定义为优化后的结构体积与设计域体积的比值；是第个单元的体积；是设计域的体积，当单元的所有边长都等于1时，设计域的体积等于；是体积分数的上限。

本文采用伴随方法推导目标函数对设计变量的灵敏度，使用移动渐近线法求解上述优化问题，采用能量插值方案来规避低密度单元引起的数值收敛问题。

四、数值算例

图 2 悬臂梁有限元模型

悬臂梁的有限元模型如图2所示，结构尺寸为120mm× 30mm×1mm，它采用120×30个六面体单元进行离散。结构的左侧全约束，右侧中间区域施加载荷F。

不同载荷F（即0.01N、0.05N、0.1N、0.15N）下的优化结果如图3-图6所示。该结构是近似不可压缩的，体积模量设置为1000MPa。当载荷较小时，优化后的结构在上下方向上近似对称。随着载荷大小的增加，优化后的结构失去了对称性。当载荷大小为0.1N时，变形后结构右侧用于传递载荷的杆是垂直的。当载荷水平为0.15N时，优化结构的右侧也是一根杆，且该杆的长度大于载荷水平为0.1N时杆的长度。从这些结果可以看出，载荷的大小对优化结构的材料分布有显著影响。

(a)未变形结构

(b)变形结构

图 3F = 0.01N时的优化结果

(a)未变形结构

(b)变形结构

图 4F = 0.05N时的优化结果

(a)未变形结构

(b)变形结构

图 5F = 0.1N时的优化结果

(a)未变形结构

(b)变形结构

图 6F = 0.15N时的优化结果

参考文献：

Huiqiang Guo, Xinyu Xie, Fei Cheng, Zhengguang Li, Ran Zhang, Jiantao Bai* & 左文杰*. 299-line topology optimization code of nearly incompressible hyperelastic materials using three-field mixed finite element. Computational Mechanics, 2026, DOI: 10.1007/s00466-025-02741-y.

文章主页：

https://link.springer.com/article/10.1007/s00466-025-02741-y #citeas

来源：固体结构CAE工业软件开发

接触有限元分析及其开源Matlab代码

一、摘要接触现象广泛存在于工程中，它是一个高度非线性问题。然而，大多数开源接触有限元代码是用C++编写的，研究人员难以理解和使用。因此，本文提供了完整的摩擦接触有限元方法的详细步骤和Matlab实施代码。本中介绍了接触投影方法、接触节点力和接触切线刚度矩阵的形式，非线性方程组使用牛顿-拉夫逊法求解。数值算例将计算结果与开源软件FEBIO进行比较，验证了Matlab程序的正确性和有效性。吉林大学左文杰教授团队撰写了“An open source MATLAB solver for contact finite element analysis”论文，并发表于Advances In Engineering Software期刊。二、摩擦接触有限元方法图 1 两个弹性接触体对于两个弹性接触体，如图 1所示，粘结和滑移状态的接触牵引力向量可写为：其中，为罚函数法的罚因子；为上时刻主面接触点在当前时刻的坐标；为从面接触点坐标；为接触穿透量；为从面外法线方向向量；为摩擦系数；为摩擦力方向。利用虚功原理，接触虚功可写为：其中，和分别为从面和主面上的虚位移。将上式进行线性化和离散化可得：其中，为线性化算子；为节点虚位移向量；为接触刚度阵；为节点位移。随后，基于牛顿-拉夫逊法进行非线性有限元求解，如下式所示：其中，k为牛顿-拉夫逊法的迭代次数；为结构刚度阵；为位移增量；为接触节点力向量；为载荷向量；为内力向量。上式不断迭代求解直至收敛，即可得到计算结果。三、数值算例第一个算例为两梁接触算例，工况如图2所示，两根梁左端全约束，上梁右侧施加向下的力。上梁下表面和下梁上表面为潜在接触面。图2 两梁接触工况当摩擦系数为0.6时，本软件接触计算位移云图如图3所示。图3 摩擦系数为0.6时的梁合位移云图本软件与FEBIO软件的对比结果如表1所示。可看出，本软件计算结果与FEBIO软件相同，验证了算法和代码的正确性。表1 MATLAB和FEBIO软件的梁合位移结果对比第二个算例为半圆环接触算例，工况如图4所示。半圆环的左侧全约束，右侧施加外力。圆环下表面为潜在接触面。图4 半圆环接触算例摩擦系数设置为0.6，算例合位移结果如5所示。图5 摩擦系数为0.6时的半圆环合位移云图本软件与FEBIO软件的对比结果如表2所示，误差低于1%，验证了算法与代码的正确性。表2 MATLAB和FEBIO软件的半圆环合位移结果对比参考文献：Wang, Bin, Jiantao Bai, Shanbin Lu, and Wenjie Zuo*. An open source MATLAB solver for contact finite element analysis. Advances In Engineering Software, 2025, 199: 103798.文章主页：https://www.sciencedirect.com/science/article/pii/S0965997824002059接触有限元Matlab代码下载地址：https://www.researchgate.net/publication/386332459_ContactFEA_codesrar来源：结构设计CAE工业软件研发