变刚度复合板中精确的3D应力恢复：一种高阶等几何方法

《COMPOSITE STRUCTURES》：Accurate 3D stress recovery in variable stiffness composite plates: A higher-order isogeometric approach

　　萨贾德·詹格拉维（Sajad Jangravi）|亚历山德罗·雷阿里（Alessandro Reali）|佩德拉姆·卡涅赫·马什杰迪（Pedram Khaneh Masjedi）
帕维亚大学土木工程与建筑系，Adolfo Ferrata街3号，帕维亚，27100，意大利

**摘要**
具有曲线形纤维的可变刚度复合材料由于其在空间上变化的刚度特性、强烈的各向异性效应以及复杂的三维应力梯度，给计算带来了相当大的挑战。因此，需要适当的理论来准确描述这种非均匀性引起的局部非线性横向剪切分布，尤其是在剪切变形效应显著的厚层板中。为了解决这一问题，本研究介绍了一种用于可变刚度复合板的三维应力恢复方法，该方法结合了三阶剪切变形理论（TSDT）和等几何分析（IGA）。通过利用非均匀有理B样条（NURBS）基函数的高阶连续性，所提出的方法能够准确地恢复层间的法向和剪切应力。与各种可变刚度层板配置下的3D有限元解进行综合比较，显示出极好的一致性，突显了该方法在预测平面内和横向应力分量方面的稳健性。结果表明，基于TSDT的等几何方法尤其在具有非均匀刚度分布的厚板中具有更高的精度和收敛速度。因此，这一框架为先进可变刚度复合结构的设计和分析提供了一个强大且计算效率高的工具。

**1. 引言**
由于层压复合材料相对于传统各向同性材料具有众多优势，它们已在许多工程领域得到了广泛应用。它们出色的强度重量比，以及通过层叠和纤维方向调整刚度特性的能力，使其在高性能应用中非常有价值[1]、[2]、[3]、[4]。随着复合结构制造技术的进步，可变刚度复合材料（VSCs）作为一类具有前沿性的定制层压材料出现，能够通过策略性地调节纤维方向、层厚和材料组成来实现空间变化的机械性能。制造技术，特别是自动纤维放置（AFP）和拖拽导向，使得在VSCs中实现曲线形纤维路径成为可能，从而提高了承载能力、损伤耐受性以及抗屈曲和振动性能，优于直纤维设计[5]、[6]、[7]、[8]、[9]。然而，现有方法往往无法捕捉具有较大横向剪切变形的板中的垂直于平面方向的应力。此外，VSCs中的空间刚度梯度导致高度非均匀的行为和较大的应力集中，这突显了需要更高效和准确的应力分析方法，这对于这些结构的可靠设计至关重要。
传统的板理论，如经典层压板理论（CLPT）和一阶剪切变形理论（FSDT），在预测厚板或高度非均匀层的局部应力场时可能存在局限性。CLPT适用于薄且刚度恒定的板；然而，它忽略横向剪切变形，导致较厚层的应力预测不够准确。FSDT通过包含一阶剪切效应来解决这个问题；但是，它仍然依赖于剪切修正因子来近似局部剪切应变分布[10]。虽然这种方法对标准层板的整体响应有效，但它不足以捕捉可变刚度复合材料中的复杂非线性剪切应变分布，这类材料的特征是空间刚度梯度和高非均匀性，而恒定剪切假设无法表示这些特性。VSCs与传统层压材料的一个根本分析区别在于，其控制平衡方程的特征系数是可变的，这是由于纤维方向的连续空间变化所致。与刚度恒定的板不同，在后者中材料属性在微分过程中被视为标量，而VSCs的非均匀性需要显式纳入刚度的空间导数以满足局部平衡。CLPT和FSDT在解决厚板和空间刚度变化方面的基本局限性表明，需要更先进的计算方法，如高阶剪切变形理论（HSDT），后者无需剪切修正因子即可准确模拟非线性横向剪切变形[11]。Akhavan等人[12]研究了具有曲线形纤维的VSC层板的大变形和应力，使用TSDT展示了通过优化纤维方向能够改善应力重新分布并减少高应力集中。Yazdani和Ribeiro[13]引入了一种p-version层析有限元方法，用于具有曲线形纤维的非对称复合板的静态分析。Amirpour等人[14]提出了用于功能梯度厚板的弹性变形的分析解，应用了HSDT。Gupta和Pradyumna[15]开发了一种基于HSDT的几何非线性有限元方法，用于VSC壳板的弯曲分析，并展示了曲线形纤维方向和层压配置对非线性变形和应力分布的影响。Sharma等人[16]将分层分析与剪切变形理论结合，评估了智能VSC板的静态性能，展示了曲线形纤维层压板在减少振动引起的变形和提高损伤耐受性方面的优势。Moreira等人[17]通过采用高阶剪切变形理论改进了层析有限元模型，分析了具有曲线形纤维的VSC层的静态响应，并证明了其方法在预测复杂应力分布方面的准确性。此外，Groh等人在[18]、[19]中提出了一种计算效率高的二维等效单层模型，能够提供异质层压板中的内在平衡三维应力预测。该模型提高了VSCs分析的准确性和计算效率。最近，Iannotta等人[20]和Zamani等人[21]扩展了Carrera的统一公式，以有效模拟可变刚度复合材料，证明了高阶运动学对这些复杂材料系统的重要性。

