CN108108578A

CN108108578A - 基于无网格法的fg-grc板屈曲载荷因子的数值算法

Info

Publication number: CN108108578A
Application number: CN201810089235.0A
Authority: CN
Inventors: 曾煌棚; 雷祖祥
Original assignee: Nanjing University of Science and Technology
Current assignee: Nanjing University of Science and Technology
Priority date: 2018-01-30
Filing date: 2018-01-30
Publication date: 2018-06-01

Abstract

本发明公开了一种基于无网格法的FG‑GRC板屈曲载荷因子的数值算法，包括以下步骤：构建模型，确定几何参数与材料参数；求解线性本构方程与总能量函数；构建无网格法的形函数并施加边界条件；计算刚度矩阵并通过特征值方程得出屈曲载荷因子。通过本发明可以准确的计算在外部热载荷及轴向压力作用下的功能梯度石墨烯增强复合材料层合板的屈曲载荷因子。在研究各类参数如：边界温度、石墨烯分布形式、层合板长宽比、宽高比和层数等，对功能梯度石墨烯增强复合材料层合板屈曲特性的影响上具有重要意义。

Description

基于无网格法的FG-GRC板屈曲载荷因子的数值算法

技术领域

本发明涉及功能梯度复合材料力学性能技术，具体涉及一种基于无网格法的FG-GRC板屈曲载荷因子的数值算法。

背景技术

目前在功能梯度复合材料的力学性能研究领域，应用最广的方法是有限元法。许多国内外的学者都使用有限元法来进行材料的各项力学性能分析。例如：《用有限元法研究复合材料层合板的热屈曲特性》一文便是采用有限元法来进行复合材料层合板的热屈曲分析。但有限元法在处理大变形问题时，存在着一些不足与缺陷。有限元法以离散的网格作为研究对象，使得有限元法的近似函数具有网格依赖性，网格畸变容易引起一些计算问题；有限元法的前处理不仅需要节点信息，还需要网格信息，不易进行复杂三维结构的分析；有限元法的计算结果是非光滑连续的，需要进行应力光滑等后处理；由于有限元法划分的网格类型的不同，使得计算结果有所偏差；使用有限元法进行计算时，需要处理网格与节点的信息，因此计算量大，计算时间长。而无网格法是在建立问题域的***代数方程时，不需要利用预定义的网格信息，或者只利用更容易生成的更灵活、更自由的网格进行域离散的方法。

发明内容

本发明的目的是提供一种基于无网格法的功能梯度石墨烯增强复合材料(FG-GRC)板屈曲载荷因子的数值算法，可以克服传统方法中对网格的依赖，用以研究各参数变量对功能梯度石墨烯增强复合材料层合板屈曲特性的影响。

实现本发明目的的技术解决方案为：一种基于无网格法的FG-GRC板屈曲载荷因子的数值算法，包括如下步骤：

步骤1、构建复合材料层合板的几何模型，定义功能梯度石墨烯增强复合材料层合板的几何参数和材料参数，并采用扩展Halpin-Tsai模型求解功能梯度石墨烯增强复合材料的整体材料参数；

步骤2、将复合材料层合板的几何模型进行离散，通过一阶剪切变形理论求得复合材料层合板离散节点的位移场和应变场，并求解线性本构方程：

u(x,y,z)＝u₀(x,y)+zΦ_x(x,y)

v(x,y,x)＝v₀(x,y)+zΦ_y(x,y)

w(x,y,z)＝w₀(x,y)

其中u,、v和w分别对应表示x,y和z方向上的位移分量，(u₀,v_o,w₀,Φ_x,Φ_y)是中平面z＝0的位移分量，Φ_x和Φ_y分别对应表示y轴正方向和x轴负方向对其横向法线的转动角度；

其中(ε_xx，ε_yy，γ_xy，γ_yz，γ_xz)为应变分量，ε₀、κ和γ₀均为中间变量；

其中(σ_xx，σ_yy，σ_xy，σ_yz，σ_xz)为应力分量，Q_ij为中间变量，α₁₁和α₂₂分别对应表示复合材料层合板在纵向和横向上的热膨胀系数，ΔT是外部温度的变化量；

步骤3、求解功能梯度石墨烯增强复合材料层合板的应变能U_ε和由平面应力引起的势能W_g，通过应变能与势能的差值得出总能量函数Π；

Π＝U_ε-W_g

其中ε为应变场，和为复合材料层合板所受载荷，S、γ₁和γ₂均为中间变量；

步骤4、基于再生核质点法，构建影响域内有效节点的无网格法的形函数；

步骤5、为了使无网格法的形函数满足Kronecker delta属性，采用完全变换法来施加本质边界条件，并通过转换矩阵来重构形函数；

步骤6、构建离散***方程，计算刚度矩阵，并采用特征值方程计算临界屈曲载荷λ，通过临界屈曲载荷得出屈曲载荷因子

(K+λK_g)u＝0

其中K和K_g均为刚度矩阵，u为位移场，a为复合材料层合板的长，h为复合材料层合板的高，E₂为复合材料层合板的杨氏模量。

本发明与现有的技术相比，有益效果在于：(1)针对传统的以网格为对象的数值计算方法的不足，提出了以节点为研究对象的数值模拟方法，有效避免了因网格畸变引起的问题以及带状稀疏系数矩阵引起的计算量过大的问题；(2)极大的提高了计算的速度与结果的准确性；(3)可以准确有效地进行不同参数条件下的功能梯度石墨烯增强复合材料层合板屈曲载荷因子的数值计算；(4)实践中发现，无网格法的形函数光滑性良好，数值模拟计算结果的精确性良好，对分析功能梯度复合材料的力学性能具有十分重要的意义。

