有限元理论基础及Abaqus内部实现方式研究系列47：约束关系（3）-船舶规范约束导致的Max Ratio问题

SnowWave02 2024年8月6日浏览：1659 评论：3 收藏：3

（原创，转载请注明出处）

1 概述

本系列文章研究成熟的有限元理论基础及在商用有限元软件的实现方式，通过

（1）基础理论

（2）商软操作

（3）自编程序

三者结合的方式将复杂繁琐的结构有限元理论通过简单直观的方式展现出来，同时深层次的学习有限元理论和商业软件的内部实现原理。

有限元的理论发展了几十年已经相当成熟，商用有限元软件同样也是采用这些成熟的有限元理论，只是在实际应用过程中，商用CAE软件在传统的理论基础上会做相应的修正以解决工程中遇到的不同问题，且各家软件的修正方法都不一样，每个主流商用软件手册中都会注明各个单元的理论采用了哪种理论公式，但都只是提一下用什么方法修正，很多没有具体的实现公式。商用软件对外就是一个黑盒子，除了开发人员，使用人员只能在黑盒子外猜测内部实现方式。

有限元理论基础及Abaqus内部实现方式研究系列47：约束关系（3）-船舶规范约束导致的Max Ratio问题的图1

一方面我们查阅各个主流商用软件的理论手册并通过进行大量的资料查阅猜测内部修正方法，另一方面我们自己编程实现结构有限元软件iSolver，通过自研CAE软件和商软的结果比较来验证我们的猜测，如同管中窥豹一般来研究的修正方法，从而猜测商用有限元软件的内部计算方法。我们关注CAE中的结构有限元，所以主要选择了商用结构有限元软件中文档相对较完备的Abaqus来研究内部实现方式，同时对某些问题也会涉及其它的Nastran/Ansys等商软。为了理解方便有很多问题在数学上其实并不严谨，同时由于水平有限可能有许多的理论错误，欢迎交流讨论，也期待有更多的合作机会。iSolver包括完整的前后处理和有限元求解器，功能如下，有兴趣可直接在下面网址下载：

百度网盘链接: https://pan.baidu.com/s/10d6jHdZ01SBY2JxiS6bffw 提取码: 6fdf

有限元理论基础及Abaqus内部实现方式研究系列47：约束关系（3）-船舶规范约束导致的Max Ratio问题的图2

2 约束关系

我们在系列文章第35章介绍了接触分析及常用的基本算法。现在Abaqus、LS-DYNA、Ansys等结构商软都说可以处理复杂的上万零部件接触的整车、整机等模型仿真，没做过实际的这种仿真分析，很好奇，接触分析算法往往涉及大变形、边界不连续，只要输入条件或者算法稍微变化一些，两个零部件算出来的接触结果就可能差异很大，更不用说上万个零部件的接触结果了，对这种大规模组装模型的仿真结果不知如何来判断它的可靠性，像普通的只校核一下材料的应力还是看一下动画是否和试验一致？毕竟仿真只有简单的标准来判断结果的正确性才能在企业中起到真正辅助设计的作用，如果你恰好做过，不妨也简单介绍一下你的经验。对自研CAE软件开发者来说，自研结构CAE软件是否真的要和商软去比拼接触等复杂算法还是花更多时间在精雕细琢那些常用功能上，这也是开发者需要慎重考虑的问题，而且很多自主CAE软件连常规线性问题都算的不对，或者都没法用鼠标稳定的走完那些材料、属性、边界、加载等流程，用户又怎么会相信你能算对接触这种复杂问题的？

有限元理论基础及Abaqus内部实现方式研究系列47：约束关系（3）-船舶规范约束导致的Max Ratio问题的图3

不管怎么样，从有限元实现的角度来讲，如果想做真正实际工程中的接触分析，那么首先需要去研究约束关系，接触分析在有限元中也仅是约束关系的一种。有限元中的约束很多场景大家用的是边界中的简支、固支等约束，但从更广泛的角度上讲，只要表示一个节点的某个自由度依赖于其它的节点自由度或者取某个特定值，就可以称为约束关系。只不过对固支、简支等直接自由度=0，在有限元中直接减缩刚度阵就行，很容易求，但对节点自由度相互依赖的约束关系就比较复杂了。约束关系主要有两类。

