基于Carrera统一公式的任意多边形网格自适应有限元评估

　　

　　

　　网格技术代表了离散感兴趣的领域的能力，以最好的方式拟合真实的物理连续体。最常用的方法是有限元法(FEM)，一种求解偏微分方程的数值方法。为了克服有限元法提出的经典问题，研究了其他模型。目标是允许用任意多边形表示的元素离散化问题域，这些多边形可以是凹的，也可以是凸的。此外，在这些方法中寻求不同的多项式一致性，并可能处理非一致性离散化，主要是局部细化等。本工作旨在提出新的自适应元素，即基于Carrera统一公式的有限元，以证明这些新元素可以完成所有先前的功能，并且易于实现相关模型。首先，通过经典的贴片测试来研究网格畸变的敏感性。然后，给出了不同的研究案例，这些案例包含了非常扭曲的凹、凸单元。

　　网格技术代表了离散感兴趣的领域的能力，以最好的方式拟合真实的物理连续体。在工程领域，网格划分是建立待解决问题的模型的一个重要里程碑，它是进行仿真的起点，重要的是它能够更准确地表示所分析的领域。一般来说，所调查的物理区域不那么规则，网格划分应该考虑到这一点。

　　最常用的技术是有限元法(FEM)，它依赖于一个非常简单的数学模型，能够描述所调查的领域(结构和/或流体动力学)达到一定的精度。该方法允许构建2D和3D规则元素，每个元素具有规定数量的节点[1]。该方法的主要问题表现为:对元素规则预设的约束;管理同一域中不同类型元素的困难，可能有不同数量的节点，因此，有不同的近似;由于使用的元素数量多，计算成本高，当实际域不是几何规则或当需要丰富网格时，因为它在感兴趣的区域。

　　为了减轻有限元法的网格划分问题，可以使用更先进的技术，采用从一开始就设计的数值方法，在更一般的形状元素上提供任意顺序的精度。例如，这些技术是基于虚元法(VEM)或多面元法(PEM)[2]。事实上，与FEM中二维单元通常是三角形和四边形，三维单元通常是四面体和六面体不同，VEM允许任意二维多边形和三维多面体单元[3,4]，其许多方面与FVM密切相关，特别是与模拟有限差分(MFD)[5]。这允许用任意多边形表示的元素来离散问题域，这些多边形可以是凹的，也可以是凸的。此外，方法内允许有不同的多项式一致性，可以处理不一致的离散化，主要是局部细化等，这是该方法的关键方面[6]。

　　新提出的自适应单元是基于Carrera统一公式(CUF)的有限元，它代表了更灵活和可管理的单元，基于一种众所周知的技术(FEM)，易于实现，并且可以产生与使用VEM相同的效益。事实上，CUF既适用于一维有限元，也适用于二维有限元，可以提供非常精确的三维解，且计算成本远低于经典的三维单元，成为分析结构问题的重要工具。

　　最近，为了进一步提高模型的数值效率，在CUF框架下，提出了一种新的方法，称为节点依赖运动学(NDK)[7,8,9,10,11]。通过将运动学假设与所选节点联系起来，可以方便地建立可变节点运动学的有限元模型。当然，采用高阶假设的临界区域可以与采用充分的低阶理论建模的非临界区域相衔接[12]，这样可以进一步简化数值模拟过程。

　　NDK方法与非正交几何的3D建模一起，为创建具有非常规节点数量、非常规几何形状和沿不同方向的不同多项式近似的元素提供了可能性，从而在构建适当的元素时具有一定的自由度。考虑到控制方程(刚度和质量矩阵)矩阵内部的积分，其中涉及近似函数，以及元素的雅可比矩阵以3D形式计算，尽管仍然使用1D和2D CUF模型，但这些元素的形式是公式化的。总之，当几何约束或运动学行为要求时，3D元素可以恢复为不规则形状，可以适应所考虑的域的边缘或共享一个或多个边的相邻部分。

　　本文组织如下:第一节介绍;第2节简要回顾了基于CUF的1D和2D模型;第3节介绍了通过使用有限元方法和CUF从虚拟位移原理推导控制方程的过程，其中也在NDK技术的框架中解释了目前的方法以及3D建模;第4节包含通过对一些有用的研究案例进行静态和模态分析获得的结果，这些研究案例展示了当前元件的能力;最后，第五节总结了本文的结论。