附图说明

图1为本发明基于无网格法的FG-GRC板屈曲载荷因子的数值算法的流程图。

图2为本发明功能梯度石墨烯增强复合材料层合板的几何模型以及石墨烯的不同分布形式图，其中图a为石墨烯分布形式为UD的层合板几何图，图b为石墨烯分布形式为FG-O的层合板几何图，图c为石墨烯分布形式为FG-X的层合板几何图。

图3为外界温度不同的功能梯度石墨烯增强复合材料层合板的屈曲载荷因子曲线图。

图4为石墨烯分布形式不同的功能梯度石墨烯增强复合材料层合板的屈曲载荷因子曲线图。

图5为长宽比不同的功能梯度石墨烯增强复合材料层合板的屈曲载荷因子曲线图。

图6为宽高比不同的功能梯度石墨烯增强复合材料层合板的屈曲载荷因子曲线图。

图7为层数不同的功能梯度石墨烯增强复合材料层合板的屈曲载荷因子曲线图。

具体实施方式

为了使本发明的技术方案更加清楚明白，以下结合附图对本发明进行详细的的描述与解释：

结合图1，本发明所述的一种基于无网格法的FG-GRC板屈曲载荷因子的数值算法，包括如下步骤：

步骤1、构建复合材料层合板的几何模型，定义功能梯度石墨烯增强复合材料层合板的几何参数和材料参数，并采用扩展Halpin-Tsai模型求解功能梯度石墨烯增强复合材料的整体材料参数。

定义复合材料层合板的几何参数，即长为a、宽为b和高为h。确定复合材料层合板增强材料与基体材料的材料参数。复合材料层合板的增强材料为石墨烯，和均为石墨烯的杨氏模量，G^G为石墨烯的剪切模量，ν^G为石墨烯的泊松比。E^m为基体材料的杨氏模量，G^m为基体材料的剪切模量，ν^m是基体材料的泊松比。并采用扩展Halpin-Tsai模型来计算复合材料层合板的杨氏模量E₁₁、E₂₂和剪切模量G₁₂，计算过程如下所示：

其中，a_G为石墨烯片的长,b_G为石墨烯片的宽，h_G为石墨烯片的高；η_j(j＝1,2,3)为石墨烯的有效比例系数；

γ₁₁、γ₂₂和γ₁₂为中间变量；

其中，和分别对应表示石墨烯片的杨氏模量和剪切模量。E^m和G^m分别对应表示基体材料的杨氏模量和剪切模量。V_G和V_m＝1-V_G分别对应表示石墨烯与基体材料的体积分数。温度的变化对泊松比的影响很小，复合材料层合板的泊松比ν₁₂为：

ν₁₂＝V_Gν^G+V_mv^m

公式中的v^G和v^m分别对应表示石墨烯与基体材料的泊松比。根据Schapery模型，我们可以得出复合材料层合板在纵向和横向上的热膨胀系数α₁₁和α₂₂。

其中的均为石墨烯的热膨胀系数，α^m为基体材料的热膨胀系数。

步骤2、通过一阶剪切变形理论求得位移场和应变场，并求解线性本构方程。

通过一阶剪切变形理论可以得到位移场。

u(x,y,z)＝u₀(x,y)+zΦ_x(x,y)

v(x,y,x)＝v₀(x,y)+zΦ_y(x,y)

w(x,y,z)＝w₀(x,y)

其中u,v和w分别对应表示x,y和z方向上的位移分量，(u₀,v_o,w₀,Φ_x,Φ_y)是中平面z＝0的位移分量，Φ_x和Φ_y分别对应表示y轴正方向和x轴负方向对其横向法线的转动角度。我们得出应变分量如下。

其中(ε_xx，ε_yy，γ_xy，γ_yz，γ_xz)为应变分量，ε₀、κ和γ₀为中间变量；

线性本构关系如下所示：

其中(σ_xx，σ_yy，σ_xy，σ_yz，σ_xz)为应力分量，线性本构方程中的中间变量Q_ij的值为：

Q₆₆＝G₁₂,Q₄₄＝G₂₃,Q₅₅＝G₁₃

上述公式中的ΔT是外部温度的变化量，E₁₁,和E₂₂均为复合材料层合板的杨氏模量，G₁₂、G₁₃和G₂₃均为复合材料层合板的剪切模量，v₁₂和v₂₁均为复合材料层合板的泊松比。α₁₁和α₂₂均为复合材料层合板的热膨胀系数。

步骤3、求解功能梯度石墨烯增强复合材料层合板的应变能U_ε以及由平面应力引起的势能W_g，通过应变能与势能的差值可以得出总能量函数Π。

功能梯度石墨烯增强复合材料层合板的应变能为：

其中ε为应变场，S为中间变量；

拉伸刚度A、耦合刚度弯曲刚度D和横向剪切刚度A^s分别为：

其中K为横向剪切修正系数，功能梯度复合材料的横向剪切修正系数为：K＝5/(6-(v₁V₁+v₂V₂))，其中v₁和v₂分别对应表示增强材料与基体材料的泊松比，V₁和V₂分别对应表示增强材料与基体材料的体积分数。