（1）一类是MPC点之间的约束。Nastran的MPC的灵活度要远远超过Abaqus，Nastran的主节点可以选择123自由度，也可以对每个从节点设置不同的自由度，还能主节点和从节点互相包含，Abaqus更多的是只负责80%的常用应用场景，复杂功能让你编子程序，但事实上一线仿真工程师又有多少人愿意编子程序呢？这种做法导致虽然Abaqus无论从用户体验、非线性还是商业化都比Nastran好很多，但很多线性的工程复杂问题还是没法替代Nastran。

（2）另一类是Contact、Tie等的面之间的约束关系。在这方面Abaqus要明显强于Nastran了。

我们将用统一的公式来求解这两类关系，同时也从软件实现层面说明一下针对这两类情况的各自差异。分几篇文章来介绍约束关系，本篇是约束关系（3）- 船舶规范约束导致的Max Ratio问题，这是我们碰到的1个实际的工程问题，当自主CAE软件往外推广时，只要用，就会有各式各样的问题，最基本也是最重要的一条是自主CAE软件算出来的结果只要不符合预期或者商软的结果，就必须要你解释why？不会有人觉得商软或者建模等等有问题，无一例外都默认是自主软件的错。不过这也正常，一开始商软推广也是这么过来的，就是现在，如果商软客户提出问题，一般商软技术支持的响应速度也是必须要在24个小时内回复。很多时候用户用自主软件算一个大模型告诉你计算不正确，查到最后都是某个单元、某个载荷、某个边界的单独或者组合的问题，但从大模型找到这个小问题往往需要调试很久。在本章我们仅举我们调试过的一个最基本的MPC导致的MAX RATIO工程问题，如果你经常用商业CAE软件查看后台计算打出的消息，这个MAX RATIO可能就会经常碰到，导致MAX RATIO出错的情况很多，在这里仅介绍我们实际碰到的由MPC约束不够导致的一种情况。

3 船舶规范约束导致的Max Ratio问题

3.1 模型建模规范

很多行业都有自身的建模规范，在船舶行业，油、散、集三种船型一直是民船市场的三大主流船型，从下面2023年新接单可知油船和散货船占据前两位。而油船和散货船的结构强度校核必须符合船舶行业的共同结构规范(Common Structure Rule，CSR)，该规范在2006年由国际船级社协会(IACS)颁布，一直沿用到现在。

有限元理论基础及Abaqus内部实现方式研究系列47：约束关系（3）-船舶规范约束导致的Max Ratio问题的图4

为配合规范的发布，IACS也同步发布了钢质海船入级规范，所有需要入级的船舶结构有限元的建模都必须符合这个规范。其中，舱段结构有限元分析一般只建中间三个舱段，而两边剩余舱段对这个三舱段的约束按规范加载。

有限元理论基础及Abaqus内部实现方式研究系列47：约束关系（3）-船舶规范约束导致的Max Ratio问题的图5

在船舶CSR规范直接计算章节规定了非船艏舱段有限元建模时两个端面的边界条件如下：

有限元理论基础及Abaqus内部实现方式研究系列47：约束关系（3）-船舶规范约束导致的Max Ratio问题的图6

对船舶规范不是很了解，只是按规范文档简单说明一下：由于舱段前后端面类似，我们以前端面为例，中和轴与纵中剖面相交处建一个独立点（MPC的主节点），该独立点仅是三维空间内的一个参考点，这个点按规范约束U2、U3和UR1。

有限元理论基础及Abaqus内部实现方式研究系列47：约束关系（3）-船舶规范约束导致的Max Ratio问题的图7

同时，前端面的横剖面内所有节点（MPC从节点）和该独立点用刚性连接RBE2绑定U2、U3和UR1自由度。

有限元理论基础及Abaqus内部实现方式研究系列47：约束关系（3）-船舶规范约束导致的Max Ratio问题的图8

3.2 反馈的工程问题

按这个绑定，理论上讲该横剖面y和z方向不会平动，且不会绕x轴（船长方向）转动，而可以在x方向平动，且可以绕y和z轴转动，我们对规范不了解，但和我们主观上其它舱段对中间三个舱段的作用还是比较吻合的。那么正常来说，有限元计算出的结果也应该是这样。但事实上用户对不同的三舱段模型，采用iSolver或者Nastran计算，部分结果正确，部分模型结果不对。