　　CUF(参见[13,14])以统一的方式描述运动场，然后利用该方式以紧凑的方式推导控制方程，这与VEM不同，易于实现。与经典理论相比，这种紧凑的形式具有巨大的优势，经典理论有一些局限性(正如之前许多与CUF相关的作品所解释的那样)。这些限制的主要是，运动学应该在理论上丰富了无限数量的术语(见Washizu[15])，以克服物理不一致性并预测高阶效应。在研究实际问题时，这是不可行的。有必要截断主要力学变量的展开，沿着结构域的最小维度，到一个给定的顺序，这通常是在制定结构理论时做的。此外，运动学中的项越多，问题的表述和求解就越复杂。CUF代表了这些问题的解决方案，以紧凑和统一的形式描述控制方程，允许制定1D和2D模型(分别是梁和板模型)。

　　让我们考虑一个一般的梁结构，它的纵轴，相对于笛卡尔坐标系，位于坐标y上，作为它在xz平面上定义的横截面。CUF框架下一维模型的位移场被描述为广义位移(在基于位移理论的情况下)通过横截面坐标的任意函数的一般展开:

　　（1）

　　式中为三维位移矢量，为一般位移矢量，M为展开式中的项数，表示求和，函数定义要使用的一维模型。事实上，根据函数的选择，可以实现不同类型的梁理论，如Euler-Bernoulli或Timoshenko梁模型等等。

　　在板理论的框架下，考虑板的中平面位于xy平面内，CUF可以用类似的方式表示(见[16]):

　　（2）

　　在上式中，广义位移是板的平面中坐标的函数，在厚度方向z上展开。同样，根据展开函数的选择，可以形成各种理论，如Kirchoff-Love或Reissner-Mindlin板模型等。本文只考虑拉格朗日展开，但如Carrera等人[8,10,14]所讨论的，可以方便地使用基于Taylor或Legendre多项式的其他类型展开[8,10,14]。

　　拉格朗日展开理论是基于使用拉格朗日型多项式作为跨截面或沿厚度的一般函数。

　　对于梁模型，截面被划分为若干局部展开子域，其多项式度取决于所采用的拉格朗日展开类型。四节点双线性Q4、九节点三次Q9和十六节点四次Q16多项式可以用来表述精细的梁理论(见Carrera和Petrolo[17])。例如，Q9展开式的插值函数定义为

　　（3）

　　其中和在和之间的横截面上变化，和和表示根在自然平面上的位置。因此，单q9梁理论的运动场为

　　（4）

　　在多域截面上采用高阶拉格朗日多项式或拉格朗日多项式的组合可得到精细化的梁模型。关于拉格朗日类梁模型的更多细节可以在[17,18,19,20]中找到。

　　基于拉格朗日展开的板理论，采用一维拉格朗日多项式作为厚度函数，最高阶为p:

　　（5）

　　其中I表示求和，取值范围是2到。这样就可以使用二节点线性多项式B2、三节点抛物多项式B3和四节点三次多项式B4。

　　在二维CUF理论的框架内，类似上述的位移场已被用于实现层合结构的分层模型(例如，参见[14,21,22])。

　　拉格朗日多项式被Kulikov和他的同事广泛地应用于统一框架下的变运动学板壳理论的表述中。关于基于拉格朗日型展开的二维模型的更多细节，请参考原文[23,24]。

　　CUF的主要优点是它允许以紧凑和统一的方式编写控制方程和相关的有限元数组，这在形式上是函数的不变量。在下面的章节中，详细介绍了CUF一维和二维模型下刚度矩阵基本核(不变量)的数学推导。

　　在本节中，采用与第2节中定义的相同的符号和参考系统。设三维位移向量定义为

　　（6）

　　根据经典弹性理论，应力张量和应变张量可以组织成六项向量，具有一定的通用性。他们分别写道:

　　（7）

　　根据这个表达式，应变与位移之间的几何关系可以用紧向量符号来定义为