由于曲线形纤维路径会导致特定层的刚度变化并产生复杂的三维应力场，从而控制承载性能和失效起始，因此在VSC板中进行准确的应力恢复至关重要。通过重构完整的应力状态并计算位移，可以识别裂纹、分层或基体失效频繁发生的局部峰值区域。Sousa等人[22]、Falcó等人[23]、Thomas等人[24]和Zoghipour等人[25]使用LaRC03准则解决了这一问题，这些准则是一组基于物理的方程，假设平面应力并使用六个准则来预测平面内的基体和纤维失效模式。然而，这些方法依赖于理想化的纤维配置，假设材料属性均匀，或者需要大量的计算资源，这些资源往往会平滑掉拖拽下降区域或无纤维区域等特征。在应力恢复方法中，基于平衡的恢复方法将平衡形式直接纳入重构步骤。它需要平衡每层内部以及整个板厚度的内部力[26]。这种集成的力平衡改善了关键区域的应力预测，提供了更准确、高分辨率的厚层、非均匀VSC层压板在载荷下的行为表示。

等几何分析（IGA）因其能够准确表示这些先进材料中复杂的几何形状和空间刚度变化而被认为是分析VSCs的有效计算方法。Le-Manh等人[27]使用IGA对可变厚度复合板进行了非线性弯曲分析，并展示了其在表示复杂变形行为方面的准确性。Hao等人[28]开发了一个结合了准确建模、等几何分析和优化的综合模型，用于VSC面板。此外，Hao等人[29]开发了一个自适应梯度增强克里金模型，使用IGA对VSC面板进行全局优化，并指示了设计效率和计算成本的显著改进。NURBS的平滑性和高连续性使它们适合准确表示与曲线形纤维结构相关的高刚度梯度，因为NURBS基函数可以有效地捕获这些几何形状，而无需进行大量的网格细化。最近的贡献进一步巩固了IGA在这一领域的作用。例如，Ding等人[30]最近提出了一种适用于具有任意切口层压板的自适应IGA方案，而Chen等人[31]成功利用IGA和更高阶理论捕捉了可变刚度夹层板中的起皱现象。