3.2.1 1号模型正确结果

譬如同样都是重力载荷下，1号模型计算出的位移为43.64，且端面位移明显很小，相对合理：

有限元理论基础及Abaqus内部实现方式研究系列47：约束关系（3）-船舶规范约束导致的Max Ratio问题的图9

3.2.2 2号模型错误结果

但2号模型计算出的位移为3.602e13，且后端z方向平动也非常大，这明显不合理。

有限元理论基础及Abaqus内部实现方式研究系列47：约束关系（3）-船舶规范约束导致的Max Ratio问题的图10

而且用户反馈采用Nastran计算，这两个模型也是同样的结果，Nastran对2号模型在.f06文件报告PIVOT RATIO过大的致命错误，Nastran计算失败不再输出结果。

这就需要回答两个问题：

（1）对同样的约束，为何有限元软件有些模型计算正常，有些模型计算失败

（2）如果连Nastran也存在这个问题，难不成规范有问题或者Nastran又是怎么来处理这个问题呢？

3.3 原理分析和猜测

3.3.1 约束关系的原理分析

有限元理论基础及Abaqus内部实现方式研究系列47：约束关系（3）-船舶规范约束导致的Max Ratio问题的图11

即两个主节点的UR2，UR3与船体端面的节点的位移无关。

同时，主节点只有U2，U3、UR1约束。导致主节点（即独立点）剩下的UR2，UR3自由度没有约束或者绑定在舱段，因为主节点这两个自由度度没有和任何其它节点组成单元或者约束关系，按我们前面的讨论，刚度阵都是两个节点之间的关系，所以导致刚度阵UR2和UR3相关的列或者行都是0，一个矩阵只要有一列或者一排都是0，那么理论上刚度阵K就无法求逆，所以正常来说，按船舶CSR规范的约束理论上都会导致两个MPC主节点的UR2和UR3都约束不够。

3.3.2 实际编程的分析

实际编程和理论相比存在差异：

3.3.2.1 模型或者精度误差

实际编程不会直接按上面少了简化自由度的（2）式来编写，依然按全部自由度的（1）来编写，这样UR3前的rx1依然要写到程序代码中，由于创建模型、计算机存储浮点数存在截断误差或者rx1可能也是某些向量叉乘得到的误差，那么rx1不一定是绝对的0，有限元求解器不一定会报这个错误了。

譬如上面正确的1号模型，查看前端面的所有节点和独立节点的x坐标，正常来说都应是121500这个x坐标，但发现70648号和72048号节点却是121504，和主节点差了一个小量，相对误差0.003%。这个小偏差导致主节点的UR2和UR3和舱段的U2、U3有约束关系，间接的绑定了这个主节点的UR2和UR3自由度。后端面也有两个类似存在小误差的节点。

有限元理论基础及Abaqus内部实现方式研究系列47：约束关系（3）-船舶规范约束导致的Max Ratio问题的图12

将这四个点的X值人为修改为所在MPC同一平面后，iSolver计算也出现了2号模型的错误结果，此时明显位移极大不正常。同时，Nastran也是计算失败，报告PIVOT RATIOS过大的致命错误。

有限元理论基础及Abaqus内部实现方式研究系列47：约束关系（3）-船舶规范约束导致的Max Ratio问题的图13

3.3.2.2 高斯消元法求解问题

第二个问题是结构有限元一般采用直接法来求解矩阵，高斯消元法最常见的直接法，高斯消元法将矩阵分解成单位下三角矩阵L和上三角矩阵U，高斯消元法不是直接求K的逆矩阵来求解。刚度阵每一行都对应某个自由度，这一行的刚度*未知量=这个自由度上的载荷（由文章第45篇可知，其实也是一个约束方程），如果这个行的刚度都是0，只要这个自由度对应的载荷项也是0，那么在高斯分解成LU矩阵后反向求解自由度对应的位移时，可以跳过这个方程，同时将这个自由度对应的位移设为为0就行。

3.3.3 商软的做法

3.3.3.1 PIVOT RATIO的定义