　　（8）

　　式中，在变形和转角较小的情况下，为线性微分算子:

　　（9）

　　另一方面，对于各向同性材料，应力与应变之间的关系可以通过著名的胡克定律得到:

　　（10）

　　其中C是各向同性刚度矩阵

　　（11）

　　刚度矩阵的系数仅取决于杨氏模量E和泊松比，在各向同性材料的情况下，它们为:

　　（12）

　　对于一维模型，采用有限元方法对梁的纵轴进行离散。将广义位移描述为未知节点矢量的函数，以及一维形状函数:

　　（13）

　　每个元素的节点数是多少，未知节点向量定义为

　　（14）

　　同样，板中表面广义位移的有限元离散化可表示为:

　　（15）

　　其中使用二维形状函数。

　　不同的多项式集合可以用来定义有限元单元。本文选择拉格朗日插值多项式来生成一维和二维元素。为简洁起见，本文没有给出它们的表达式，但在Carrera等人[13]的书中，对两节点(B2)、三节点(B3)、四节点(B4)梁单元和四节点(Q4)、九节点(Q9)、16节点(Q16)板单元进行了描述。

　　将有限元逼近与CUF的运动学假设相结合，三维位移场可表示为

　　（16）

　　其中函数和根据单元类型(梁或板)定义。

　　在Eq. 16中，形函数和展开函数是独立的。最近，Carrera等人[13,25,26]通过以下形式将展开函数与形状函数联系起来，引入了一种耦合:

　　(17)

　　方程17和方程16的不同之处在于附加的上标i，它现在也是函数的一个指标。该定义引入了运动学假设对有限元节点的依赖性，即NDK。

　　如前所述，增加阶数的扩展函数为更好地逼近结构响应提供了机会。利用NDK可以在特定节点上对运动学模型进行局部细化，便于进行局部自适应的运动学细化。不同的结构理论可以通过单元域内的节点型函数自然地融合在一起，而无需任何特殊的耦合方法。同时，不需要对不同的节点运动学进行兼容性要求。

　　NDK技术可以有效地应用于结构的全局-局部建模。在临界区域的运动学模型可以细化，直到达到理想的精度，而留下的外围区域与适当的低阶理论建模。该方法可以在不改变网格的情况下方便地构建全局局部有限元模型，并且可以重复使用同一组网格网格来构建一组模型，用于并发全局局部分析。该方法已在一维[25,26]和二维[13,17,27]情况下用于层压结构的高效建模。

　　对于变截面或断面与轴线不正交的梁，以及变厚度或边缘与中表面不正交的板，经典的梁板理论在建模中不能得到明确的应用。在NDK方法的框架下，通过将CUF运动学假设和FEM离散化结合在一个独特的三维近似中，可以将CUF模型扩展到这种几何形状的建模，如下所示:

　　(18）

　　其中为x, y, z对应的自然坐标。在此表达式中，表示一个非传统的三维形状函数，其中沿某个空间方向的展开顺序可以不同。类似地，位移的虚变化，将用于下面的控制方程的推导，可以近似为

　　（19）

　　其中引入了新的求和索引s和j。

　　考虑到这种形式，控制方程中涉及的体积积分将不会像通常那样分成一维和二维积分，而是将函数和作为规则的三维形状函数处理;因此，相对于从自然坐标到全局坐标x, y, z的变换，雅可比矩阵将以三维形式计算。

　　结果将表明，这种特殊的方法使我们能够节省自由度，而不是使用基于经典3D有限元的网格。事实上，传统的3D元素在三个空间方向上使用相同的扩展顺序，这通常意味着元素长宽比的选择受到一些限制，导致使用的元素数量有害增加。

　　利用虚位移原理得到了控制方程。这个变分表述设定了一个结构平衡的必要条件，即内部功的虚变等于外部功的虚变，或者

　　（20）

　　内功等于弹性应变能:

　　（21）

　　其中V代表定义域的体积。采用几何关系[Eq.(8)]、本构律[Eq.(10)]、CUF运动场[Eq.(1)]和有限元离散化[Eq.(13)]，内功可改写为

　　（22)

　　其中为刚度矩阵的基本核。