基于IGA的平衡应力恢复方法能够利用NURBS基函数的平滑性和高连续性，以强形式应用平衡方程，从而提高层压复合板分析的准确性。Dufour等人[32]开发了一种计算效率高的基于平衡的方法，使用IGA框架中的3D实心元素来准确恢复垂直于平面方向的应力，证明了其对薄多层复合层的重要性。IGA使用在元素边界处具有C1或更高连续性的平滑NURBS基函数，从而提高了位移导数计算的准确性。利用这种连续性，最近的基于IGA的平衡应力恢复方法[33]、[34]、[35]、[36]、[37]集成了均质层特性和单一贯穿厚度的NURBS元素，以较低的计算成本实现准确的平面内和垂直于平面方向的应力预测。此外，Jangravi等人[38]将该方法扩展到各向异性复合板，解决了各向异性耦合效应，并准确重构了交叉层和角度层层压板中的三维垂直于平面方向的应力。然而，这些公式仅针对刚度恒定的薄层板进行了验证，未考虑显著的横向剪切应变。尽管最近的基于TSDT的等几何方法[39]成功解决了厚且刚度恒定的各向异性板的这一限制，但它没有考虑刚度变化。将TSDT结合到基于IGA的框架中，如本文所提出的，解决了这些综合挑战，有效地捕捉了显著的剪切变形、复杂的曲线形纤维路径和VSC层板中的完整三维应力状态。

本文首次提出了一种基于平衡的应力恢复方法，该方法在IGA框架内结合了TSDT，以满足大剪切效应下VSC层板的准确且高效的应力恢复需求。TSDT通过准确模拟非线性横向剪切分布，无需剪切修正因子，提高了VSCs中应力重建的准确性。IGA中NURBS基函数的连续性使得能够有效使用强形式平衡方程，从而实现准确且可靠的三维应力重构，同时考虑了纤维路径的变化。需要注意的是，将这一框架扩展到可变刚度复合材料涉及到与标准直纤维层压板相比的重大理论挑战。与刚度恒定的板不同，VSCs需要平衡方程考虑刚度矩阵的空间导数，这在应力恢复过程中引入了额外项。捕捉这些由梯度驱动的应力场需要点态本构更新，这在传统层压板公式中是缺失的。

以下是本文内容的总结：第2节回顾了Reddy的TSDT，并定义了层压复合板的运动学和本构关系。第3节介绍了双变量NURBS基函数并推导了等几何Galerkin离散化。第4节提出了一种基于平衡的后处理方法，该方法直接从平面TSDT解重构完整的三维应力场，特别是垂直于平面方向的分量。第5节使用涉及可变刚度层板的广泛基准问题验证了所提出的TSDT-IGA框架。与完整的3D有限元解的比较证明了该方法的准确性和收敛性。最后，第6节提供了结论性的评论。

**2. 可变厚度层压板的三阶剪切变形理论**
根据Reddy的三阶剪切变形理论[40]、[41]，层压复合板的位移场可以描述如下：
(1)
u(x,y,z) = u0 + zβx + cz
(2)
v(x,y,z) = v0 + zβy + cz
(3)
w(x,y,z) = w0

其中u0、v0和w0分别表示x、y和z方向的位移分量，c = ?4/3t2，t是层压板的总厚度。

**2.1. 弹性和运动学**
平面内应变向量?=[?xx,?yy,γxy]T可以表示为膜应变?0、弯曲曲率κ1和高阶曲率κ2的和：
(4)
? = ?0 + zκ1 + z3κ2
其中：
(4a) ?0 = ?u0/?x + ?v0/?y + ?v0/?x
(4b) κ1 = ?βx/?x + ?βy/?y + ?βx/?y
(4c) κ2 = ?βx/?x + ?2w0/?x2 + ?2w0/?y2 + ?2w0/?y2 + ?βy/?x + 2?2w0/?x?y

该公式涵盖了平面内和横向变形，z3依赖项包括高阶效应。横向剪切应变向量γ=[γxz,γyz]T定义为：
(5)
γ = ?s + z2κs
其中?s = βx + ?w0/?xβy + ?w0/?y，κs = 3cβx + ?w0/?xβy + ?w0/?y