无论上面两种方法中的由于截断误差导致的rx1不为0还是第二种刚度阵某一行都是0的情况，在实际工程中都必须判断一个数是不是和0的关系，计算机也只能严格按照一个定死的判据判断这个数是否为0。由于刚度在不同单位不同材料情况下绝对值差异很大，所以商软都采用相对量来判断近似0的情况，实际中Nastran或者Abaqus都采用高斯消元法的变形，将刚度阵直接LDLT分解为下三角L、对角阵D、上三角阵L’。

K=L*D*L’

它们采用Pivot Ratio来说明这个近似0到什么程度，也就是原始的刚度阵对角元Kii和LDLT得到的D矩阵对角元Dii的比值：

Pivot Ratio (i)= K(i,i)/D(i,i)

3.3.3.2 MAX PIVOT RATIO不同软件的修正

无论是Nastran还是Abaqus，PIVOT RATIO都是上面的定义，但我们研究发现，Nastran和Abaqus还是略微有些不同的，Nastran做了进一步的修正。

我们做了一个两头按船舶CSR规范加载的盒子3号模型。这个模型的两个端面和主节点严格在一个平面，且模型相对简单。

有限元理论基础及Abaqus内部实现方式研究系列47：约束关系（3）-船舶规范约束导致的Max Ratio问题的图14

iSolver、Abaqus、Nastran计算都可发现能得到正常结果。

有限元理论基础及Abaqus内部实现方式研究系列47：约束关系（3）-船舶规范约束导致的Max Ratio问题的图15

Abaqus计算提示MPC主点的1/5/6自由度没有约束，它的PIVOT RATIO为3.7e14。

有限元理论基础及Abaqus内部实现方式研究系列47：约束关系（3）-船舶规范约束导致的Max Ratio问题的图16

Nastran计算不再提示MPC主节点的1/5/6自由度没有约束，Nastran PIVOT RATIO的阈值为1e7，如果超过这个阈值会报错的。修改bdf，加入PARAM MAXRATIO 1.E2：

有限元理论基础及Abaqus内部实现方式研究系列47：约束关系（3）-船舶规范约束导致的Max Ratio问题的图17

发现Nastran只会报节点5的Max Ratio错误，反推MPC两主节点64、65的Max Ratio算出来还比正常节点5的3.4e2更小，猜测Nastran在计算MAX PIVOT RATIO时，如果Kii=0，那么直接将MAX PIVOT RATIO设为了0，但Abaqus估计不是这样的。

有限元理论基础及Abaqus内部实现方式研究系列47：约束关系（3）-船舶规范约束导致的Max Ratio问题的图18

3.3.3.3 AUTOSPC的作用

PIVOT RATIO超过阈值的情况是非常常见的，譬如这个船舶CSR规范正常来说都会使得PIVOT RATIO超过阈值，所以一般商软会尝试自动处理这种情况。

Nastran默认情况下采用自动约束处理(具体在bdf文件中设置autospc关键词=yes)在处理完 1）从属自由度(from MPC entries or RIGID elements)，2）用户指定的拘束自由度(from SPC/SPC1 entries)后，在进行高斯消元的LDLT分解前或分解过程中，对明显的奇异节点自由度进行自动约束处理。

很奇怪，Nastran的AUTOSPC后台没有直接将MAX RATIO和阈值比较，而是将原始刚度阵的对角元Kii和该节点类似最大主刚度（采用3X3的特征值求解）的比值Ri和EPZERO比较，如果Ri小于EPZERO，那么自动约束该自由度，否则不自动约束，这种求特征值和Ri的方式显然要比MAX RATIO更复杂。

因为Nastran如果bdf没有特别设置，那么默认采用AUTOSPC=YES，且一般情况Ri总是小于EPZERO，使得船舶CSR规范的约束一般都能在Nastran正确计算。

3.3.3.4 AUTOSPC失败时怎么办

上面只能说一般情况下设置AUTOSPC能解决这种规范约束不够的问题，但当测试模型足够多后，我们发现有些模型即使Nastran设置了AUTOSPC=yes，依然会报告UR2或者UR3超过MAX RATIO导致Nastran计算失败。

譬如：4号模型为一个工程模型，没加AUTOSPC时，f06报Max Ratio错误，注意此时的PIVOT RATIO是-2.0e15。