层合板刚度的计算公式为：

表示的是层合板第k层变换后的刚度矩阵，如下所示：

变换矩阵[T]为：

其中，θ表示的是层合板的层交角度。

由平面应力引起的势能为：

其中和为复合材料层合板所受载荷大小，γ₁和γ₂均为中间变量，由应变能与势能的差值得出总能量函数为：

Π＝U_ε-W_g

步骤4、基于再生核质点法，构建影响域内有效节点的无网格法的形函数。

离散化的一般位移场为：

其中，ψ_I(x)、u_I和NP分别对应表示形函数、节点参数与离散节点数量。

基于再生核质点法，得到的形函数为：

ψ_I(x)＝C(x；x-x_I)Φ_a(x-x_I)

Φ_a(x-x_I)为满足再生条件的核函数，C(x；x-x_I)为满足再生条件的修正函数。

修正函数C(x；x-x_I)是由多项式基函数的线性组合表示而成。

C(x；x-x_I)＝H^T(x-x_I)b(x)

其中H^T(x-x_I)和b(x)均为中间变量：

b(x)＝[b₀(x,y),b₁(x,y),b₂(x,y),b₃(x,y),b₄(x,y),b₅(x,y)]^T

H^T(x-x_I)＝[1,x-x_I,y-y_I,(x-x_I)(y-y_I),(x-x_I)²,(y-y_I)²]

形函数Ψ_I(x)可以表示为：

Ψ_I(x)＝b^T(x)H(x-x_I)Φ_a(x-x_I)

b(x)为：

b(x)＝M^-1(x)H(0)

其中

H(0)＝[1,0,0,0,0,0,]^T

其中M(x)为中间变量，H^T(x-x_I)和H(x-x_I)互为转置矩阵，Φ_a(x-x_I)为二维核函数：

Φ_a(x-x_I)＝Φ_a(x)Φ_a(y)

其中Φ_a(x)和Φ_a(y)为中间变量：

选用了三次样条函数作为权函数：

其中z_I为中间变量：

其中d_max表示的是比例因子，d_I表示的是扩张系数。为了避免矩阵M的奇异性，在选择距离c_I的时候必须遍历足够多的节点，(x_I，y_I，z_I)为节点I的坐标。。

可以将形函数改写为：

Ψ_I(x)＝H^T(0)M^-1(x)H(x-x_I)Φ_a(x-x_I)

M(x)b(x)＝H(0)

对矩阵M(x)进行LU分解可以用来计算b(x)，可以得出

M_x(x)b(x)+M(x)b_x(x)＝H_x(0)

M(x)b_x(x)＝H_x(0)-M_x(x)b(x)

步骤5、为了使无网格法的形函数满足Kronecker delta属性，采用完全变换法来施加本质边界条件，并通过转换矩阵来重构形函数。

在离散位移场的基础上，广义位移u表示为：

其中是节点x_J处的位移，L_IJ＝Ψ_I(x_J)为形函数。

由上述两个公式可以推导出：

其中

其中δ_IJ为中间变量，为满足Kronecker delta属性的形函数。

步骤6、构建离散***方程并计算刚度矩阵，并采用特征值方程计算临界屈曲载荷λ，通过临界屈曲载荷得出屈曲载荷因子a为复合材料层合板的长，h为复合材料层合板的高，E₂为复合材料层合板的杨氏模量。

对总能量函数执行Ritz变换可得到离散***方程：

(K+λK_g)u＝0

其中λ为临界屈曲载荷，u为位移场，刚度矩阵K和K_g的计算公式为：

Κ＝Κ^b+Κ^m+Κ^s

其中，Κ^b，Κ^m，Κ^s、和均为中间变量，在上述刚度矩阵的计算公式中使用稳定节点积分和直接节点积分可以得到：

其中A_L为离散域面积；

其中n_x(x_L)和n_y(x_L)为节点法向量的值；和为轴向载荷；和均为中间变量。

实施例

为了验证该算法方案的可行性，采用MATLAB编程软件对该算法编写计算程序，代入具体的参数，求解功能梯度石墨烯增强复合材料层合板的屈曲载荷因子。

步骤1、构建模型，定义功能梯度石墨烯增强复合材料层合板的几何参数和材料参数，并采用扩展Halpin-Tsai模型求解功能梯度石墨烯增强复合材料的整体材料参数。

功能梯度石墨烯增强复合材料层合板的几何参数为长a＝0.2m、宽b＝0.2m和高h＝0.02m。石墨烯增强材料的分布形式为均匀分布和功能梯度分布，分别为：UD、FG-O和FG-X。复合材料层合板的几何参数和石墨烯的分布形式如图2所示。为了实现石墨烯的功能梯度分布，我们将复合材料板在厚度方向上均匀分割为十层，通过改变每一层的石墨烯增强材料的体积分数，我们就可以实现石墨烯在厚度方向上的功能梯度分布。对于UD分布形式，我们定义石墨烯在每一层中的体积分数为：[(0.07)/(0.07)/(0.07)/(0.07)/(0.07)/(0.07)/(0.07)/(0.07)/(0.07)/(0.07)]；对于FG-O分布形式，我们定义石墨烯在每一层中的体积分数为：[(0.03)/(0.05)/(0.07)/(0.09)/(0.11)/(0.11)/(0.09)/(0.07)/(0.05)/(0.03)]；对于FG-X分布形式，我们定义石墨烯在每一层中的体积分数为：[(0.11)/(0.09)/(0.07)/(0.05)/(0.03)/(0.03)/(0.05)/(0.07)/(0.09)/(0.11)]。复合材料层合板存在四种边界条件，分别为四边简支(SSSS)、四边固支(CCCC)、两边简支两边固支(SCSC)以及两边简支两边自由(SFSF)。选定三种载荷类型分别为：单轴压缩(γ₁＝-1，γ₂＝0，γ₃＝0)、双轴压缩(γ₁＝-1，γ₂＝-1，γ₃＝0)和三轴压缩(γ₁＝-1，γ₂＝-1，γ₃＝-1)。