第k层正交各向异性层的应力-应变关系由刚度矩阵Q?控制，将应力与应变关联起来：
(6)
σxx σyy σxy σxz σyz
k = Q?11Q?12Q?1600Q?21Q?22Q?2600Q?61Q?62Q?6600
000Q?55Q?54000Q?45Q?44
k ?xx ?yy γxy γxz γyz

变换后的刚度系数Q?ij取决于每层中的纤维方向θ，计算如下：
(7a) Q?11 = Q11cos?θ + Q22sin?θ + 2Q12 + 2Q66sin2θcos2θ
(7b) Q?12 = Q11 + Q22 ? 4Q66sin2θcos2θ + Q12sin?θ + cos?θ
(7c) Q?16 = Q11 ? Q12 ? 2Q66sinθcos3θ ? Q22 ? Q12 ? 2Q66cosθsin3θ
(7d) Q?22 = Q11sin?θ + 2Q12 + 2Q66sin2θcos2θ + Q22cos?θ
(7e) Q?26 = Q11 ? Q12 ? 2Q66cos3θsinθ ? Q22 ? Q12 ? 2Q66sin3θcosθ
(7f) Q?66 = Q11 + Q22 ? 2Q12 ? 2Q66sin2θcos2θ + Q66sin?θ + cos?θ
(7g) Q?44 = Q44cos2θ + Q55sin2θ
(7h) Q?45 = Q55 ? Q44sinθcosθ
(7i) Q?55 = Q55cos2θ + Q44sinθ

其中Q11 = E11 ? ν12ν21，Q12 = ν12E21 ? ν12ν21，Q22 = E21 ? ν12ν21，Q66 = G12，Q55 = G13，Q44 = G23。

设θ表示纤维方向角；E1和E2分别表示平行和垂直于纤维方向的杨氏模量；ν12和ν21是相应的泊松比；G12、G23和G13表示剪切模量。

**2.2. 变角拖拽（VAT）结构中的纤维方向变化**
由于制造过程中的实际限制，纤维通常会相互平行排列，因此纤维方向通常只沿着其中一个坐标变化。一种常见的技术是沿该路径使用线性纤维方向变化。在VAT结构中，纤维取向角度θ沿着曲线路径变化，导致在不同位置的层压板具有不同的刚度特性。Gurdal和Olmedo [42] 提供了以下方程来表示纤维角度的变化：(8)θ(x′)=Φ+T0+(T1?T0)d|x′|，其中Φ是关于参考方向的旋转角度，T0是参考点的纤维取向，T1是路径上某点的纤维取向。参数d表示纤维取向变化的特征距离，通常如图1所示等于a/2。坐标x′定义为：x′=xcosΦ+ysinΦ。

下载：下载高分辨率图像（191KB）
下载：下载全尺寸图像
图1. 具有不同纤维取向T0和T1的VAT模型的示意图。该方程允许纤维路径围绕某个参考方向旋转角度Φ，其中x′沿纤维的曲线路径测量。随着x′的增加，纤维取向根据方程中给出的线性变化逐渐从T0变为T1。