　　惯性功可以写成

　　(23)

　　采用CUF运动场和有限元离散化(方程)。(18)和(19))，则变为:

　　(24)

　　其中是质量矩阵的基本核。

　　外部功可以写成

　　（25）

　　按照对内部做功的相同步骤，也可以将外部做功用节点向量表示为

　　（26）

　　其中，载荷矢量以形状函数和展开函数表示。

　　应该注意的是，基本核组成部分的形式表达式与扩展函数的选择无关，扩展函数决定了结构理论，而形状函数决定了有限元近似的数值精度。这意味着任何梁或板单元都可以根据指标s、i和j适当展开基本核，从而自动形成。关于刚度矩阵、质量矩阵和载荷矩阵基本核的详细表达式，可参考著作[28]。

　　本文分别使用了线性、二次和三次拉格朗日型形状函数。请注意，缩写词B2(线性)，B3(二次)和B4(三次)既用于梁模型中轴的一维离散化，也用于板模型中通过厚度的近似;以及缩写词Q4(线性)，Q9(二次)和Q16(三次)都用于梁模型中截面的二维近似和板模型中中表面的离散化。根据这些特点，设置展开阶数作为模型的自由输入，决定待解未知量的个数。同时，通过多域截面的定义，可以在特定的兴趣区域丰富运动场。

　　摘要。

　　1 介绍

　　2 卡雷拉统一公式:一维和二维模型

　　3 有限元公式

　　4 数值结果

　　5 结论

　　参考文献。

　　作者信息

　　# # # # #

　　在本节中，考虑了一些基准问题来评估基于CUF的自适应元素的鲁棒性。首先，一组涉及均匀和各向同性板的问题被用于评估元件的性能，如4.1节所示:考虑畸变元件形状的位移敏感性和自由振动分析。本文还讨论了边界条件的影响。

　　然后，对具有变形单元的平板进行了静力和模态分析。扭曲的网格是从规则模式开始以随机方式创建的:节点的新位置使用一个函数获得，该函数从正态分布或高斯分布中生成一个随机标量，其平均值等于先前的位置，并且具有不同的标准差。

　　对于所有列出的问题，都考虑了四节点和九节点元素。

　　该标准试验通常用于研究板弯曲问题中的网格变形敏感性。所考虑的方形板的平面内尺寸为，厚度为。所应用的边界条件为:clamped-clamped (CC)，所有边都是固定的;简支(SS)，所有边缘在厚度方向和平行于边缘的方向上没有位移。荷载表示为顶面中心的集中力等于。所考虑的材料是各向同性的，具有以下特性:，和。起始网格是平面上的Q4元素和Q9元素，在这两种情况下沿厚度都有B3近似。

　　由于问题的对称性，在几何对称和边界条件对称的情况下，为了简化分析，考虑四分之一的板，将整个板分成四个相等的子板，只对其中一个子板进行研究，如图1所示。其中，选择左下角四分之一，并应用以下边界条件，分别考虑u, v和w沿x, y和z方向的位移:底部和左侧边缘固定或简支，在顶部边缘，右侧边缘和右上角施加力。

　　图1

　　

　　中心点横向集中荷载的各向同性板。只有四分之一的板与元件啮合

　　得到的板四分之一网格基于Q4或Q9单元，沿厚度进行B3离散，如图2所示。

　　图2

　　

　　四分之一板原始网格

　　注意，从现在开始，在网格中引入一些梁单元，以获得沿边缘具有不同节点数的扭曲单元。然后，分析了三种不同的配置:

　　2xQ4 1xB2:在这种构型中，顶部的前两个元素被一个具有2个节点的一维元素取代，在截面上具有Q9近似;

　　2xQ4 1xB3:与前一种情况一样，顶部的两个元素被单个1D元素取代，但沿轴有3个节点，其中截面表示为Q9元素;

　　2xQ4 1xQ9:节点数与前一种情况相同，但顶部两个元素的离散化不同，因为单个Q9元素沿厚度具有B3离散化。