功能梯度石墨烯增强复合材料是将石墨烯添加进甲基丙烯酸甲酯聚合物基体中复合形成的，为了对复合材料的力学性能进行分析，必须先计算出复合材料层合板的整体材料参数。石墨烯与基体的材料参数是已知的。聚合物基体的泊松比为v^m＝0.34，聚合物基体的杨氏模量为E^m＝(3.52-0.0034T)Gpa(其中T＝T₀+ΔT，T₀＝300K),聚合物基体的剪切模量为G^m＝E^m/(2*(1+v^m))GPa。石墨烯材料的几何参数为：长a_G＝14.76nm、宽b_G＝14.77nm和高h_G＝0.188nm，泊松比为ν^G＝0.177，在不同温度下石墨烯材料的杨氏模量和剪切模量如表1所示。

表1不同温度下石墨烯的材料参数

使用扩展Halpin-Tsai模型来实现功能梯度石墨烯增强复合材料有效材料参数的计算：

其中η_j(j＝1,2,3)为石墨烯的有效比例系数。η_j的值如表2所示。

表2不同温度和体积分数下，石墨烯的有效比例系数

位移场为：

u(x,y,z)＝u₀(x,y)+zΦ_x(x,y)

v(x,y,x)＝v₀(x,y)+zΦ_y(x,y)

w(x,y,z)＝w₀(x,y)

应变场为：

线性本构关系如下所示：

步骤3、求解功能梯度石墨烯增强复合材料层合板的应变能U_ε以及由平面应力引起的势能W_g，通过应变能与势能的差值得出总能量函数П。

П＝U_ε-W_g

步骤4、基于再生核质点法，构建影响域内有效节点的无网格法的形函数Ψ_I(x)。

Ψ_I(x)＝H^T(0)M^-1(x)H(x-x_I)Φ_a(x-x_I)

步骤5、采用完全变换法来施加本质边界条件，使无网格法的形函数满足Kronecker delta属性，并通过转换矩阵来重构形函数

步骤6、对总能量函数执行Ritz变换得到离散***方程：

(K+λK_g)u＝0

采用特征值方程计算临界屈曲载荷λ，通过临界屈曲载荷得出屈曲载荷因子

图3-7分别表示的是不同参数条件下功能梯度石墨烯增强复合材料层合板的屈曲载荷因子。结果表明，该发明能有效的进行功能梯度石墨烯增强复合材料层合板屈曲载荷因子的计算，并且通过改变参数变量可以用于分析不同的参数对功能梯度石墨烯增强复合材料层合板屈曲特性的影响，验证了该数值计算方法的可行性与准确性。

综上所述，相比于有限元法将离散的单元作为研究对象，无网格法的研究对象为离散节点。无网格法的近似函数没有网格依赖性，减少了因网格畸变而引起的困难；无网格法的前处理只要结点位置信息，不用网格信息，容易分析复杂的三维结构；无网格法计算的结果是光滑连续的，不必再进行应力光滑等后处理。无网格方法相对于有限元在处理大变形问题方面有很大的优势。因此，采用无网格法可以高效且准确地进行功能梯度复合材料的力学分析。

Claims

1.一种基于无网格法的FG-GRC板屈曲载荷因子的数值算法，其特征在于，包括如下步骤：

u(x,y,z)＝u₀(x,y)+zΦ_x(x,y)

v(x,y,x)＝v₀(x,y)+zΦ_y(x,y)

w(x,y,z)＝w₀(x,y)

<mrow> <mfenced open = "{" close = "}"> <mtable> <mtr> <mtd> <msub> <mi>&epsiv;</mi> <mrow> <mi>x</mi> <mi>x</mi> </mrow> </msub> </mtd> </mtr> <mtr> <mtd> <msub> <mi>&epsiv;</mi> <mrow> <mi>y</mi> <mi>y</mi> </mrow> </msub> </mtd> </mtr> <mtr> <mtd> <msub> <mi>&gamma;</mi> <mrow> <mi>x</mi> <mi>y</mi> </mrow> </msub> </mtd> </mtr> </mtable> </mfenced> <mo>=</mo> <msub> <mi>&epsiv;</mi> <mn>0</mn> </msub> <mo>+</mo> <mi>z</mi> <mi>&kappa;</mi> <mo>,</mo> <mfenced open = "{" close = "}"> <mtable> <mtr> <mtd> <msub> <mi>&gamma;</mi> <mrow> <mi>y</mi> <mi>z</mi> </mrow> </msub> </mtd> </mtr> <mtr> <mtd> <msub> <mi>&gamma;</mi> <mrow> <mi>x</mi> <mi>z</mi> </mrow> </msub> </mtd> </mtr> </mtable> </mfenced> <mo>=</mo> <msub> <mi>&gamma;</mi> <mn>0</mn> </msub> </mrow>