有限元理论基础及Abaqus内部实现方式研究系列47：约束关系（3）-船舶规范约束导致的Max Ratio问题的图19

加入AUTOSPC，发现依然为MAXRATIO错误，且RATIO值还不一样了，从-2.0e15变为了-1.5e15，同样也不清楚为何。

有限元理论基础及Abaqus内部实现方式研究系列47：约束关系（3）-船舶规范约束导致的Max Ratio问题的图20

在该例子中，应该是得到的Ri>EPZERO导致没有自动约束。按照Nastran文档，EPZERO默认为1e-8，可以修改采用PARAM, EPZERO关键词修改默认值。将bdf的PARAM, EPZERO修改为3e-7。

有限元理论基础及Abaqus内部实现方式研究系列47：约束关系（3）-船舶规范约束导致的Max Ratio问题的图21

可发现结果正常。查看f06，可发现GPST中有4个自由度Ri=0，而另两个是2.3e-7，此时Ri<EPZERO(3e-7)，自动约束。

有限元理论基础及Abaqus内部实现方式研究系列47：约束关系（3）-船舶规范约束导致的Max Ratio问题的图22

3.4 总结

本章是一个实际MPC约束不够导致的MAXRATIO工程问题，我们发现同样按照船舶CSR规范创建的模型计算部分成功，部分失败，可能原因有两个：

（1）模型中MPC主从节点不是严格在同一个X平面，有部分小的偏差，这个小偏差导致主节点的UR2和UR3和舱段的U2、U3有约束关系，间接的绑定了这个主节点的UR2和UR3自由度。

（2）即使MPC主从节点依然在同一个X平面，实际计算机存在截断误差或者计算精度误差，导致部分数据不一定是绝对的0，需要和近似0的阈值比较判断，部分能通过判断，部分不能通过。

同时，我们提出了一些MPC原理分析和商软的猜测，也遇到了一些问题，也希望有大神能解答一下：

（1） Nastran采用PARAM，MAXRATIO来设置MAX RATIO的阈值只能针对正数，但Nastran打出的PIVOT RATIO有时还有负数，这时有类似的负数的阈值设置吗？我们没有找到。

（2）如果程序来实现自动约束的功能，按我们的想法，应该直接判断MAX RATIO是否超过阈值，超过就直接自动约束，但实际Nastran的AUTOSPC自动约束重新计算了一个类似MAX RATIO的Ri参数，这个Ri还与特征值有关，为什么Nastran明明简单的MAX RATIO，而要舍近求远的再求一个更复杂的Ri呢？

（3）在此我们只讨论Nastran针对这个问题的做法，Nastran的关键词更加灵活，且bdf各种参数说明文档相对更详细，Abaqus暂时没有研究过，如哪位大神也有类似研究，欢迎不吝赐教，我们也很想知道Abaqus后台矩阵求解的逻辑和针对高斯消元法的实际工程处理。

（4）不了解船舶专业规范，但仅从原理来讲，如果报错，一种修改方式是主节点直接都约束UR2和UR3，这种方式以前能算的的结果不影响，但报UR2和UR3没约束的就可以消除UR2和UR3，得到正常的结果。从一个外行的角度觉得船舶专业规范太依赖Nastran的AUTOSPC功能了，上面也证明了Nastran的AUTOSPC有时也可能导致结果不正确，实际上可以不用管后台到底是哪个软件，而在规范中规定两个端面主节点UR2和UR3也强制约束。

这些猜测具体什么原因最好的方式还是调试矩阵求解的代码。只可惜我们没有自己编写矩阵求解，为了追求稳定和高效，我们选择了现成的成熟的矩阵解算库。很希望有哪位大神通过调试底层矩阵求解来说明和实际有限元软件报的某个自由度PIVOT RATIOS过大的关系，然后能确切的告诉我们如何设置矩阵解算库类似Nastran的autospc功能。

同时，对约束导致的MAX PIVOT RATIO和AUTOSPC研究的越深，一方面感叹商软在这个问题上处理的工程修正和精细化程度相当的深入，另一方面也让我们心生敬意，觉得比起商软来说，自主软件无论是底层矩阵求解方面还是高层的工程应用真的还有很长一段路要走，而且很多东西是系统化的耦合在一起的，只有从底层和应用层一起努力才可能达到商软的水平。