　　在图3和图4中，可以将所描述的配置可视化:黑点是实现过程中使用的节点，而白点分别表示考虑板单元和梁单元的厚度和截面而添加的节点。为不同的元素引入了颜色代码:B2元素为天堂色，B3元素为黄色，Q4元素为浅蓝色，Q9元素为橙色。

　　图3

　　

　　四分之一板- 2xq4 1xB2

　　图4

　　

　　四分之一板，两种不同的结构

　　现在的重点将集中在梁元素上，以及由于其近似函数的三维积分而自由移动其节点的可能性，如第3.4节所讨论的。

　　网格畸变用参数来表征，该参数定义了板截面中心节点的坐标。请注意，参数s可以同时具有负值和正值:这在文献中并不常见，特别是在本工作中，它用于控制元素的属性:事实上，如图5所示，如果参数s为负值，则所考虑的元素是凸的;如果参数s为正，则该元素变为非凸。中节点的y坐标根据s变化，定义a为四分之一板的长度，坐标变化为，在图5a中，其中，则坐标为;在图5b中，坐标为。

　　图5

　　

　　由参数s定义的四分之一板网畸变

　　4.2.1 静力与模态分析

　　如文献[29]所述，使用三种不同的配置进行静态和模态分析。在静力分析中，变形网格的计算结果是按横向位移(整个板的中表面中心点)计算的，并根据规则网格计算的参考值进行归一化()。为了验证目的，将不同网格的值与Nastran中计算的3D解决方案进行比较。可以注意到，2x2Q9解决方案是最准确的，自由度最多，2x2Q4是最不准确的，其他解决方案根据其相关的近似程度处于中间。网格2xQ4 1xQ9和2xQ4 1xB3非常相似，它们被比较以突出它们的解决方案完全相同，尽管第一个使用板元素(Q9)和另一个使用梁元素(B3)。此外，可以看到2xQ4 1xB2解决方案比2x2Q4解决方案更准确，尽管自由度的数量是相同的:这是由于在B2梁单元的横截面上使用二次逼近。我们可以得出结论，目前的自适应元素允许提高精度相对于具有相同自由度的规则网格(表1)。

　　表1叶尖位移

　　此外，通过改变参数s计算的结果如图6和7所示，其中也考虑了不同的边界条件(分别为CC和SS)。在这两种情况下，这些图形都表明，具有梁元素的混合网格对畸变的敏感性与常规网格(2x2Q4和2x2Q9)相当。2xQ4、1xQ9和2xQ4、1xB3网格的行为也非常相似;此外，可以注意到，在网格2xQ4 1xB2中使用B2梁单元并不会影响网格变形时结果的准确性，尽管它允许通过减少单元中节点的总数来区分单元边缘的节点数量(见图3b)。

　　图6

　　

　　CC四分之一板比例的变化

　　图7

　　

　　SS四分之一板的变化率

　　最后，考虑参数s的不同值，进行自由振动分析。图8和图9显示了与前16个模态相关的固有频率的结果:在此分析中，没有应用边界条件，因此只报告了非零模态。显然，之前的结论得到了证实:混合网格的解介于规则网格、2x2Q4和2x2Q9之间;混合网格2xQ4 1xQ9和2xQ4 1xB3除SS边界条件外，结果相同;解决方案2xQ4 1xB2比2x2Q4规则网格提供更低的固有频率，因此可以假设它更准确。

　　图8

　　

　　不同s值的四分之一板cc -模态分析

　　图9

　　

　　不同s值的四分之一板ss模态分析

　　现在考虑用一个完整的方形板来进行静力和模态分析。板的面内尺寸为，厚度为10 m。所应用的边界条件是固定的边缘(如图10a所示的上边缘)，载荷表示为右下角的集中力等于，同样在图10a中表示为交叉圆。所考虑的材料是各向同性的，具有以下特性:，和。请注意，我们使用商业软件Nastran分析了一个参考问题，如图10b所示。开始网格是沿厚度B3近似的平面上的Q4个元素(dof=225)，以及沿厚度B3近似的平面上的Q9个元素(dof=729)。

　　图10

　　

　　方形板-规则网