<mrow> <mfenced open = "{" close = "}"> <mtable> <mtr> <mtd> <msub> <mi>&sigma;</mi> <mrow> <mi>x</mi> <mi>x</mi> </mrow> </msub> </mtd> </mtr> <mtr> <mtd> <msub> <mi>&sigma;</mi> <mrow> <mi>y</mi> <mi>y</mi> </mrow> </msub> </mtd> </mtr> <mtr> <mtd> <msub> <mi>&sigma;</mi> <mrow> <mi>x</mi> <mi>y</mi> </mrow> </msub> </mtd> </mtr> <mtr> <mtd> <msub> <mi>&sigma;</mi> <mrow> <mi>y</mi> <mi>z</mi> </mrow> </msub> </mtd> </mtr> <mtr> <mtd> <msub> <mi>&sigma;</mi> <mrow> <mi>x</mi> <mi>z</mi> </mrow> </msub> </mtd> </mtr> </mtable> </mfenced> <mo>=</mo> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <msub> <mi>Q</mi> <mn>11</mn> </msub> </mtd> <mtd> <msub> <mi>Q</mi> <mn>12</mn> </msub> </mtd> <mtd> <mn>0</mn> </mtd> <mtd> <mn>0</mn> </mtd> <mtd> <mn>0</mn> </mtd> </mtr> <mtr> <mtd> <msub> <mi>Q</mi> <mn>12</mn> </msub> </mtd> <mtd> <msub> <mi>Q</mi> <mn>22</mn> </msub> </mtd> <mtd> <mn>0</mn> </mtd> <mtd> <mn>0</mn> </mtd> <mtd> <mn>0</mn> </mtd> </mtr> <mtr> <mtd> <mn>0</mn> </mtd> <mtd> <mn>0</mn> </mtd> <mtd> <msub> <mi>Q</mi> <mn>66</mn> </msub> </mtd> <mtd> <mn>0</mn> </mtd> <mtd> <mn>0</mn> </mtd> </mtr> <mtr> <mtd> <mn>0</mn> </mtd> <mtd> <mn>0</mn> </mtd> <mtd> <mn>0</mn> </mtd> <mtd> <msub> <mi>Q</mi> <mn>44</mn> </msub> </mtd> <mtd> <mn>0</mn> </mtd> </mtr> <mtr> <mtd> <mn>0</mn> </mtd> <mtd> <mn>0</mn> </mtd> <mtd> <mn>0</mn> </mtd> <mtd> <mn>0</mn> </mtd> <mtd> <msub> <mi>Q</mi> <mn>55</mn> </msub> </mtd> </mtr> </mtable> </mfenced> <mrow> <mo>(</mo> <mrow> <mfenced open = "{" close = "}"> <mtable> <mtr> <mtd> <msub> <mi>&sigma;</mi> <mrow> <mi>x</mi> <mi>x</mi> </mrow> </msub> </mtd> </mtr> <mtr> <mtd> <msub> <mi>&sigma;</mi> <mrow> <mi>y</mi> <mi>y</mi> </mrow> </msub> </mtd> </mtr> <mtr> <mtd> <msub> <mi>&sigma;</mi> <mrow> <mi>x</mi> <mi>y</mi> </mrow> </msub> </mtd> </mtr> <mtr> <mtd> <msub> <mi>&sigma;</mi> <mrow> <mi>y</mi> <mi>z</mi> </mrow> </msub> </mtd> </mtr> <mtr> <mtd> <msub> <mi>&sigma;</mi> <mrow> <mi>x</mi> <mi>z</mi> </mrow> </msub> </mtd> </mtr> </mtable> </mfenced> <mo>-</mo> <mfenced open = "{" close = "}"> <mtable> <mtr> <mtd> <msub> <mi>&alpha;</mi> <mn>11</mn> </msub> </mtd> </mtr> <mtr> <mtd> <msub> <mi>&alpha;</mi> <mn>22</mn> </msub> </mtd> </mtr> <mtr> <mtd> <mn>0</mn> </mtd> </mtr> <mtr> <mtd> <mn>0</mn> </mtd> </mtr> <mtr> <mtd> <mn>0</mn> </mtd> </mtr> </mtable> </mfenced> <mi>&Delta;</mi> <mi>T</mi> </mrow> <mo>)</mo> </mrow> </mrow>

步骤3、求解功能梯度石墨烯增强复合材料层合板的应变能U_ε和由平面应力引起的势能W_g，通过应变能与势能的差值得出总能量函数П；

<mrow> <msub> <mi>U</mi> <mi>&epsiv;</mi> </msub> <mo>=</mo> <mfrac> <mn>1</mn> <mn>2</mn> </mfrac> <msub> <mo>&Integral;</mo> <mi>&Omega;</mi> </msub> <msup> <mi>&epsiv;</mi> <mi>T</mi> </msup> <mi>S</mi> <mi>&epsiv;</mi> <mi>d</mi> <mi>&Omega;</mi> </mrow>