2.3. 变刚度复合材料的弱形式
在VSC中，由于纤维取向、材料分布或其他变量的变化，刚度矩阵在域内是不同的。因此，必须通过在每个伽辽金积分点计算这些矩阵来考虑这些变化。在域内的给定点，刚度矩阵A、B、D、E、F、H由局部材料属性和纤维取向决定。因此，弱形式被修改为包括这些空间依赖的刚度矩阵。内部虚功表达式现在写为：(9)????δWint=∫Ω(δ?mTD?m(x,y)?m+δ?sTD?s(x,y)?s)dΩ，其中(10)???m=?0κ1κ2T，?s=?sκsT。伽辽金点的坐标用(x,y)表示，刚度矩阵D?m(x,y)和D?s(x,y)是在每个伽辽金点评估的位置的函数，定义如下：D?m(x,y)=A(x,y)B(x,y)E(x,y)B(x,y)D(x,y)F(x,y)E(x,y)F(x,y)H(x,y)，(11)D?s(x,y)=As(x,y)Bs(x,y)Bs(x,y)Ds(x,y)。这些刚度矩阵取决于每个伽辽金点的材料属性和纤维取向，通常从局部纤维角度和修改后的材料刚度矩阵Q?获得。对于变刚度复合材料，刚度矩阵可以写为：(12)(Aij,Bij,Dij,Eij,Fij,Hij)=∫?h/2h/2(1,z,z2,z3,z4,z6)Q?ij(x,y)dz(i,j=1,2,6)，(13)(Aijs,Bijs,Dijs)=∫?h/2h/21,z2,z4Q?ijx,ydz(i,j=4,5)，其中Q?ij(x,y)表示(x,y)处的修改后材料刚度系数。刚度矩阵的空间依赖性意味着控制方程的特征是变系数。因此，这里需要特别注意的是在每个伽辽金积分点逐点更新这些刚度矩阵，以准确捕捉VSC固有的空间刚度梯度。横向载荷P和牵引载荷Tx和Ty的外部虚功表示为：(14)δWext=∫ΩTxδu0dΩ+∫ΩTyδv0dΩ+∫ΩPδw0dΩ。

3. 双变量结构的NURBS基函数
3.1. 双变量基函数
双变量NURBS基函数可以使用两个单变量节点向量Ξd来定义，这些向量包括参数平面每个方向d的参数域。节点向量写为：(15)Ξd={ξ1d,…,ξmd+pd+1d},d=1,2，其中pd是多项式度数，md表示相应方向d中的基函数数量。单变量B样条基函数可以使用Cox–de Boor递归算法来构建。对于多项式度数pd=0，基函数表示为：(16)Nid,0d(ξd)=1,ξidd≤ξd<ξid+1d,0，否则。对于更高的多项式度数pd>0，基函数递归计算为：(17)Nid,pdd(ξd)=ξd?ξiddξid+pdd?ξiddNid,pd?1d(ξd)+ξid+pd+1d?ξid+1dNid+1,pd?1d(ξd)。在构建单变量基函数后，通过计算它们的张量积获得双变量基函数：(18)Bi,p(ξ)=∏d=12Nid,pdd(ξd)，其中i={i1,i2}，p={p1,p2}，ξ={ξ1,ξ2}。这使得可以使用控制点Pi来积分双变量基函数，从而表示复杂的双变量NURBS几何形状。所得几何形状定义为：(19)S(ξ)=∑iBi,p(ξ)Pi。控制点Pi指定几何形状，而双变量基函数提供了平滑的参数映射。这利用NURBS灵活且连续地表示二维几何形状。
3.2. 离散方程
每个元素内的位移场使用NURBS基函数来计算，包括平面内位移、横向挠度和旋转。这种近似是通过将位移表示为控制点位移的线性组合来获得的，每个位移乘以相应的形状函数：(20)Uh=∑A=1CPRAqA。位移向量定义为Uh=[u0,v0,w0,βx,βy]T，其中u0、v0和w0表示平面内位移和横向挠度，βx和βy表示旋转。控制点A的位移向量表示为qA=[uA,vA,wA,βxA,βyA]T。与控制点A相关的有理基函数表示为RA，CP是NURBS基函数的总数。通过将方程(20)中的位移场代入应变-位移关系，我们得到平面内应变、曲率和剪切应变的公式：(21)???m?s=∑A=1CPBAqA。应变-位移矩阵BA由几个子矩阵组成，表示不同的变形模式：(22)BA=BAmBAb1BAb2BAs0BAs1。这些子矩阵的定义如下：(23)BAm=RA,x00000RA,y000RA,yRA,x0000BAb1=000RA,x00000RA,y BAb2=cRA,xx00000RA,yy000RA,xyRA,xy000，(24)BAs0=00RA,xRA0000RARA,y (25)BAs1=3c00RA,xRA0000RARA,y
这里，BAm与膜应变相关，BAb1和BAb2对应于弯曲应变，而BAs0和BAs1考虑剪切应变。
3.3. 刚度矩阵
通过结合内部和外部虚功贡献，并使用等几何Galerkin方法来制定能量泛函的变化。将方程(21)代入内部虚功表达式，得到全局刚度矩阵K：(26)K=∫ΩBmBb1Bb2TD?m(x,y)BmBb1Bb2dΩ+∫ΩBs0Bs1TD?s(x,y)Bs0Bs1dΩ。全局刚度矩阵K包括来自所有系统组件的膜、弯曲和剪切刚度的组合贡献。