　　从这些规则网格开始，元素逐渐扭曲，以展示所呈现的自适应元素的能力。考虑到高斯分布的均值等于原坐标，标准差逐渐增大，一些元素的节点坐标是随机变化的。设为节点的原始坐标、标准差和均值，得到新的坐标如下:，其中Matlab函数randn返回从标准正态分布中抽取的伪随机值，标准差在等于的范围内变化。

　　此外，在其他示例中，有意移动节点以获得特定形状(例如，五边形，这是虚拟元素的基准示例[3,5,30])。最后一个例子是考虑一个完全的不规则变形网格。对于所有这些情况，利用板和梁元素的组合。

　　前面的网格报告如图11所示，其中:

　　例A:图11a中网格9xQ9 1xQ4 6xB2 (DOFs=387) -网格的标准偏差等于;

　　案例B:图11b中的Mesh 3xQ9 3xQ4 10xB2 (DOFs=330) -创建一些五边形元素获得网格;

　　案例C:图11c中Mesh 3xQ9 3xQ4 10xB2 v2 (DOFs=330)——得到先前网格的升级版本，突出五边形元素的存在，使得相同DOFs数的网格更加扭曲;

　　图11d中的情形D: 9xQ9 8xB3 (DOFs=765) -一个完整的变形网格被获得，并在运动学近似(二次插值)中使用高阶精度建模;

　　案例E:在图11e中，网格10xQ4 7xB2 (DOFs=315) -一个与前一个类似的网格被考虑，使用较低阶精度的运动学近似(线性插值)，以研究是否有可能获得较少自由度的可比结果。

　　图11

　　

　　方形板变形网格

　　4.3.1 静力与模态分析

　　通过进行静力分析对不同网格进行了研究，并根据底部自由边缘的位移报告了结果。特别是在图12中，显示了上述不同网格的位移幅度。对每个网格进行分析，并与Q4和Q9规则网格相关的曲线族进行比较。特别地，我们考虑:情况A，其中高阶精度元素(Q9)和低阶精度元素(Q4, B2)之间几乎处于平衡状态;情形B，其主要元素与Q4规则网格相似;情况C与前一种情况类似，不同之处在于现在五边形沉浸在更扭曲的网格中;情况D和E，由于单元的严重畸变，第一种情况与Q9畸变网格相关的参考解相匹配，第二种情况与Q4畸变网格相关的参考解相匹配。综上所述，变形网格的结果明显介于Q4和Q9规则网格的两条参考曲线之间，随着相关自由度的增加，精度也在增加。

　　图12

　　

　　方板静力分析

　　与静力分析类似，模态研究也将频率与Q4和Q9规则网格的解进行比较，如图13所示。注意，扭曲网格的频率与规则网格的频率完美地结合在一起，特别是在D非常接近完整的Q9解的情况下。

　　图13

　　

　　方板-模态分析

　　本文提出了一种基于CUF的变形网格自适应单元评估方法。在CUF框架中，这是一种先进的工具，能够以统一的方式制定结构理论，NDK方法与3D建模方法一起被利用。通过允许在同一元素上选择不同的展开，可以对某些节点进行运动学细化。这个方面代表了一个关键点，可以减少元素内的节点数量，从而通过保持相同的精度来减少计算成本。

　　在使用有限元法处理网格划分技术时，主要问题是由于计算成本的原因，单元在二维中必须是三角形或四边形，而在三维中必须是四面体或六面体，因此难以用单元逼近物理域。所提出的自适应元件可以克服这些问题，提供了一种方便的解决方案。

　　这些元素的能力已经在其他作品中进行了研究，参见[28,31]。这项工作的目的是证明不同的多边形和多面体可以使用自适应元素构建，具有相对较低的自由度，并评估它们在严重网格畸变出现时的性能。进行了静力和模态分析，给出了位移和固有频率方面的结果。结合具有不同运动近似的单元，给出了不同变形网格的结果。标准试验和算例分析表明，从计算的角度来看，可以实现方便的变形网格，变形效果与经典单元相当。

　　根据最近作者的工作[32]，进一步的研究应该关注自适应单元在壳结构分析中的应用，依靠这些单元在曲线坐标中的表述。

　　ccDownload: /内容/ pdf / 10.1007 / s42496 - 023 - 00165 - 6. - pdf