<mrow> <msub> <mi>W</mi> <mi>g</mi> </msub> <mo>=</mo> <msub> <mo>&Integral;</mo> <mi>&Omega;</mi> </msub> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <mfrac> <mrow> <mo>&part;</mo> <mi>w</mi> </mrow> <mrow> <mo>&part;</mo> <mi>x</mi> </mrow> </mfrac> </mtd> <mtd> <mfrac> <mrow> <mo>&part;</mo> <mi>w</mi> </mrow> <mrow> <mo>&part;</mo> <mi>y</mi> </mrow> </mfrac> </mtd> </mtr> </mtable> </mfenced> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <mrow> <msub> <mi>&gamma;</mi> <mn>1</mn> </msub> <msubsup> <mi>N</mi> <mi>x</mi> <mn>0</mn> </msubsup> </mrow> </mtd> <mtd> <mn>0</mn> </mtd> </mtr> <mtr> <mtd> <mn>0</mn> </mtd> <mtd> <mrow> <msub> <mi>&gamma;</mi> <mn>2</mn> </msub> <msubsup> <mi>N</mi> <mi>y</mi> <mn>0</mn> </msubsup> </mrow> </mtd> </mtr> </mtable> </mfenced> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <mfrac> <mrow> <mo>&part;</mo> <mi>w</mi> </mrow> <mrow> <mo>&part;</mo> <mi>x</mi> </mrow> </mfrac> </mtd> </mtr> <mtr> <mtd> <mfrac> <mrow> <mo>&part;</mo> <mi>w</mi> </mrow> <mrow> <mo>&part;</mo> <mi>y</mi> </mrow> </mfrac> </mtd> </mtr> </mtable> </mfenced> <mi>d</mi> <mi>&Omega;</mi> </mrow>

Π＝U_ε-W_g

(K+λK_g)u＝0

<mrow> <mover> <mi>k</mi> <mo>&OverBar;</mo> </mover> <mo>=</mo> <msup> <mi>&lambda;a</mi> <mn>2</mn> </msup> <mo>/</mo> <mrow> <mo>(</mo> <msub> <mi>E</mi> <mn>2</mn> </msub> <msup> <mi>h</mi> <mn>3</mn> </msup> <mo>)</mo> </mrow> </mrow>

2.根据权利要求1所述的基于无网格法的FG-GRC板屈曲载荷因子的数值算法，其特征在于：上述步骤1中，复合材料板的几何参数包括长a、宽b和高h；复合材料层合板由石墨烯与基体材料构成，确定复合材料层合板中石墨烯与基体材料的材料参数；和均为石墨烯的杨氏模量，G^G为石墨烯的剪切模量，ν^G为石墨烯的泊松比；E^m为基体材料的杨氏模量，G^m为基体材料的剪切模量，ν^m是基体材料的泊松比；并采用扩展Halpin-Tsai模型来计算复合材料层合板的杨氏模量E₁₁、E₂₂和剪切模量G₁₂：

<mrow> <msub> <mi>E</mi> <mn>11</mn> </msub> <mo>=</mo> <msub> <mi>&eta;</mi> <mn>1</mn> </msub> <mfrac> <mrow> <mn>1</mn> <mo>+</mo> <mrow> <mo>(</mo> <mn>2</mn> <msub> <mi>a</mi> <mi>G</mi> </msub> <mo>/</mo> <msub> <mi>h</mi> <mi>G</mi> </msub> <mo>)</mo> </mrow> <msub> <mi>&gamma;</mi> <mn>11</mn> </msub> <msub> <mi>V</mi> <mi>G</mi> </msub> </mrow> <mrow> <mn>1</mn> <mo>-</mo> <msub> <mi>&gamma;</mi> <mn>11</mn> </msub> <msub> <mi>V</mi> <mi>G</mi> </msub> </mrow> </mfrac> <msup> <mi>E</mi> <mi>m</mi> </msup> </mrow>

<mrow> <msub> <mi>E</mi> <mn>22</mn> </msub> <mo>=</mo> <msub> <mi>&eta;</mi> <mn>2</mn> </msub> <mfrac> <mrow> <mn>1</mn> <mo>+</mo> <mrow> <mo>(</mo> <mn>2</mn> <msub> <mi>b</mi> <mi>G</mi> </msub> <mo>/</mo> <msub> <mi>h</mi> <mi>G</mi> </msub> <mo>)</mo> </mrow> <msub> <mi>&gamma;</mi> <mn>22</mn> </msub> <msub> <mi>V</mi> <mi>G</mi> </msub> </mrow> <mrow> <mn>1</mn> <mo>-</mo> <msub> <mi>&gamma;</mi> <mn>22</mn> </msub> <msub> <mi>V</mi> <mi>G</mi> </msub> </mrow> </mfrac> <msup> <mi>E</mi> <mi>m</mi> </msup> </mrow>

<mrow> <msub> <mi>G</mi> <mn>12</mn> </msub> <mo>=</mo> <msub> <mi>&eta;</mi> <mn>3</mn> </msub> <mfrac> <mn>1</mn> <mrow> <mn>1</mn> <mo>-</mo> <msub> <mi>&gamma;</mi> <mn>12</mn> </msub> <msub> <mi>V</mi> <mi>G</mi> </msub> </mrow> </mfrac> <msup> <mi>G</mi> <mi>m</mi> </msup> </mrow>

中间变量γ₁₁、γ₂₂和γ₁₂分别为：

<mrow> <msub> <mi>&gamma;</mi> <mn>11</mn> </msub> <mo>=</mo> <mfrac> <mrow> <msubsup> <mi>E</mi> <mn>11</mn> <mi>G</mi> </msubsup> <mo>/</mo> <msup> <mi>E</mi> <mi>m</mi> </msup> <mo>-</mo> <mn>1</mn> </mrow> <mrow> <msubsup> <mi>E</mi> <mn>11</mn> <mi>G</mi> </msubsup> <mo>/</mo> <msup> <mi>E</mi> <mi>m</mi> </msup> <mo>+</mo> <mn>2</mn> <msub> <mi>a</mi> <mi>G</mi> </msub> <mo>/</mo> <msub> <mi>h</mi> <mi>G</mi> </msub> </mrow> </mfrac> </mrow>