4. 应力恢复过程
为了准确恢复平面外应力，需要一个基于平衡的后处理步骤。这种方法改进了等几何Galerkin方法的二维框架，以实现三维应力场的恢复。三维平衡方程可以表示为：(27)??σ+b=0，其中??σ是应力张量的散度，b是体力向量。平衡方程可以按每个应力分量展开：(28a)?σxx?x+?σxy?y+?σxz?z=?bx，(28b)?σxy?x+?σyy?y+?σyz?z=?by，(28c)?σxz?x+?σyz?y+?σzz?z=?bz。

4.1. 平面外剪切应力的恢复
平面外剪切应力σxz(z)和σyz(z)定义为z的函数，z是穿过厚度的坐标。通过在厚度维度上积分平衡方程来进行评估：(29a)σxz(z)=?∫zζ(?σxx(ξ)?x+?σxy(ξ)?y+bx(ξ))dξ+σxz(ζ)，(29b)σyz(z)=?∫zζ(?σxy(ξ)?x+?σyy(ξ)?y+by(ξ))dξ+σyz(ζ)。这里，ζ表示z的参考值，积分常数使用复合梯形规则根据板材上下表面的边界条件计算。为了准确评估这些积分中的应力导数，必须考虑 constitutive matrix Q?(x,y) 的空间变化。与恒定刚度层压板不同，这里的平衡方程涉及变系数。这需要在应力恢复过程中明确应用乘积规则和链式法则，以严格考虑空间刚度梯度。具体来说，评估应力分量的空间导数需要计算：(30)???σ?x=?Q?(x,y)?x?+Q?(x,y)???x。此外，由于刚度梯度由连续的曲线纤维路径θ(x,y)控制，明确计算这个导数需要链式法则：(31)?Q??x=?Q??θ?θ?x。捕捉这些由梯度驱动的应力场需要在每个伽辽金积分点进行逐点constitutive更新，这是直线纤维公式所没有的数学复杂性。为了计算σzz(z)，通过对之前导出的剪切应力分量进行积分来执行厚度方向的积分：(32)σzz(z)=∫zζ?2σxx(ξ)?x2+?2σyy(ξ)?y2+2?2σxy(ξ)?x?y+?bx(ξ)?x+?by(ξ)?ydξ。

5. 数值结果
5.1. 与统一公式的验证
在这里，提出的高阶IGA公式通过与Iannotta等人的[20]的结果进行比较来验证，其中包括3D FEM参考解和各种用于完全夹紧的正方形VAT板的板理论，该板的尺寸为a=b=1m，厚度为t，如图2所示，受到均匀压力P=10kPa的作用。所有结果按以下方式归一化：(33)w?=wa,σ?xx=σxxP。我们考虑两个不同的示例，即案例1和案例2。案例1和案例2的正交各向异性属性在表1中呈现。案例1使用沿y轴线性变化的单层纤维路径，Φ=90°，T0=0°，T1=90°，因此θ(y)从y=0处的0°变化到y=b处的90°。案例2由双层层压板组成，底部层的Φ=0°，T0=90°，T1=45°，顶部层的Φ=0°，T0=0°，T1=45°。我们将我们的结果与[20]中定义的多种板理论进行比较：等效单层(ESL)位移模型ED2和ED6，具有二阶和六阶的厚度方向多项式展开；逐层位移模型2LD2、3LD4、LD2和LD4，它们在每层独立应用二阶和四阶展开；逐层(LW)混合变分模型2LM2、3LM4、LM2和LM4，结合每层的二阶和四阶位移和应力近似；Chebyshev–Lobatto增强混合变体2LM2*和LM4*，它们将Chebyshev–Lobatto基函数纳入混合公式；以及FSDT，它使用一阶剪切变形假设，以及经典的层压板理论CLT，后者忽略了横向剪切。如表2和表3所示，当前公式在所有厚度比下都优于每个二阶ESL模型(ED2)，并且在a/t=100和10时优于六阶ESL模型(ED6)，仅在最厚的情况(a/t=5)下略低于ED6。它始终比所有逐层位移方案(2LD2、3LD4、LD2、LD4)更准确，并且在a/t=5时与它们相当。相对于经典理论，当前方法在每个比率下都优于CLT，并且在最薄的情况下也优于FSDT。

下载：下载高分辨率图像（223KB）
下载：下载全尺寸图像
图2. 受均匀压载作用的夹紧变刚度复合层压板，具有变化的纤维取向。对于横向位移，当前公式的准确性与混合变分家族(LM2、2LM2、3LM4)在a/t=100和10时相当或更高，并且在a/t=5时与Chebyshev–Lobatto增强变体(LM4*、2LM2*)处于相同的误差范围内。对于平面内正应力，当前公式在a/t=100时是最准确的解决方案。在a/t=10和5时，仅被最高阶LW和混合变分模型(LD2、LD4、LM4、LM4*)以及ED6超过，但仍优于LM2、LM2*、ED2、FSDT和CLT。这些趋势证实了所提出的高阶IGA公式在整个厚度范围内的稳健性和高预测质量。关于数值收敛性，当前结果是使用100 × 100元素的网格（p=5）获得位移的结果，以及110 × 110元素的网格（p=6）获得σxx的结果，而参考3D FEM解决方案需要每侧140个元素和厚度方向18个元素。此外，为了完整性，图3展示了厚板、中等厚度板和薄板的长度方向的归一化横向位移剖面（案例1）和厚度方向的平面内正应力分布（案例2）。

表1. 案例1和案例2的材料属性。
案例 E1 E2=E3 G12 G13 G23
ν12 ν13 ν23
空单元 (MPa) (MPa) (MPa) (MPa) (MPa)
空单元 空单元
案例1 50.00 10.00 5.00 5.00 0.25 0.25 0.25
案例2 137.90 8.96 7.10 7.10 6.21 0.30 0.30 0.49

表2. 在(a2,b2,t2)处归一化的横向位移w?×10?6（案例1）。
方法 a/t=100 a/t=10 a/t=5
3DFEM 65 44 8.76 81.87 73
LM4 65 27 8.76 01.87 52
LM2 65 30 8.78 11.89 03
LD4 65 26 8.75 91.87 62
LD2 65 26 8.75 71.87
ED6 65 22 8.43 51.74
FSDT 65 42 8.48 11.71
CLT 65 21 6.52 10.81
当前 65 53 8.76 91.80 4

表3. 在(a4,b2,t2)处归一化的平面内正应力σ?xx×10?1（案例2）。
方法 a/t=100 a/t=10 a/t=5
3DFEM 96 7.58 012.02 04.25
3LM4* 98 7.39 41 2.16 34.33
1LM4 98 7.50 012.16 34.33
1LM2* 98 3.39 71 11.58 73.81
LD4 100 8.79 21 2.29 14.33
LD2 98 98 48.29 34.29
ED6 99 49 94.29 12.56 5
CLT 99 2.82 09.92 42.48
当前 96 8.21 11.41 23.99 2

图3. 案例1的梁长度方向的归一化横向位移，以及案例2的平面内正应力的厚度方向归一化分布（当前TSDT-IGA）。