<mrow> <msub> <mi>&gamma;</mi> <mn>22</mn> </msub> <mo>=</mo> <mfrac> <mrow> <msubsup> <mi>E</mi> <mn>22</mn> <mi>G</mi> </msubsup> <mo>/</mo> <msup> <mi>E</mi> <mi>m</mi> </msup> <mo>-</mo> <mn>1</mn> </mrow> <mrow> <msubsup> <mi>E</mi> <mn>22</mn> <mi>G</mi> </msubsup> <mo>/</mo> <msup> <mi>E</mi> <mi>m</mi> </msup> <mo>+</mo> <mn>2</mn> <msub> <mi>b</mi> <mi>G</mi> </msub> <mo>/</mo> <msub> <mi>h</mi> <mi>G</mi> </msub> </mrow> </mfrac> </mrow>

<mrow> <msub> <mi>&gamma;</mi> <mn>12</mn> </msub> <mo>=</mo> <mfrac> <mrow> <msubsup> <mi>G</mi> <mn>12</mn> <mi>G</mi> </msubsup> <mo>/</mo> <msup> <mi>G</mi> <mi>m</mi> </msup> <mo>-</mo> <mn>1</mn> </mrow> <mrow> <msubsup> <mi>G</mi> <mn>12</mn> <mi>G</mi> </msubsup> <mo>/</mo> <msup> <mi>G</mi> <mi>m</mi> </msup> </mrow> </mfrac> </mrow>

其中，V_G为石墨烯的体积分数，V_m＝1-V_G为基体材料的体积分数；温度的变化对泊松比的影响很小，复合材料层合板的泊松比ν₁₂为：

ν₁₂＝V_Gν^G+V_mν^m

根据Schapery模型，得出复合材料在纵向和横向上的热膨胀系数；

<mrow> <msub> <mi>&alpha;</mi> <mn>11</mn> </msub> <mo>=</mo> <mfrac> <mrow> <msub> <mi>V</mi> <mi>G</mi> </msub> <msubsup> <mi>E</mi> <mn>11</mn> <mi>G</mi> </msubsup> <msubsup> <mi>&alpha;</mi> <mn>11</mn> <mi>G</mi> </msubsup> <mo>+</mo> <msub> <mi>V</mi> <mi>m</mi> </msub> <msup> <mi>E</mi> <mi>m</mi> </msup> <msup> <mi>&alpha;</mi> <mi>m</mi> </msup> </mrow> <mrow> <msub> <mi>V</mi> <mi>G</mi> </msub> <msubsup> <mi>E</mi> <mn>11</mn> <mi>G</mi> </msubsup> <mo>+</mo> <msub> <mi>V</mi> <mi>m</mi> </msub> <msup> <mi>E</mi> <mi>m</mi> </msup> </mrow> </mfrac> </mrow>

<mrow> <msub> <mi>&alpha;</mi> <mn>22</mn> </msub> <mo>=</mo> <mrow> <mo>(</mo> <mn>1</mn> <mo>+</mo> <msup> <mi>v</mi> <mi>G</mi> </msup> <mo>)</mo> </mrow> <msub> <mi>V</mi> <mi>G</mi> </msub> <msubsup> <mi>&alpha;</mi> <mn>22</mn> <mi>G</mi> </msubsup> <mo>+</mo> <mrow> <mo>(</mo> <mn>1</mn> <mo>+</mo> <msup> <mi>v</mi> <mi>m</mi> </msup> <mo>)</mo> </mrow> <msup> <mi>V</mi> <mi>m</mi> </msup> <msup> <mi>&alpha;</mi> <mi>m</mi> </msup> <mo>-</mo> <msub> <mi>v</mi> <mn>12</mn> </msub> <msub> <mi>&alpha;</mi> <mn>11</mn> </msub> </mrow>

3.根据权利要求1所述的基于无网格法的FG-GRC板屈曲载荷因子的数值算法，其特征在于，步骤4中，基于再生核质点法，得到的形函数ψ_I(x)为：

ψ_I(x)＝C(x；x-x_I)Φ_a(x-x_I)

其中Φ_a(x-x_I)为满足再生条件的核函数，C(x；x-x_I)为满足再生条件的修正函数；

<mfenced open = "" close = ""> <mtable> <mtr> <mtd> <mrow> <msubsup> <mo>&Sigma;</mo> <mrow> <mi>I</mi> <mo>-</mo> <mn>1</mn> </mrow> <mrow> <mi>N</mi> <mi>P</mi> </mrow> </msubsup> <msub> <mi>&Psi;</mi> <mi>I</mi> </msub> <mrow> <mo>(</mo> <mi>x</mi> <mo>)</mo> </mrow> <msubsup> <mi>x</mi> <mi>I</mi> <mi>p</mi> </msubsup> <msubsup> <mi>y</mi> <mi>I</mi> <mi>q</mi> </msubsup> <mo>=</mo> <msup> <mi>x</mi> <mi>p</mi> </msup> <msup> <mi>y</mi> <mi>q</mi> </msup> </mrow> </mtd> <mtd> <mrow> <mi>f</mi> <mi>o</mi> <mi>r</mi> </mrow> </mtd> <mtd> <mrow> <mo>(</mo> <mi>p</mi> <mo>+</mo> <mi>q</mi> <mo>=</mo> <mn>0</mn> <mo>,</mo> <mn>1</mn> <mo>,</mo> <mn>2</mn> <mo>)</mo> </mrow> </mtd> </mtr> </mtable> </mfenced>