5.2. 变刚度复合层的数值验证
在本节中，我们分析了一个正方形VSC层压板，其平面尺寸为a=b=1m，总厚度为t=100mm。板沿边缘夹紧并受到均匀横向压力P=100kPa的作用。层压板的配置在表5中详细说明，根据方程(8)，纤维取向表示为[Φ+〈T0∣T1〉]s，此外，表6报告了在z=?t/5处恢复的垂直于平面应力分量的百分比误差，而表7则展示了z=?t/10处的相应结果。TSDT结果与使用149×149个平面元素和80个厚度元素的3D FEM结果之间的良好一致性，证实了所提出的方法能够准确捕捉到具有各种纤维角度引起的陡峭梯度的变刚度层压板中的平面内和垂直于平面的应力分布。表4列出了变刚度层压板的材料力学性能。

表5显示了变刚度复合材料的层压板厚度比和堆叠顺序：

| 层压板厚度比 | 堆叠顺序 |
|---------|---------|
| VAT A | [0.254][0±〈0|90〉]s |
| VAT B | [0.254][90±〈0|90〉]s |
| VAT C | [0.1258][0±〈0|45〉/0±〈90|45〉]s |
| VAT D | [0.1258][90±〈0|70〉/90±〈90|20〉]s |
| VAT E | [0.1258][90±〈0|45〉/0±〈45|0〕s |

图4展示了使用NURBS基函数（度数为p=5）时，VAT配置A-E的横向位移与元素数量（Nel）的收敛情况。

表6比较了z=?t/5处VAT配置（3D FEM vs. 现有IGA）的无量纲应力分量：

| 组配 | FEM | IGA | 差异（%） |
|---------|------------|------------|---------|
| VAT A | ?1.96 | ?1.92 | 1.76 | ?0.79 |
| VAT B | ?2.76 | ?2.65 | 3.82 | ?1.47 |
| VAT C | ?2.40 | ?2.30 | 4.03 | ?0.95 |
| VAT D | ?2.38 | ?2.40 | 0.78 | ?1.38 |
| VAT E | ?2.81 | ?2.78 | 1.29 | ?0.92 |

表7比较了z=t/10处VAT配置（3D FEM vs. 现有IGA）的无量纲应力分量：

| 组配 | FEM | IGA | 差异（%） |
|---------|------------|------------|---------|
| VAT A | ?2.49 | ?2.48 | 0.26 | ?0.94 |
| VAT B | ?3.12 | ?3.08 | 1.41 | ?1.87 |
| VAT C | ?2.47 | ?2.42 | 1.93 | ?1.73 |
| VAT D | ?2.44 | ?2.45 | 0.60 | ?1.76 |
| VAT E | ?3.19 | ?3.20 | 0.47 | ?1.22 |

图5和图6分别展示了VAT A和VAT B的厚度方向上的归一化应力分量。

图7至图10展示了VAT C至VAT E的厚度方向上的归一化应力分量。

6. 结论

本研究提出了一种三阶剪切变形理论，采用等几何Galerkin公式和基于平衡的恢复程序，用于三维分析变刚度复合板。所提出的方法利用NURBS基函数的高阶连续性，在存在空间刚度梯度的情况下强制实现强形式平衡。我们通过与其他现有数值结果及3D有限元解进行对比，验证了该框架的预测能力。结果表明，该方法从细板到厚板都能保持高精度，有效解决了仅靠单层理论难以捕捉的剪切应力分布问题。即使在纤维取向梯度和强各向异性耦合的情况下，该方法也能准确重建层间应力。这些发现表明，TSDT-IGA框架为先进变刚度复合板的应力分析提供了一种计算效率高且准确的工具。

CRediT作者贡献声明：
Sajad Jangravi：撰写 – 审稿与编辑、撰写 – 原稿、可视化、验证、软件开发、方法论、概念化。
Alessandro Reali：撰写 – 审稿与编辑、监督、方法论、概念化。
Pedram Khaneh Masjedi：撰写 – 审稿与编辑、监督、方法论、概念化。