其中NP为离散节点数量，修正函数C(x；x-x_I)是由多项式基函数的线性组合表示而成；

C(x；x-x_I)＝H^T(x-x_I)b(x)

其中

b(x)＝[b₀(x,y),b₁(x,y),b₂(x,y),b₃(x,y),b₄(x,y),b₅(x,y)]^T

H^T(x-x_I)＝[1,x-x_I,y-y_I,(x-x_I)(y-y_I),(x-x_I)²,(y-y_I)²]

形函数表示为：

Ψ_I(x)＝b^T(x)H(x-x_I)Φ_a(x-x_I)

b(x)为：

b(x)＝M^-1(x)H(0)

其中

<mrow> <mi>M</mi> <mrow> <mo>(</mo> <mi>x</mi> <mo>)</mo> </mrow> <mo>=</mo> <msubsup> <mo>&Sigma;</mo> <mrow> <mi>I</mi> <mo>=</mo> <mn>1</mn> </mrow> <mrow> <mi>N</mi> <mi>P</mi> </mrow> </msubsup> <mi>H</mi> <mrow> <mo>(</mo> <mi>x</mi> <mo>-</mo> <msub> <mi>x</mi> <mi>I</mi> </msub> <mo>)</mo> </mrow> <msup> <mi>H</mi> <mi>T</mi> </msup> <mrow> <mo>(</mo> <mi>x</mi> <mo>-</mo> <msub> <mi>x</mi> <mi>I</mi> </msub> <mo>)</mo> </mrow> <msub> <mi>&Phi;</mi> <mi>a</mi> </msub> <mrow> <mo>(</mo> <mi>x</mi> <mo>-</mo> <msub> <mi>x</mi> <mi>I</mi> </msub> <mo>)</mo> </mrow> </mrow>

H(0)＝[1,0,0,0,0,0,]^T

Φ_a(x-x_I)是二维核函数：

Φ_a(x-x_I)＝Φ_a(x)Φ_a(y)

其中Φ_a(x)和Φ_a(y)为中间变量；

选用了三次样条函数作为权函数：

其中

其中d_max表示的是比例因子，d_I表示的是扩张系数，为了避免矩阵M的奇异性，在选择距离c_I的时候必须遍历足够多的节点；(x_I，y_I，z_I)为节点I的坐标；

将形函数改写为：

Ψ_I(x)＝H^T(0)M^-1(x)H(x-x_I)Φ_a(x-x_I)

M(x)b(x)＝H(0)

对矩阵M(x)进行LU分解用来计算b(x)，得出

M_x(x)b(x)+M(x)b_x(x)＝H_x(0)

M(x)b_x(x)＝H_x(0)-M_x(x)b(x)

其中z_I、b(x)、M(x)、M(x)和H^T(x-x_I)均为中间变量，H^T(x-x_I)和H(x-x_I)互为转置矩阵。

4.根据权利要求1所述的基于无网格法的FG-GRC板屈曲载荷因子的数值算法，其特征在于，步骤5中，通过完全变换法来施加本质边界条件使得无网格方法的形函数满足Kronecker delta属性：

广义位移u表示为：

<mrow> <msub> <mover> <mi>u</mi> <mo>~</mo> </mover> <mi>J</mi> </msub> <mo>=</mo> <msubsup> <mo>&Sigma;</mo> <mrow> <mi>I</mi> <mo>=</mo> <mn>1</mn> </mrow> <mrow> <mi>N</mi> <mi>P</mi> </mrow> </msubsup> <msub> <mi>&Psi;</mi> <mi>I</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mi>J</mi> </msub> <mo>)</mo> </mrow> <msub> <mi>u</mi> <mi>I</mi> </msub> <mo>=</mo> <msubsup> <mo>&Sigma;</mo> <mrow> <mi>I</mi> <mo>=</mo> <mn>1</mn> </mrow> <mrow> <mi>N</mi> <mi>P</mi> </mrow> </msubsup> <msub> <mi>L</mi> <mrow> <mi>I</mi> <mi>J</mi> </mrow> </msub> <msub> <mi>u</mi> <mi>I</mi> </msub> </mrow>

其中是节点x_J处的位移，u_I为节点参数，L_IJ＝Ψ_I(x_J)为形函数；

<mrow> <msub> <mi>u</mi> <mi>I</mi> </msub> <mo>=</mo> <msubsup> <mo>&Sigma;</mo> <mrow> <mi>K</mi> <mo>-</mo> <mn>1</mn> </mrow> <mrow> <mi>N</mi> <mi>P</mi> </mrow> </msubsup> <msubsup> <mi>L</mi> <mrow> <mi>I</mi> <mi>K</mi> </mrow> <mrow> <mo>-</mo> <mi>T</mi> </mrow> </msubsup> <msub> <mover> <mi>u</mi> <mo>~</mo> </mover> <mi>I</mi> </msub> </mrow>

由上述两个公式推导出：

其中

为满足Kronecker delta属性的形函数，δ_IJ为中间变量。