A ROBUST ISOGEOMETRIC MESHFREE COLLOCATION METHOD FOR TARGETED REFINEMENT ANALYSIS
-
摘要: 现有等几何基函数的无网格表示理论虽然搭建起了无网格法和等几何分析方法的连接桥梁, 但该方法并未解决高阶梯度计算复杂且耗时的问题. 文章在无网格形函数再生梯度理论基础上, 构建了由等几何基函数再生点定义的混合梯度基向量, 提出了一种等几何基函数梯度的无网格等价表述形式, 丰富和完善了无网格法与等几何分析的无缝衔接. 在此基础上, 提出了一种鲁棒的等几何无网格配点法, 该方法保证了全域内近似函数的一致性、高效性和模型局部靶向精细化的简捷性. 与传统等几何基函数梯度的递推算法相比, 所提方法可同步构造相应阶次等几何基函数梯度, 数值实现更加简捷. 文中以二维二次基函数为例, 数值验证了特定影响域条件下所构造的无网格再生梯度与等几何基函数的标准梯度之间的等价性. 同时对比分析了所提无网格再生梯度与等几何无网格形函数直接梯度的再生条件, 结果表明, 等几何无网格形函数直接梯度的再生条件并不稳定, 而所提等几何基函数的无网格再生梯度在不同节点序列精确满足不同阶次的再生条件, 保证了配点法分析的结果收敛性. 数值算例结果表明, 所提配点方法相较于传统等几何配点法在局部细化问题分析方面具有更高的计算精度.Abstract: The reproducing kernel meshfree representation of isogeometric basis functions provides a meaningful way to link the meshfree methods and isogeometric analysis, while the costly gradient computation of isogeometric basis functions has not been well addressed. In this work, based upon the reproducing kernel gradient formulation of meshfree shape functions, a mixed gradient reproducing basis vector is defined by the reproducing points of isogeometric basis functions. Subsequently, a meshfree framework is developed to construct the gradients of isogeometric basis functions in a meshfree fashion, which further unifies the meshfree methods and isogeometric analysis. Under this framework, a robust isogeometric meshfree collocation method is proposed in this study, which ensures the consistency and efficiency of the approximation function in the whole domain and the flexibility of the local model refinement. Unlike the recursive construction process of isogeometric basis functions, the proposed approach is a direct and one-step procedure for the gradient evaluation, which is trivial for numerical implementation. It turns out that the proposed meshfree formulation of isogeometric gradients yields the same gradients as the standard gradients of isogeometric basis functions under certain conditions regarding the influence domains of meshfree shape functions. Meanwhile, the reproducing conditions of the direct gradients of the isogeometric basis functions are unstable, but the reconstructed meshfree reproducing gradients meet complete the gradient reproducing conditions, which ensures the convergence of the collocation methods. Numerical results demonstrate that the proposed collocation methodology exhibits much higher computational accuracy than the conventional isogeometric collocation in the analysis of local refinement problems.
-
引 言
等几何分析方法[1] (isogeometric analysis, IGA)实现了有限元分析与计算机辅助设计(CAD)的无缝衔接, 该方法利用非均匀有理B样条(non-uniform rational B-spline, NURBS)作为有限元分析的形函数, 具有几何精确和高阶光滑的特性. 与采用C0连续性形函数的传统有限元法相比, 表现出更高计算精度. 因此, 在解决薄板壳高阶问题[2-3]、流固耦合问题[4]、大变形分析[5]和裂纹扩展模拟[6]等复杂工程问题时展现出广阔的应用前景. 此外, 由于等几何分析和无网格法[7-11]采用的基函数或形函数均具有高阶光滑特性, 使得需要计算基函数或形函数高阶梯度的配点法得到了快速发展. 例如最小二乘与再生核无网格配点法[12-13]、径向基无网格配点法[14-16]、无网格稳定配点法[17-18]和等几何配点法[19-22]等. 总的来说, 配点法具有构造简单、计算效率高以及边界条件易施加的优点, 近年来逐渐被广泛应用于各科学工程领域.
在实际工程应用中, 时常需要通过模型局部细化来提高关键位置的数值解精度, 但是基于张量积运算建立多维基函数的等几何分析方法在进行模型局部细化时, 会导致相关方向区域的网格被加密, 大幅增加额外工作量. 图1给出了1/4圆环曲面采用等几何分析方法局部细化结果, 从图中可以看出, 在初始细化单元基础上进行域内局部细化时, 相邻区域也得到了非必要细化, 增加了额外自由度数量.
![]()
图 1 1/4圆环曲面的等几何局部细化Figure 1. Local refinement for a quarter of the annulus in isogeometric analysis鉴于上述等几何分析局部细化的缺陷, Sederberg等[23]提出了T样条等几何分析方法. 王悦等[24]基于T样条单元的局部细化方法, 有效提升了薄板动力学分析的精度和计算效率. 此外, Bazilevs等[25]、Scott等[26]、Chen等 [27]和Casquero等[28]诸多学者对T样条等几何分析方法的改进和应用做了许多研究工作. 同时, 分层B样条[29]、PHT样条[30]、LR样条[31]和THB样条[32]也被应用于提升等几何局部细化过程的计算效率, 但上述细化方法中节点的插入形式较为复杂. 为了更有效地实现局部精细化处理, Zhang等[33]构建了等几何无网格形函数统一表示理论, 使得局部节点的加密操作变得更为简捷. 基于该理论框架, Li等[34]发展了一种弹性接触断裂问题的等几何无网格配点法, 验证了该理论局部细化的高效性. Huang等[35]将等几何无网格配点法推广到多层复合板的静力、自由振动和屈曲问题分析, 揭示了等几何无网格配点法良好的计算精度. 需要指出的是, 该理论仅建立了等几何基函数与无网格形函数的内在联系, 并未解决形函数高阶梯度计算复杂且耗时的问题. 同时, 等几何无网格形函数中矩量矩阵的逆矩阵求导运算对配点法稳定性影响也是一个待研究问题.
为了规避形函数高阶梯度计算时矩量矩阵的逆矩阵求导复杂运算问题, 提升配点法中高阶梯度计算效率, 本文在无网格法再生梯度理论[36]框架下, 通过分析无网格形函数与等几何基函数再生条件, 提出了由等几何基函数再生点定义的混合梯度基向量, 建立了无网格形函数和等几何基函数两者之间梯度的内在联系和等效性. 与等几何无网格形函数直接梯度[34-35]不同, 所构造的无网格再生等几何基函数梯度, 满足稳定且完备的再生多项式条件, 保证了配点法分析的收敛性. 在此基础上, 文中提出了一种鲁棒的等几何无网格配点法. 所提方法兼具了等几何分析几何精确的特性和无网格法灵活建模的优势, 可以实现局部特征区域多层细化, 在保证计算精度的前提下自由控制实际问题离散过程中非必要自由度的数量, 从而有效降低计算成本. 最后, 文中通过系列算例验证了所提方法的收敛精度.
1. 等几何配点法离散方程
本文考虑场变量为$ {\boldsymbol{u}}({\boldsymbol{x}}) $的二阶边值问题, 其控制方程的强形式可表示为
$$\left.\begin{split} &\mathcal{L}{\boldsymbol{u}}({\boldsymbol{x}}) + {\boldsymbol{b}}({\boldsymbol{x}}) = {\boldsymbol{0}},\quad {\boldsymbol{x}}\in \varOmega \\ &\mathcal{B}{\boldsymbol{u}}({\boldsymbol{x}}) = {\boldsymbol{g}}({\boldsymbol{x}}),\quad {\boldsymbol{x}}\in \varGamma \end{split} \right\}$$ (1) 式中, $ \mathcal{L} $表示二阶变分算子, 在稳态势问题和弹性力学问题中, $ \mathcal{L} $分别为拉普拉斯算子和弹性变分算子; ${\boldsymbol{b}}({\boldsymbol{x}})$为体力项或源项; $ \mathcal{B} $代表边界算子, $ {\boldsymbol{g}}({\boldsymbol{x}}) $为预先给定的$ \varGamma $边界值.
在等几何分析中, 通常采用具有几何精确性的NURBS基函数$ {R_A}({{\boldsymbol{\xi}} }) $来离散场变量$ {\boldsymbol{u}}({\boldsymbol{x}}) $, 相应近似解$ {{\boldsymbol{u}}^h}({\boldsymbol{x}}) $表示为
$$ {{\boldsymbol{u}}^h}({\boldsymbol{x}}){\text{ = }}\sum\limits_{A = 1}^{{NP} } {R_A^{}({{\boldsymbol{\xi}} })} {{\boldsymbol{d}}_A} $$ (2) 几何物理构型由NURBS基函数描述为
$$ {\boldsymbol{x}}({{\boldsymbol{\xi}} }){\text{ = }}\sum\limits_{A = 1}^{{NP} } {R_A^{}({{\boldsymbol{\xi}} })} {{\boldsymbol{P}}_A} $$ (3) 式中, $ {{\{}}{{\boldsymbol{P}}_A}{{\} }}_A^{{NP} } $为几何描述的控制点, NP为控制点的个数, $ {{\boldsymbol{d}}_A} $为控制点$ {{\boldsymbol{P}}_A} $所对应的场变量系数.
在等几何配点法框架下, 将场变量近似式(2)和几何构型描述式(3)代入控制方程式(1), 可得到如下离散配点方程
$$ {\boldsymbol{Kd}} = {\boldsymbol{f}} $$ (4) 其中, $ {\boldsymbol{K}} $, $ {\boldsymbol{d}} $和$ {\boldsymbol{f}} $分别表示刚度矩阵、位移系数列向量和力列向量, 有
$$ {\boldsymbol{K}} = \left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{K}}^i}} \\ {{{\boldsymbol{K}}^b}} \end{array}} \right],\quad {\text{ }}{\boldsymbol{f}} = \left\{ {\begin{array}{*{20}{c}} {{{\boldsymbol{f}}^i}} \\ {{{\boldsymbol{f}}^b}} \end{array}} \right\} $$ (5) $$ \left. \begin{aligned} & {{\boldsymbol{K}}_{AB}^i = \mathcal{L}[{R_B}({{\boldsymbol{x}}_A})],\;\;{\boldsymbol{f}}_A^i = - {\boldsymbol{b}}({{\boldsymbol{x}}_A}),}\;\; {{{\boldsymbol{x}}_A} \in \varOmega } \\ & {{\boldsymbol{K}}_{AB}^b = \mathcal{B}[{R_B}({{\boldsymbol{x}}_A})],\;\;{\boldsymbol{f}}_A^b = {\boldsymbol{g}}({{\boldsymbol{x}}_A}),}\;\; {{{\boldsymbol{x}}_A} \in \varGamma } \end{aligned} \right\} $$ (6) 式中$ {{\boldsymbol{x}}_A} = {\boldsymbol{x}}({{{\boldsymbol{\xi}} }_A}) $为参数空间点$ {{{\boldsymbol{\xi}} }_A} $在物理空间中的投影点, $ {{{\boldsymbol{\xi}} }_A} $代表等几何参数空间上配点的位置.
值得注意的是, NURBS基函数是由建立在参数空间中B样条基函数加权得到. 因此, NURBS基函数$ {R_A}({{\boldsymbol{\xi}} }) $的物理空间导数与等参有限元求解相类似, 需要引入物理空间与参数空间变换的雅克比矩阵求得. 在二维情况下, NURBS基函数的空间导数为
$$ {\mathcal{G}_{\boldsymbol{x}}}[R_A^{}] = {{\boldsymbol{J}}^{ - {\mathrm{T}}}}{\mathcal{G}_{\xi }}[R_A^{}] $$ (7) $$ {\mathcal{G}_{{\boldsymbol{xx}}}}[R_A^{}] = {{\boldsymbol{H}}^{ - 1}}({\mathcal{G}_{{\xi \xi }}}[R_A^{}] - {{\boldsymbol{L}}^{\mathrm{T}}}{\mathcal{G}_{\boldsymbol{x}}}[R_A^{}]) $$ (8) 其中, $ {\mathcal{G}_{\xi }}[{R_{{A}}}] $和$ {\mathcal{G}_{{\xi \xi }}}[{R_A}] $分别表示$ {R_A}({{\boldsymbol{\xi}} }) $的参数空间一阶导数和二阶导数, $ {\mathcal{G}_{\boldsymbol{x}}}[{R_{{A}}}] $和$ {\mathcal{G}_{{\boldsymbol{xx}}}}[{R_A}] $分别表示$ {R_A}({{\boldsymbol{\xi}} }) $的物理空间一阶导数和二阶导数, 有
$$ {\mathcal{G}_{\xi }}[R_A^{}] = {\{ \begin{array}{*{20}{c}} {R_{A,\xi }^{}}&{R_{A,\eta }^{}} \end{array}\} ^{\mathrm{T}}} $$ (9) $$ {\mathcal{G}_{\boldsymbol{x}}}[R_A^{}] = {\{ \begin{array}{*{20}{c}} {R_{A,x}^{}}&{R_{A,y}^{}} \end{array}\} ^{\mathrm{T}}} $$ (10) $$ {\mathcal{G}_{{\xi \xi }}}[R_A^{}] = {\{ \begin{array}{*{20}{c}} {R_{A,\xi \xi }^{}}&{R_{A,\eta \eta }^{}}&{R_{A,\xi \eta }^{}} \end{array}\} ^{\mathrm{T}}} $$ (11) $$ {\mathcal{G}_{{\boldsymbol{xx}}}}[R_A^{}] = {\{ \begin{array}{*{20}{c}} {R_{A,xx}^{}}&{R_{A,yy}^{}}&{2R_{A,xy}^{}} \end{array}\} ^{\mathrm{T}}} $$ (12) 相应的雅克比矩阵${\boldsymbol{J}}$, ${\boldsymbol{L}}$和${\boldsymbol{H}}$分别为
$$ {\boldsymbol{J}} = \left[ {\begin{array}{*{20}{c}} {{x_{,\xi }}}&{{x_{,\eta }}} \\ {{y_{,\xi }}}&{{y_{,\eta }}} \end{array}} \right],\quad {\boldsymbol{L}} = \left[ {\begin{array}{*{20}{l}} {{x_{,\xi \xi }}}&{{x_{,\eta \eta }}}&{{x_{,\xi \eta }}} \\ {{y_{,\xi \xi }}}&{{y_{,\eta \eta }}}&{{y_{,\xi \eta }}} \end{array}} \right] $$ (13) $$ {\boldsymbol{H}} = \left[ {\begin{array}{*{20}{c}} {x_{,\xi }^2}&{y_{,\xi }^2}&{{x_{,\xi }}{y_{,\xi }}} \\ {x_{,\eta }^2}&{y_{,\eta }^2}&{{x_{,\eta }}{y_{,\eta }}} \\ {{x_{,\xi }}{x_{,\eta }}}&{{y_{,\xi }}{y_{,\eta }}}&{({x_{,\xi }}{y_{,\eta }} + {x_{,\eta }}{y_{,\xi }})/2} \end{array}} \right] $$ (14) 式(13)和式(14)雅克比矩阵均可由式(3)计算得到.
2. 等几何分析的再生核无网格表示理论
2.1 等几何基函数的无网格表述
为了便于分析等几何B样条基函数的无网格表述形式, 可将无网格形函数在参数空间进行构造. 基于再生核无网格近似理论[8], 无网格形函数的多项式再生条件为
$$ \sum\limits_A {{\varPsi _A}({{\boldsymbol{\xi}} })} \xi _a^i\eta _b^j = {\xi ^i}{\eta ^j} $$ (15) 其中, $ \xi $和$ \eta $表示参数空间维度的方向. 将式(15)改写成向量形式, 有
$$ \sum\limits_A {{\varPsi _A}({{\boldsymbol{\xi}} })} {\boldsymbol{p}}({{\boldsymbol{\xi }}_A}) = {\boldsymbol{p}}({\boldsymbol{\xi}} ) $$ (16) 其中, $ {\boldsymbol{p}}({\boldsymbol{\xi}} ) $为基向量. 二维情况下, $p$次单项式基向量为
$$ {\boldsymbol{p}}({\boldsymbol{\xi}} ) = \{ {\psi _{00}},{\psi _{10}},{\psi _{01}},{\psi _{20}},{\psi _{11}},{\psi _{02}}, \cdots , {\psi _{ij}}, \cdots ,{\psi _{0p}}{\} ^{\mathrm{T}}} $$ (17) 式中, ${\psi _{ij}} = {\xi ^i}{\eta ^j}$, $0 \leqslant i + j \leqslant p$.
由式(16)可见, 无网格形函数在固定离散点满足完备再生条件. 与无网格形函数的再生条件不同, 等几何B样条基函数的再生条件依赖于各阶再生点序列[37], 有
$$ {\displaystyle \sum _{A}{N}_{A}^{p}({\boldsymbol{\xi}} ){\boldsymbol{p}}({{\boldsymbol{\xi}} }_{A}^{[\cdot]}) = {\boldsymbol{p}}({{\boldsymbol{\xi}} }^{[\cdot]})} $$ (18) 其中, $ {{\boldsymbol{\xi}} }_{A}^{[\cdot]} $为B样条再生点, ${\boldsymbol{ p}}({{\boldsymbol{\xi}} }_{A}^{[\cdot]}) $为由B样条再生点定义的混合基向量. 例如在二维情况下, $p$次混合基向量定义为
$$\begin{split} & \boldsymbol{p}({\boldsymbol{\xi}} _{A}^{[\cdot ]})\text{=}\Big\{1,\xi _{a}^{[1]},\eta _{b}^{[1]},\xi _{a}^{[1]}\eta _{b}^{[1]},{{(\xi _{a}^{[2]})}^{2}},{{(\eta _{b}^{[2]})}^{2}}, \\ &\qquad \xi _{a}^{[1]}{{(\eta _{b}^{[2]})}^{2}},{{(\xi _{a}^{[2]})}^{2}}\eta _{b}^{[1]},{{(\xi _{a}^{[2]})}^{2}}{{(\eta _{b}^{[2]})}^{2}}, \cdots ,\\ &\qquad {{(\xi _{a}^{[p]})}^{p}},{{(\eta _{b}^{[p]})}^{p}},\cdots ,{{(\xi _{a}^{[p]})}^{p}}{{(\eta _{b}^{[p]})}^{p}}{{\Big\}}^{{\mathrm{T}}}} \end{split}$$ (19) 在传统再生核无网格近似中, 式(19)中再生点$ {{\boldsymbol{\xi}} }_{A}^{[\cdot]} $为无网格节点, 即$ {{\boldsymbol{\xi}} }_{A}^{}\text{ = }{{\boldsymbol{\xi}} }_{A}^{[\cdot]} $. 在等几何分析中, 多维情况下的等几何基函数再生点可由一维再生点通过张量积运算获得, 如二维再生点$ {{\boldsymbol{\xi}} }_{A}^{[\cdot]}\text{ = (}{\xi }_{a}^{[\cdot]}\text{, }{\eta }_{b}^{[\cdot]}\text{)} $可由一维$\xi $和$\eta $方向的相应阶次再生点的张量积形式得到, 一维情况下l阶再生点${(\xi _a^{[l]})^l}$的表达式为[37]
$$ \xi _a^{[l]} = \sqrt[l]{{S_p^l[G_{a + 1}^{a + p}]/C_p^l}},\;\;C_p^l = p!/[l!(p - l)!] $$ (20) 其中, $ G_{a + 1}^{a + p} = [a + 1,a + 2, \cdots ,a + p] $, $ S_p^l[G_{a + 1}^{a + p}] $表示集合$ G_{a + 1}^{a + p} $中任意l个元素相乘项的求和. 等几何配点法通常采用Greville点作为配点位置, 而Greville点本质上等价于一阶再生点. 因此, 二维情况下Greville配点$ {{{\boldsymbol{\xi}} }_A} $可根据式(20)获得, 有
$$ {{\boldsymbol{\xi}} _A} = (\xi _a^{[1]},\eta _b^{[1]}) = \frac{1}{p}\left(\sum\limits_{i = 1}^p {{\xi _{a + i}}} ,\sum\limits_{i = 1}^p {{\eta _{b + i}}} \right) $$ (21) 由式(16)和式(18)可知, 等几何基函数和无网格形函数均满足类似的多项式再生条件. 将再生点定义的混合基向量$ {\boldsymbol{p}}({{\boldsymbol{\xi}} }_{A}^{[\cdot]}) $引入到无网格形函数中, 便建立了等几何B样条基函数的再生核无网格形式[33]
$$ {\varPsi }_{A}^{[p]}({\boldsymbol{\xi}}) = {{\boldsymbol{p}}}^{{\mathrm{T}}}({{\boldsymbol{\xi}}}^{[\cdot]}){{\boldsymbol{M}}}^{-1}({\boldsymbol{\xi}}){\boldsymbol{p}}({{\boldsymbol{\xi}}}_{A}^{[\cdot]}){\varphi }_{s}({{\boldsymbol{\xi}}}_{A}^{[1]}-{\boldsymbol{\xi}}) $$ (22) 其中, 矩量矩阵${{\boldsymbol{M}}}({\boldsymbol{\xi}} )$为
$$ {{\boldsymbol{M}}}({\boldsymbol{\xi}}) = {\displaystyle \sum _{A}{\boldsymbol{p}}({{\boldsymbol{\xi}}}_{A}^{[\cdot]}){{\boldsymbol{p}}}^{{\mathrm{T}}}({{\boldsymbol{\xi}}}_{A}^{[\cdot]}){\varphi }_{s}({{\boldsymbol{\xi}}}_{A}^{[1]}-{\boldsymbol{\xi}})} $$ (23) 式中, 核函数${\varphi _s}({{\boldsymbol{\xi}} }_A^{[1]} - {{\boldsymbol{\xi}} })$是以一阶再生点${{\boldsymbol{\xi}} }_A^{[1]}$作为参照点, 进而确定影响域范围内任意计算点的权重值, 且一阶再生点作为无网格离散节点.
需要指出的是, 传统无网格形函数计算过程中必须保证矩量矩阵${ {\boldsymbol{M}}}({\boldsymbol{\xi}} )$的可逆性, 即所选核函数$ {\varphi _s}({{\boldsymbol{\xi}} }_A^{[1]} - {{\boldsymbol{\xi}} }) $的覆盖域应该足够大. 在一维情况下, 需保证不少于$(p + 1)$条形函数覆盖任意计算点${\boldsymbol{\xi}} $; 在多维情况下, 需保证不少于$ (p + {n_{sd}})!/(p!{n_{sd}}!) $条形函数覆盖任意区间内的计算点$ {\boldsymbol{\xi}} $, $ {n_{sd}} $表示求解问题的维数[8]. 而当核函数的影响域取为等几何基函数的覆盖域, 即$(p + 1)/2$时, 由式(22)可以精确表示相应阶次的等几何B样条基函数, 且这种等价关系只与B样条基函数的阶次$p$有关, 与核函数的选取无关[33].
在等几何分析中, 传统多维等几何B样条基函数可通过一维等几何基函数张量积运算得到, 这就导致了该方法在模型局部细化过程中会引入额外细化单元, 增加了计算成本. 而等几何B样条基函数采用再生核无网格形函数进行表述, 即等几何无网格形函数, 该方法兼具了等几何分析精确建模和无网格法易于实现计算模型局部细化的优势. 以二次混合基函数$(p = 2)$为例, 图2给出了基于初始离散模型的内部网格单元逐级靶向精细化处理过程, 其中红色点分布区域为局部细化区域, 从图中可以看出, 基于节点离散信息构建等几何无网格形函数, 使得模型局部细化变得更加简捷, 且减少了冗余自由度的计算, 提高了计算效率; 图3绘制了图2中指定细化区域内节点的等几何无网格形函数, 从图中可以直观看出, 等几何B样条基函数的再生核无网格形式具有准凸近似特性.
![]()
图 3 模型细化离散下的等几何无网格形函数Figure 3. The isogeometric meshfree shape functions under the model different level refinements此外, 基于等几何B样条基函数的无网格表示形式(22), 引入控制几何形状的权重因子, 便建立了用于等几何分析的NURBS基函数$ {R_A}({\boldsymbol{\xi}} ) $, 有
$$ {R_A}({\boldsymbol{\xi}} ) = \frac{{\varPsi _A^{[p]}({{\boldsymbol{\xi}} }){w_A}}}{{\displaystyle\sum\limits_{B = 1}^{{NP} } {\varPsi _B^{[p]}({{\boldsymbol{\xi}} }){w_B}} }} $$ (24) 其中, ${w_A}$和${w_B}$为权重因子.
2.2 等几何B样条基函数梯度的再生核无网格重构
在配点法分析中, 往往需要形函数或基函数的高阶梯度. 虽然等几何B样条基函数可以采用再生核无网格形式进行表示, 但采用无网格形函数形式直接计算其一阶梯度和二阶梯度时, 需要计算矩量矩阵的逆矩阵${{ {\boldsymbol{M}}}^{ - 1}}({\boldsymbol{\xi}} )$的导数, 增加了其高阶梯度计算难度 [13]. 在图2所示的逐层细化离散模型基础上, 图4给出了二次混合基函数下等几何无网格形函数直接梯度的再生条件数值验证结果, 其中黑色圆点为一阶再生点, 黑色菱形点为二阶再生点. 从图中不难发现, 等几何无网格形函数直接梯度不能精确满足相应再生条件, 即无法保证配点法分析结果的收敛性.
![]()
图 4 二次混合基函数下等几何无网格形函数直接梯度的再生条件数值验证Figure 4. Numerical verification of the direct gradient reproducing conditions for 2D isogeometric meshfree shape functions with quadratic mixed basis function针对逆矩阵${{ {\boldsymbol{M}}}^{ - 1}}({\boldsymbol{\xi}} )$求导运算所导致的等几何无网格形函数梯度的再生条件不稳定问题, 本文在无网格形函数再生梯度理论[36]基础上, 重构了等几何基函数梯度的无网格表述形式. 类似于无网格形函数的建立过程, 该重构形式基于形函数梯度再生条件, 定义了混合基向量$ {\boldsymbol{p}}({{\boldsymbol{\xi}} }_{}^{[\cdot]}) $的梯度基向量$ {{\boldsymbol{p}}}_{,\alpha }({{\boldsymbol{\xi}} }_{}^{[\cdot]}) $. 在构造上将式(22)中的基向量$ {\boldsymbol{p}}({{\boldsymbol{\xi}} }_{}^{[\cdot]}) $换为梯度基向量$ {{\boldsymbol{p}}}_{,\alpha }({{\boldsymbol{\xi}} }_{}^{[\cdot]}) $, 便可得到等几何分析的无网格再生梯度表达式
$$ {\varPsi }_{A,\alpha }^{[p]}({\boldsymbol{\xi}} )\text{ = }{{\boldsymbol{p}}}_{,\alpha }^{{\mathrm{T}}}({{\boldsymbol{\xi}}}^{[\cdot]}){{\boldsymbol{M}}}^{-1}({\boldsymbol{\xi}}){\boldsymbol{p}}({{\boldsymbol{\xi}}}_{A}^{[\cdot]}){\varphi }_{s}({{\boldsymbol{\xi}}}_{A}^{[1]}-{\boldsymbol{\xi}}) $$ (25) 其中, 多重指标符号$ {{\boldsymbol{\alpha }}} = {({\alpha _1},{\alpha _2}, \cdots ,{\alpha _{{n_{sd}}}})^{\text{T}}} $, $ \left| {{\boldsymbol{\alpha}} } \right|\; = \displaystyle\sum\limits_{i = 1}^{{n_{sd}}} {{\alpha _i}} $, $ \;\left| {\boldsymbol{\alpha}} \right| \leqslant p $; $ \varPsi _{A,\alpha }^{[p]}({\boldsymbol{\xi}} ) $为形函数的$ \left| {{\boldsymbol{\alpha }}} \right| $阶再生梯度; $ {{\boldsymbol{p}}}_{,\alpha }({{\boldsymbol{\xi}} }_{}^{[\cdot]}) $为梯度基向量, 二维一阶和二阶梯度基向量为
$$ \left.\begin{aligned} &{{\boldsymbol{p}}}_{,\xi }^{}({{\boldsymbol{\xi}} }_{}^{[\cdot]}) = \{0,1,0,\eta ,2\xi ,0,{\eta }^{2},2\xi \eta ,2\xi {\eta }^{2},\cdots ,\\ &\qquad p{\xi }^{p-1},0,\cdots ,p{\xi }^{p-1}{\eta }^{p}{\}}^{{\mathrm{T}}}\\ &{{\boldsymbol{p}}}_{,\eta }^{}({{\boldsymbol{\xi}} }_{}^{[\cdot]}) = \{0,0,1,\xi ,0,2\eta ,2\xi \eta ,{\xi }^{2},2{\xi }^{2}\eta ,\cdots ,\\ &\qquad 0,p{\eta }^{p-1},\cdots ,p{\xi }^{p}{\eta }^{p-1}{\}}^{{\mathrm{T}}}\end{aligned}\right\} $$ (26) $$ \left.\begin{aligned} &{{\boldsymbol{p}}}_{,\xi \xi }^{}({{\boldsymbol{\xi}} }^{[\cdot]}) = \{0,0,0,0,2,0,0,2\eta ,2{\eta }^{2},\cdots ,\\ &\qquad p(p-1){\xi }^{p-2},0,\cdots ,p(p-1){\xi }^{p-2}{\eta }^{p}{\}}^{{\mathrm{T}}}\\ &{{\boldsymbol{p}}}_{,\eta \eta }^{}({{\boldsymbol{\xi}} }_{}^{[\cdot]}) = \{0,0,0,0,0,2,2\xi ,0,2{\xi }^{2},\cdots ,0,\\ &\qquad p(p-1){\eta }^{p-2},\cdots ,p(p-1){\xi }^{p}{\eta }^{p-2}{\}}^{{\mathrm{T}}}\\ & {{\boldsymbol{p}}}_{,\xi \eta }^{}({{\boldsymbol{\xi}} }_{}^{[\cdot]}) = \{0,0,0,1,0,0,2\eta ,2\xi ,4\xi \eta ,\cdots ,0,\\ &\qquad 0,\cdots ,{p}^{2}{\xi }^{p-1}{\eta }^{p-1}{\}}^{{\mathrm{T}}}\end{aligned}\right\} $$ (27) 在二维情况下, 将一阶和二阶基向量式(26)和式(27)直接代入式(25)便可获得相应阶次的等几何B样条基函数梯度, 数值实现更加简捷, 显著提高了梯度计算效率. 为了系统证明所提无网格再生梯度与标准等几何B样条基函数梯度间的等价性, 附录A给出详细的理论验证过程, 进一步说明了所提方法的正确性.
当无网格形函数的影响域与B样条基函数覆盖域一致时, 图5给出了二维情况下二次B样条基函数梯度的再生核无网格表述结果, 从图中可以看出, 采用梯度基向量构造的无网格再生梯度等价于相应B样条基函数梯度. 同时, 采用图2中逐层细化离散模型的节点布置形式, 图6绘制了二次基混合函数下无网格再生梯度的再生条件. 与图4等几何无网格形函数的直接梯度再生条件相比, 所提无网格再生梯度在不同节点序列精确满足不同阶次的多项式再生条件, 保证了配点法计算结果的收敛性.
![]()
图 5 二维二次B样条基函数梯度的再生核无网格表述Figure 5. Reproducing kernel meshfree representations of the gradients for 2D quadratic B-spline basis functions![]()
图 6 二维二次混合基函数下无网格再生梯度的再生条件数值验证Figure 6. Numerical verification of the reproducing conditions for 2D meshfree reproducing gradients with quadratic mixed basis function3. 数值算例
本小节通过若干典型局部高梯度问题及应力集中弹性力学问题验证所提方法的收敛精度, 其中“IGC” 表示传统等几何配点法(isogeometric collocation method), “IGAMC” 表示本文所提等几何无网格配点法(isogeometric meshfree collocation method). 为便于对比分析, 算例中的精度度量采用如下误差范数形式
$$ {L_2}{{\;{\mathrm{error}}}} = \sqrt {{{\displaystyle\sum\limits_{{{I}} = 1}^{{{NP}}} {({{\boldsymbol{u}}_{\text{I}}} - {\boldsymbol{u}}_{\text{I}}^h) \cdot ({{\boldsymbol{u}}_{\text{I}}} - {\boldsymbol{u}}_{\text{I}}^h)} }}\Bigg/{{\displaystyle\sum\limits_{{{I}} = 1}^{{{NP}}} {{\boldsymbol{u}}_{\text{I}}^{} \cdot {\boldsymbol{u}}_{\text{I}}^{}} }}} $$ (28) $$ {{{\mathrm{energy}}\;{\mathrm{error}}}} = \sqrt {{{\displaystyle\sum\limits_{{{I}} = 1}^{{{NP}}} {({{\boldsymbol{\varepsilon}} _{\text{I}}} - {\boldsymbol{\varepsilon}} _{\text{I}}^h):({{\boldsymbol{\sigma}} _{\text{I}}} - {\boldsymbol{\sigma}} _{\text{I}}^h)} }}\Bigg/{{\displaystyle\sum\limits_{{{I}} = 1}^{{{NP}}} {{{\boldsymbol{\varepsilon}} _{\text{I}}}:{{\boldsymbol{\sigma}} _{\text{I}}}} }}} $$ (29) $$ {H_1}{{\;{\mathrm{error}}}} = \sqrt {{{\displaystyle\sum\limits_{{{I}} = 1}^{{{NP}}} {\left[({{\boldsymbol{u}}_{\text{I}}} - {\boldsymbol{u}}_{\text{I}}^h) \cdot ({{\boldsymbol{u}}_{\text{I}}} - {\boldsymbol{u}}_{\text{I}}^h) + \displaystyle\sum\limits_{r = 1}^{{n_{sd}}} {({{\boldsymbol{u}}_{{\text{I}},{x_r}}} - {\boldsymbol{u}}_{{\text{I}},{x_r}}^h) \cdot ({{\boldsymbol{u}}_{{\text{I}},{x_r}}} - {\boldsymbol{u}}_{{\text{I}},{x_r}}^h)} \right]} }}\Bigg/{{\displaystyle\sum\limits_{{{I}} = 1}^{{{NP}}} {\left({\boldsymbol{u}}_{\text{I}}^{} \cdot {\boldsymbol{u}}_{\text{I}}^{} + \displaystyle\sum\limits_{r = 1}^{{n_{sd}}} {{{\boldsymbol{u}}_{{\text{I}},{x_r}}} \cdot {{\boldsymbol{u}}_{{\text{I}},{x_r}}}} \right)} }}} $$ (30) 式中, $ {{\boldsymbol{u}}_{\text{I}}} $和$ {\boldsymbol{u}}_{\text{I}}^h $表示场变量$ {\boldsymbol{u}} $的精确解和数值解, $ {{\boldsymbol{\varepsilon}} _{\text{I}}} $和$ {{\boldsymbol{\sigma}} _{\text{I}}} $表示应变和应力精确解列向量, $ {\boldsymbol{\varepsilon}} _{\text{I}}^h $和$ {\boldsymbol{\sigma}} _{\text{I}}^h $表示应变和应力数值解列向量.
3.1 二维热传导问题
考虑二维方形区域$ \varOmega = [0,1] \otimes [0,1] $下的热传导问题, 该问题在中心处存在局部集中热源
$$\begin{split} &s({\boldsymbol{x}}) = - 4 - \left\{ 2({\alpha ^2} + {\beta ^2}) - 4{[{\alpha ^2}(x - 0.5)]^2} - 4{[{\beta ^2}(y - 0.5)]^2} \right\}\cdot\\ &\qquad {{\mathrm{e}}^{ - {{[\alpha (x - 0.5)]}^2} - {{[\beta (y - 0.5)]}^2}}}\end{split} $$ (31) 该问题的温度场真实解为$ T({\boldsymbol{x}}) = x(x - 1) + y(y - 1) + {{\mathrm{e}}^{ - {{[\alpha (x - 0.5)]}^2} - {{[\beta (y - 0.5)]}^2}}} $, 其中参数$ \alpha $和$ \beta $的取值对温度场分布影响很大. 图7 给出了不同参数取值下该问题的温度场分布情况, 从图中可以看出, 由于受局部热源影响, 参数值增大使得问题域中心处温度解$T({x})$及其梯度解${T_{,x}}({x})$分布范围越集中, 而越靠近边界的温度场分布越平缓.
![]()
图 7 不同参数下二维热传导问题温度场真实解$T({x})$和梯度解${T_{,x}}({x})$Figure 7. The exact solution $T({x})$ and gradient solution ${T_{,x}}({x})$ for the 2D heat conduction problem with different parameters为了凸显本文所提方法在局部精细化处理问题方面的优越性, 选择参数$ \alpha $和$ \beta $为20时的情况进行配点法分析, 并根据精确解对方形求解域四周施加强制边界条件. 首先, 图8给出了稀疏均布离散情况下IGC方法的计算结果, 可以发现, 采用IGC方法所得数值结果与图7真实解相比相差甚远, 不能很好求解该问题. 为了得到更为精确的计算结果, 在图8(a)均布离散基础上采用所提方法对误差较大区域进行靶向精细化处理, 如图9(a)所示, 并采用本文所提出的IGAMC方法对该问题进行求解分析. 图9(b)和图9(c)为相应局部细化离散形式的温度场和其梯度的数值计算结果, 从图中可以看出, 在局部加密情况下, IGAMC方法可以很好地逼近解析解. 图10给出了二次细化集中热源部位的数值结果, 可见采用逐层局部细化高梯度分布区域可有效控制问题离散模型的总自由度数量, 从而在保证计算精度的前提下达到降低计算成本的目的. 图11分别绘出了采用IGC和IGAMC方法的${L_2}$误差和${H_1}$误差收敛对比结果, 从图中可以看出, 相较于IGC方法而言, 所提IGAMC方法利用较少离散点便可得到较高的计算精度.
![]()
图 8 二维热传导问题的均布离散模型及IGC方法计算结果Figure 8. The results of the uniform discretization for the 2D heat conduction problem using IGC method![]()
图 9 第一次模型局部细化下IGAMC的计算结果及误差Figure 9. The results of the first level refinement for the 2D heat conduction problem using IGAMC method![]()
图 10 二维热传导问题第二次模型局部细化下IGAMC的计算结果Figure 10. The results of the second level refinement for the 2D heat conduction problem using IGAMC method![]()
图 11 二维热传导问题的误差收敛结果对比Figure 11. Convergence comparison of L2 and H1 error norms for the 2D heat conduction problem3.2 厚壁圆筒受压问题
考虑如图12所示的厚壁圆筒受压问题, 该圆筒的内径${r_i} = 2{\text{ cm}}$, 外径${r_o} = 5{\text{ cm}}$, 材料的弹性模量为$ E = 2.0 \times {10^7}{\text{ kPa}} $, 泊松比为$\nu {\text{ = }}0.3$, 其内外壁分别承受均布压力${P_1} = 1000{\text{ kPa}}$和${P_2} = 500{\text{ kPa}}$. 由于问题的对称性, 采用1/4模型进行数值分析, 该问题的应力解析解具有如下形式[38]
$$ \left.\begin{split} & {\sigma _{rr}}({\boldsymbol{x}}) = \frac{1}{{{r^2}({r_o}^2 - {r_i}^2)}}[({r_i}^2 - {r^2}){r_o}^2{P_2} - ({r_o}^2 - {r^2}){r_i}^2{P_1}] \\ & {\sigma _{\theta \theta }}({\boldsymbol{x}}) = \frac{1}{{{r^2}({r_o}^2 - {r_i}^2)}}[ - ({r_i}^2 + {r^2}){r_o}^2{P_2} + ({r_o}^2 + {r^2}){r_i}^2{P_1}] \end{split}\right\} $$ (32) 由于该问题内环压力较大, 采用图13所示的等几何均布离散模型求解该问题时, 在靠近内环区域中传统IGC法所获得的径向应力和环向应力存在比较明显的误差. 如若提高内环处应力计算精度, 采用传统IGC法需要更加密集的离散点, 而远离内环的非关键区域也同样得到了细化, 增加了计算成本. 本文所提IGAMC法继承了无网格法灵活建模的优点, 通过对应力集中区域进行靶向精细化离散分析, 可有效降低计算成本. 图14 给出了该问题逐层细化离散模型和相应应力误差结果, 从图中可以看出, IGAMC法不仅保持了等几何精确几何造型的属性, 还可实现全域差异化离散分析; 与传统IGC法相比, 有效提高了该问题内环处的应力计算精度. 同时, 图15给出了采用IGC和IGAMC方法的能量误差收敛对比结果, 进一步验证了本文所提方法在求解应力集中类型问题时具有显著的计算精度.
![]()
图 13 厚壁圆筒受压问题的IGC法均布离散求解结果Figure 13. The results of the uniform discretization for the 2D hollow cylinder pressurized problem using IGC method![]()
图 14 厚壁圆筒受压问题的IGAMC法局部细化离散求解结果Figure 14. The results of the local refinement for the 2D hollow cylinder pressurized problem using IGAMC method![]()
图 15 厚壁圆筒受压问题的能量误差收敛结果对比Figure 15. Convergence comparison of the energy error norms for the 2D hollow cylinder pressurized problem3.3 三维立方体势问题
考虑求解域为$\varOmega = [0,1] \otimes [0,1] \otimes [0,1]$的三维立方体势问题, 假定该问题的真实解具有以下形式
$$ \begin{split} & u({\boldsymbol{x}}) = \frac{{{x^3} + {y^3} + {z^3}}}{3} + {\mathrm{exp}}\Big\{ - {[\alpha (x - 0.5)]^2} - {[\beta (y - 0.5)]^2} - \\ &\qquad [\gamma (z - 0.5)]\Big\}^2,\quad \alpha = \beta = \gamma = 20 \end{split} $$ (33) 在求解过程中, 采用解析解式(33)对立方体边界面施加强制边界条件. 该问题的解析解如图16所示, 从图中可以发现, 在立方体中心位置位移解的峰值较大. 图17给出了疏密两种均布离散形式及采用传统IGC法的计算结果, 从图中可以看出, 在稀疏离散情况下, 立方体中心区域的位移解计算结果与真实解相差很大; 当离散点加密近3倍后, 中心区域位移解仍然存在一定的误差, 需要进一步细化才能得到较为精确的计算结果.
为了有效降低处理该类问题的计算成本, 在图17(a)所示的离散模型基础上, 将其中心区域分两次进行逐层细化离散, 并采用本文所提IGAMC方法进行计算分析. 图18为第一次局部细化的计算结果, 对比图17(f)可以看出, IGAMC法采用3240个离散点的计算结果要优于IGC法采用4913个离散点的计算结果. 图19为中心求解域进行二次细化后的计算结果, 显然二次细化后基本接近于真实解. 与此同时, 图20给出了两种方法分别采用均布离散和局部细化算法时中心点位移解峰值的收敛情况. 结果进一步表明, 通过对局部特殊部位进行细化处理, 本文所提方法具有更高的收敛精度. 需要说明的是, 该问题峰值出现在中心节点处, 而均布离散细化模型的节点分布不一定恰好坐落在中心位置, 这也就意味着该中心位置位移只能通过周围节点的位移值光滑近似得到. 此时, 该点峰值会产生低于精确解现象, 但会随着进一步节点细化得到缓解.
![]()
图 17 三维立方体势问题均布离散形式下IGC法的计算结果Figure 17. The results of the uniform discretization for the 3D cube potential problem using IGC method![]()
图 18 三维立方体势问题第一次局部细化下IGAMC法的计算结果Figure 18. The results of the first level local refinement for the 3D cube potential problem using IGAMC method![]()
图 19 三维立方体势问题第二次局部细化下IGAMC法的计算结果Figure 19. The results of the second level local refinement for the 3D cube potential problem using IGAMC method![]()
图 20 三维立方体势问题的中心点峰值位移收敛结果对比Figure 20. Convergence comparison of the peak displacement for the 3D cube potential problem4. 结论
本文以无网格法再生梯度的基本理论为基础, 构建了由等几何基函数再生点定义的混合梯度基向量, 提出了等几何B样条基函数梯度的无网格重构形式, 从而发展了一种鲁棒的等几何无网格配点法. 与传统等几何B样条基函数梯度的递推方法相比, 所提方法数值实现更为简捷. 同时, 与等几何无网格形函数的直接梯度相比, 所提无网格再生梯度避免了直接梯度的复杂耗时问题, 且在不同节点序列精确满足不同阶次的多项式再生条件, 保证了数值计算分析的收敛性. 更重要的是, 本文所提理论框架兼具了等几何分析精确建模和无网格法易于实现计算模型局部细化的优势, 解决了传统等几何分析精细化处理模型过程中的额外引入自由度问题, 明显提高了计算效率. 文中通过若干典型算例系统验证了所提配点法在求解局部高梯度问题及应力集中弹性力学问题方面具有显著收敛精度的优势. 未来本文所提方法将在求解接触问题、裂纹扩展和大变形等各类复杂工程问题方面具有广泛的应用前景.
附录A
由式(19)可知, 二维二次混合基向量为
$$ \begin{split} & {\boldsymbol{p}}({\boldsymbol{\xi}} _{I}^{[\centerdot ]})=\{1,\xi _{a}^{[1]},\eta _{b}^{[1]},\xi _{a}^{[1]}\eta _{b}^{[1]},{{(\xi _{a}^{[2]})}^{2}},{{(\eta _{b}^{[2]})}^{2}}, \\ & \qquad \xi _{a}^{[1]}{{(\eta _{b}^{[2]})}^{2}},{{(\xi _{a}^{[2]})}^{2}}\eta _{b}^{[1]},{{(\xi _{a}^{[2]})}^{2}}{{(\eta _{b}^{[2]})}^{2}}{{\}}^{{\mathrm{T}}}}\end{split} $$ (A1) 二维一阶梯度基向量由式(26)可表示为
$$ \left. \begin{split} & {\boldsymbol{p}}_{,\xi }^{{}}({\boldsymbol{\xi}} _{{}}^{[\centerdot ]})={{\{0,1,0,\eta ,2\xi ,0,{{\eta }^{2}},2\xi \eta ,2\xi {{\eta }^{2}}\}}^{{\mathrm{T}}}} \\ & {\boldsymbol{p}}_{,\eta }^{{}}({\boldsymbol{\xi}} _{{}}^{[\centerdot ]})={{\{0,0,1,\xi ,0,2\eta ,2\xi \eta ,{{\xi }^{2}},2{{\xi }^{2}}\eta \}}^{{\mathrm{T}}}}\end{split}\right\}$$ (A2) 相应的二阶梯度基向量由式(27)可表示为
$$ \left. \begin{split} & {\boldsymbol{p}}_{,\xi \xi }^{{}}({{{\boldsymbol{\xi}} }^{[\centerdot ]}})={{\{0,0,0,0,2,0,0,2\eta ,2{{\eta }^{2}}\}}^{{\mathrm{T}}}} \\ & {\boldsymbol{p}}_{,\eta \eta }^{{}}({\boldsymbol{\xi}} _{{}}^{[\centerdot ]})={{\{0,0,0,0,0,2,2\xi ,0,2{{\xi }^{2}}\}}^{{\mathrm{T}}}} \\ & {\boldsymbol{p}}_{,\xi \eta }^{{}}({\boldsymbol{\xi}} _{{}}^{[\centerdot ]})={{\{0,0,0,1,0,0,2\eta ,2\xi ,4\xi \eta \}}^{{\mathrm{T}}}} \end{split} \right\} $$ (A3) 考虑以$ {\boldsymbol{\xi}} \in [{\xi _{a + 1}},{\xi _{a + 2}}] \otimes [{\eta _{b + 1}},{\eta _{b + 2}}] $单元区间为研究区域, 并分别选取$\xi $和$\eta $方向$s = 1.5$倍单元尺寸为形函数影响域, 则该研究区域内共存在9条形函数. 为了简化理论证明, 本文以形函数$ \varPsi _A^{[2]}({\boldsymbol{\xi}} ) $为例, 研究其在$\xi $方向的一阶再生梯度$ \varPsi _{A,\xi }^{[2]}({\boldsymbol{\xi}} ) $及二阶再生梯度$ \varPsi _{A,\xi \xi }^{[2]}({\boldsymbol{\xi}} ) $具体表达形式, $\eta $方向上的各阶梯度简化形式证明与该推导过程类似.
形函数中矩量矩阵式(23)可写为
$$ \begin{split} &{\boldsymbol{M}}({\boldsymbol{\xi}} ) = {\displaystyle \sum _{A = a-1,b-1}^{a + 1,b + 1}{\boldsymbol{p}}({{\boldsymbol{\xi}} }_{A}^{[\cdot]}){{\boldsymbol{p}}}^{{\mathrm{T}}}({{\boldsymbol{\xi}} }_{A}^{[\cdot]}){\varphi }_{s}({{\boldsymbol{\xi}} }_{A}^{[1]}-{\boldsymbol{\xi}} )}= {\boldsymbol{p}}({\xi }_{a-1}^{[\cdot]},{\eta }_{b-1}^{[\cdot]}){{\boldsymbol{p}}}^{{\mathrm{T}}}({\xi }_{a-1}^{[\cdot]},{\eta }_{b-1}^{[\cdot]}){\varphi }_{a-1}^{}{\varphi }_{b-1}^{}+\\ &\qquad {\boldsymbol{p}}({\xi }_{a-1}^{[\cdot]},{\eta }_{b}^{[\cdot]}){{\boldsymbol{p}}}^{{\mathrm{T}}}({\xi }_{a-1}^{[\cdot]},{\eta }_{b}^{[\cdot]}){\varphi }_{a-1}^{}{\varphi }_{b}^{}+ {\boldsymbol{p}}({\xi }_{a-1}^{[\cdot]},{\eta }_{b + 1}^{[\cdot]}){{\boldsymbol{p}}}^{{\mathrm{T}}}({\xi }_{a-1}^{[\cdot]},{\eta }_{b + 1}^{[\cdot]}){\varphi }_{a-1}^{}{\varphi }_{b + 1}^{}+\\ &\qquad {\boldsymbol{p}}({\xi }_{a}^{[\cdot]},{\eta }_{b-1}^{[\cdot]}){{\boldsymbol{p}}}^{{\mathrm{T}}}({\xi }_{a}^{[\cdot]},{\eta }_{b-1}^{[\cdot]}){\varphi }_{a}^{}{\varphi }_{b-1}^{}+ {\boldsymbol{p}}({\xi }_{a}^{[\cdot]},{\eta }_{b}^{[\cdot]}){{\boldsymbol{p}}}^{{\mathrm{T}}}({\xi }_{a}^{[\cdot]},{\eta }_{b}^{[\cdot]}){\varphi }_{a}^{}{\varphi }_{b}^{}+\\ &\qquad {\boldsymbol{p}}({\xi }_{a}^{[\cdot]},{\eta }_{b + 1}^{[\cdot]}){{\boldsymbol{p}}}^{{\mathrm{T}}}({\xi }_{a}^{[\cdot]},{\eta }_{b + 1}^{[\cdot]}){\varphi }_{a}^{}{\varphi }_{b + 1}^{}+ {\boldsymbol{p}}({\xi }_{a + 1}^{[\cdot]},{\eta }_{b-1}^{[\cdot]}){{\boldsymbol{p}}}^{{\mathrm{T}}}({\xi }_{a + 1}^{[\cdot]},{\eta }_{b-1}^{[\cdot]}){\varphi }_{a + 1}^{}{\varphi }_{b-1}^{}+\\ &\qquad {\boldsymbol{p}}({\xi }_{a + 1}^{[\cdot]},{\eta }_{b}^{[\cdot]}){{\boldsymbol{p}}}^{{\mathrm{T}}}({\xi }_{a + 1}^{[\cdot]},{\eta }_{b}^{[\cdot]}){\varphi }_{a + 1}^{}{\varphi }_{b}^{}+ {\boldsymbol{p}}({\xi }_{a + 1}^{[\cdot]},{\eta }_{b + 1}^{[\cdot]}){{\boldsymbol{p}}}^{{\mathrm{T}}}({\xi }_{a + 1}^{[\cdot]},{\eta }_{b + 1}^{[\cdot]}){\varphi }_{a + 1}^{}{\varphi }_{b + 1}^{}\end{split} $$ (A4) 式中, 矩量矩阵$ {\boldsymbol{M}}({\boldsymbol{\xi}} ) $为$9 \times 9$对称矩阵, 其各项元素可参见文献[33].
基于无网格再生梯度式(25), 无网格再生一阶梯度$ \varPsi _{A,\xi }^{[2]}({\boldsymbol{\xi}} ) $和二阶梯度$ \varPsi _{A,\xi \xi }^{[2]}({\boldsymbol{\xi}} ) $分别为
$$ \begin{split} &{\varPsi }_{A,\xi }^{[2]}({\boldsymbol{\xi}} ) = {{\boldsymbol{p}}}_{,\xi }^{{\mathrm{T}}}({{\boldsymbol{\xi}} }_{}^{[\cdot]}){{\boldsymbol{M}}}^{-1}({\boldsymbol{\xi}} ){\boldsymbol{p}}({\xi }_{a}^{[\cdot]},{\eta }_{b}^{[\cdot]}){\varphi }_{a}{\varphi }_{b}=\\ &\qquad -\frac{2\left({\xi }_{a}{\xi }_{a + 1}-{\xi }_{a}\xi -{\xi }_{a + 1}\xi -{\xi }_{a + 2}{\xi }_{a + 3} + {\xi }_{a + 2}\xi + {\xi }_{a + 3}\xi\right)}{\left[{\xi }_{a}{\xi }_{a + 1}^{2}-{\xi }_{a}{\xi }_{a + 1}{\xi }_{a + 2}-{\xi }_{a}{\xi }_{a + 1}{\xi }_{a + 3} + {\xi }_{a}{\xi }_{a + 2}{\xi }_{a + 3}-{\xi }_{a + 1}^{2}{\xi }_{a + 2} + {\xi }_{a + 1}{\xi }_{a + 2}^{2} + {\xi }_{a + 1}{\xi }_{a + 2}{\xi }_{a + 3}-{\xi }_{a + 2}^{2}{\xi }_{a + 3}\right]}\cdot\\ &\qquad \frac{\left[{\eta }_{b}{\eta }^{2} + {\eta }_{b + 1}{\eta }^{2}-{\eta }_{b + 2}{\eta }^{2}-{\eta }_{b + 3}{\eta }^{2} -2{\eta }_{b}{\eta }_{b + 1}\eta + 2{\eta }_{b + 2}{\eta }_{b + 3}\eta + {\eta }_{b}{\eta }_{b + 1}{\eta }_{b + 2} + {\eta }_{b}{\eta }_{b + 1}{\eta }_{b + 3} -{\eta }_{b}{\eta }_{b + 2}{\eta }_{b + 3}-{\eta }_{b + 1}{\eta }_{b + 2}{\eta }_{b + 3}\right]}{\left[{\eta }_{b}{\eta }_{b + 1}^{2}-{\eta }_{b}{\eta }_{b + 1}{\eta }_{b + 2}-{\eta }_{b}{\eta }_{b + 1}{\eta }_{b + 3} + {\eta }_{b}{\eta }_{b + 2}{\eta }_{b + 3}-{\eta }_{b + 1}^{2}{\eta }_{b + 2} + {\eta }_{b + 1}{\eta }_{b + 2}^{2} + {\eta }_{b + 1}{\eta }_{b + 2}{\eta }_{b + 3}-{\eta }_{b + 2}^{2}{\eta }_{b + 3}\right]}=\\ &\qquad \left[\frac{{\xi }_{a + 2} + {\xi }_{a}-2\xi }{({\xi }_{a + 2}-{\xi }_{a})({\xi }_{a + 2}-{\xi }_{a + 1})} + \frac{{\xi }_{a + 3} + {\xi }_{a + 1}-2\xi }{({\xi }_{a + 3}-{\xi }_{a + 1})({\xi }_{a + 2}-{\xi }_{a + 1})}\right]\cdot \left[\frac{(\eta -{\eta }_{b})({\eta }_{b + 2}-\eta )}{({\eta }_{b + 2}-{\eta }_{b})({\eta }_{b + 2}-{\eta }_{b + 1})} + \frac{({\eta }_{b + 3}-\eta )(\eta -{\eta }_{b + 1})}{({\eta }_{b + 3}-{\eta }_{b + 1})({\eta }_{b + 2}-{\eta }_{b + 1})}\right]= {N}_{a,\xi }^{2}(\xi ){N}_{b}^{2}(\eta ) =\\ &\qquad {N}_{A,\xi }^{2}({\boldsymbol{\xi}} )\end{split} $$ (A5) $$ \begin{split} &{\varPsi }_{A,\xi \xi }^{[2]}({\boldsymbol{\xi}} ) = {{\boldsymbol{p}}}_{,\xi \xi }^{{\mathrm{T}}}({{\boldsymbol{\xi}} }_{}^{[\cdot]}){{\boldsymbol{M}}}^{-1}(\xi ){\boldsymbol{p}}({{\boldsymbol{\xi}} }_{a}^{[\cdot]},{\eta }_{b}^{[\cdot]}){\varphi }_{a}{\varphi }_{b}= -\frac{2(-{\xi }_{a}-{\xi }_{a + 1} + {\xi }_{a + 2} + {\xi }_{a + 3})}{\left[{\xi }_{a}{\xi }_{a + 1}^{2}-{\xi }_{a}{\xi }_{a + 1}{\xi }_{a + 2}-{\xi }_{a}{\xi }_{a + 1}{\xi }_{a + 3} + {\xi }_{a}{\xi }_{a + 2}{\xi }_{a + 3}-{\xi }_{a + 1}^{2}{\xi }_{a + 2} + {\xi }_{a + 1}{\xi }_{a + 2}^{2} + {\xi }_{a + 1}{\xi }_{a + 2}{\xi }_{a + 3}-{\xi }_{a + 2}^{2}{\xi }_{a + 3}\right]}\cdot\\ &\qquad \frac{\left[{\eta }_{b}{\eta }^{2} + {\eta }_{b + 1}{\eta }^{2}-{\eta }_{b + 2}{\eta }^{2}-{\eta }_{b + 3}{\eta }^{2} -2{\eta }_{b}{\eta }_{b + 1}\eta + 2{\eta }_{b + 2}{\eta }_{b + 3}\eta + {\eta }_{b}{\eta }_{b + 1}{\eta }_{b + 2} + {\eta }_{b}{\eta }_{b + 1}{\eta }_{b + 3}-{\eta }_{b}{\eta }_{b + 2}{\eta }_{b + 3}-{\eta }_{b + 1}{\eta }_{b + 2}{\eta }_{b + 3}\right]}{\left[{\eta }_{b}{\eta }_{b + 1}^{2}-{\eta }_{b}{\eta }_{b + 1}{\eta }_{b + 2}-{\eta }_{b}{\eta }_{b + 1}{\eta }_{b + 3} + {\eta }_{b}{\eta }_{b + 2}{\eta }_{b + 3}-{\eta }_{b + 1}^{2}{\eta }_{b + 2} + {\eta }_{b + 1}{\eta }_{b + 2}^{2} + {\eta }_{b + 1}{\eta }_{b + 2}{\eta }_{b + 3}-{\eta }_{b + 2}^{2}{\eta }_{b + 3}\right]}= \\ &\qquad -\left[\frac{2}{({\xi }_{a + 2}-{\xi }_{a})({\xi }_{a + 2}-{\xi }_{a + 1})} + \frac{2}{({\xi }_{a + 3}-{\xi }_{a + 1})({\xi }_{a + 2}-{\xi }_{a + 1})}\right]\cdot \left[\frac{(\eta -{\eta }_{b})({\eta }_{b + 2}-\eta )}{({\eta }_{b + 2}-{\eta }_{b})({\eta }_{b + 2}-{\eta }_{b + 1})} + \frac{({\eta }_{b + 3}-\eta )(\eta -{\eta }_{b + 1})}{({\eta }_{b + 3}-{\eta }_{b + 1})({\eta }_{b + 2}-{\eta }_{b + 1})}\right] =\\ &\qquad {N}_{a,\xi \xi }^{2}({\boldsymbol{\xi}} ){N}_{b}^{2}(\eta ) = {N}_{A,\xi \xi }^{2}({\boldsymbol{\xi}} )\end{split} $$ (A6) -
![]()
图 1 1/4圆环曲面的等几何局部细化
Figure 1. Local refinement for a quarter of the annulus in isogeometric analysis
![]()
图 3 模型细化离散下的等几何无网格形函数
Figure 3. The isogeometric meshfree shape functions under the model different level refinements
![]()
图 4 二次混合基函数下等几何无网格形函数直接梯度的再生条件数值验证
Figure 4. Numerical verification of the direct gradient reproducing conditions for 2D isogeometric meshfree shape functions with quadratic mixed basis function
![]()
图 5 二维二次B样条基函数梯度的再生核无网格表述
Figure 5. Reproducing kernel meshfree representations of the gradients for 2D quadratic B-spline basis functions
![]()
图 6 二维二次混合基函数下无网格再生梯度的再生条件数值验证
Figure 6. Numerical verification of the reproducing conditions for 2D meshfree reproducing gradients with quadratic mixed basis function
![]()
图 7 不同参数下二维热传导问题温度场真实解$T({x})$和梯度解${T_{,x}}({x})$
Figure 7. The exact solution $T({x})$ and gradient solution ${T_{,x}}({x})$ for the 2D heat conduction problem with different parameters
![]()
图 8 二维热传导问题的均布离散模型及IGC方法计算结果
Figure 8. The results of the uniform discretization for the 2D heat conduction problem using IGC method
![]()
图 9 第一次模型局部细化下IGAMC的计算结果及误差
Figure 9. The results of the first level refinement for the 2D heat conduction problem using IGAMC method
![]()
图 10 二维热传导问题第二次模型局部细化下IGAMC的计算结果
Figure 10. The results of the second level refinement for the 2D heat conduction problem using IGAMC method
![]()
图 11 二维热传导问题的误差收敛结果对比
Figure 11. Convergence comparison of L2 and H1 error norms for the 2D heat conduction problem
![]()
图 13 厚壁圆筒受压问题的IGC法均布离散求解结果
Figure 13. The results of the uniform discretization for the 2D hollow cylinder pressurized problem using IGC method
![]()
图 14 厚壁圆筒受压问题的IGAMC法局部细化离散求解结果
Figure 14. The results of the local refinement for the 2D hollow cylinder pressurized problem using IGAMC method
![]()
图 15 厚壁圆筒受压问题的能量误差收敛结果对比
Figure 15. Convergence comparison of the energy error norms for the 2D hollow cylinder pressurized problem
![]()
图 17 三维立方体势问题均布离散形式下IGC法的计算结果
Figure 17. The results of the uniform discretization for the 3D cube potential problem using IGC method
![]()
图 18 三维立方体势问题第一次局部细化下IGAMC法的计算结果
Figure 18. The results of the first level local refinement for the 3D cube potential problem using IGAMC method
![]()
图 19 三维立方体势问题第二次局部细化下IGAMC法的计算结果
Figure 19. The results of the second level local refinement for the 3D cube potential problem using IGAMC method
-
[1] Hughes TJR, Cottrell JA, Bazilevs Y. Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement. Computer Methods in Applied Mechanics and Engineering, 2005, 194: 4135-4195 doi: 10.1016/j.cma.2004.10.008
[2] Lin ZW, Wang DD, Wu ZY, et al. A simple non-conforming isogeometric formulation with superior accuracy for free vibration analysis of thin beams and plates. International Journal of Structural Stability and Dynamics, 2023, 23(11): 2350128
[3] 金灵智, 王禹, 郝鹏等. 加筋路径驱动的板壳自适应等几何屈曲分析. 力学学报, 2023, 55(5): 1151-1164 (Jin Lingzhi, Wang Yu, Hao Peng, et al. Adaptive isogeometric buckling analysis of stiffened panels driven by stiffener paths. Chinese Journal of Theoretical and Applied Mechanics, 2023, 55(5): 1151-1164 (in Chinese) Jin Lingzhi, Wang Yu, Hao Peng, et al. Adaptive isogeometric buckling analysis of stiffened panels driven by stiffener paths. Chinese Journal of Theoretical and Applied Mechanics, 2023, 55(5): 1151-1164 (in Chinese)
[4] Agrawal V, Kulachenko A, Scapin N, et al. An efficient isogeometric/finite-difference immersed boundary method for the fluid-structure interactions of slender flexible structures. Computer Methods in Applied Mechanics and Engineering, 2024, 418: 116495 doi: 10.1016/j.cma.2023.116495
[5] Magisano D, Corrado A, Leonetti L, et al. Large deformation Kirchhofflove shell hierarchically enriched with warping: Isogeometric formulation and modeling of alternating stiff/soft layups. Computer Methods in Applied Mechanics and Engineering, 2024, 418: 116556 doi: 10.1016/j.cma.2023.116556
[6] Si ZF, Hirshikesh, Yu TT, et al. An adaptive phase-field simulation for hydrogen embrittlement fracture with multi-patch isogeometric method. Computer Methods in Applied Mechanics and Engineering, 2024, 418: 116539 doi: 10.1016/j.cma.2023.116539
[7] 张雄, 刘岩. 无网格法. 北京: 清华大学出版社, 2004 (Zhang Xiong, Liu Yan. Meshless Methods. Beijing: Tsinghua University Press, 2004 (in Chinese) Zhang Xiong, Liu Yan. Meshless Methods. Beijing: Tsinghua University Press, 2004 (in Chinese)
[8] Li SF, Liu WK. Meshfree Particle Methods. Springer Science & Business Media, 2007
[9] 吴俊超, 吴新瑜, 赵珧冰等. 基于赫林格−赖斯纳变分原理的一致高效无网格本质边界条件施加方法. 力学学报, 2022, 54(12): 3283-3296 (Wu Junchao, Wu Xinyu, Zhao Yaobing, et al. A consistent and efficient method for imposing meshfree essential boundary conditions via hellinger-reissner variational principle. Chinese Journal of Theoretical and Applied Mechanics, 2022, 54(12): 3283-3296 (in Chinese) Wu Junchao, Wu Xinyu, Zhao Yaobing, et al. A consistent and efficient method for imposing meshfree essential boundary conditions via hellinger-reissner variational principle. Chinese Journal of Theoretical and Applied Mechanics, 2022, 54(12): 3283-3296 (in Chinese)
[10] Rastelli P, Vacondio R, Marongiu J. An arbitrarily Lagrangian-Eulerian SPH scheme with implicit iterative particle shifting procedure. Computer Methods in Applied Mechanics and Engineering, 2023, 410: 116015 doi: 10.1016/j.cma.2023.116015
[11] Eirfs A, Ramirez L, Couceiro I, et al. MLS-SPH-ALE: A review of meshless-FV methods and a unifying formulation for particle discretizations. Archives of Computational Methods in Engineering, 2023, 30(8): 4959-4981 doi: 10.1007/s11831-023-09965-2
[12] Jafari H, Malidareh BF, Hosseini VR. Collocation discrete least squares meshless method for solving nonlinear multi-term time fractional differential equations. Engineering Analysis with Boundary Elements, 2024, 158: 107-120 doi: 10.1016/j.enganabound.2023.10.014
[13] Qi DL, Wang DD, Deng LK, et al. Reproducing Kernel mesh-free collocation analysis of structural vibrations. Engineering Computations, 2019, 36(3): 734-764 doi: 10.1108/EC-10-2018-0439
[14] Chen JS, Hu W, Hu HY. Reproducing kernel enhanced local radial basis collocation method. International Journal for Numerical Methods in Engineering, 2010, 75(5): 600-627
[15] 王莉华, 李溢铭, 褚福运. 基于分区径向基函数配点法的大变形分析. 力学学报, 2019, 51(3): 743-753 (Wang Lihua, Li Yiming, Chu Fuyun. Finte subdomain radial basis collocation method for the large deformation analysis. Chinese Journal of Theoretical and Applied Mechanics, 2019, 51(3): 743-753 (in Chinese) Wang Lihua, Li Yiming, Chu Fuyun. Finte subdomain radial basis collocation method for the large deformation analysis. Chinese Journal of Theoretical and Applied Mechanics, 2019, 51(3): 743-753 (in Chinese)
[16] Xue ZY, Wang LH, Ren XD, et al. Weighted radial basis collocation method for large deformation analysis of rubber-like materials. Engineering Analysis with Boundary Elements, 2024, 159: 95-110 doi: 10.1016/j.enganabound.2023.11.016
[17] Wang LH, Hu MH, Zhong Z, et al. Stabilized lagrange interpolation collocation method: A meshfree method incorporating the advantages of finite element method. Computer Methods in Applied Mechanics and Engineering, 2023, 404: 115780 doi: 10.1016/j.cma.2022.115780
[18] 胡明皓, 王莉华. 基于拉格朗日插值的无网格直接配点法和稳定配点法. 力学学报, 2023, 55(7): 1526-1536 (Hu Minghao, Wang Lihua. Direct collocation method and stabilized collocation method based on Lagrange interpolation function. Chinese Journal of Theoretical and Applied Mechanics, 2023, 55(7): 1526-1536 (in Chinese) Hu Minghao, Wang Lihua. Direct collocation method and stabilized collocation method based on Lagrange interpolation function. Chinese Journal of Theoretical and Applied Mechanics, 2023, 55(7): 1526-1536 (in Chinese)
[19] Auricchio F, da Veiga LB, Hughes TJR, et al. Isogeometric collocation methods. Mathematical Models and Methods in Applied Sciences, 2010, 20: 2075-210 doi: 10.1142/S0218202510004878
[20] Schillinger D, Evans JA, Reali A, et al. Isogeometric collocation: Cost comparison with Galerkin methods and extension to adaptive hierarchical NURBS discretizations. Computer Methods in Applied Mechanics and Engineering, 2013, 13(1): 170-232
[21] Fedeli P, Frangi A, Auricchio F, et al. Phase-field modeling for polarization evolution in ferroelectric materials via an isogeometric collocation method. Computer Methods in Applied Mechanics and Engineering, 2019, 351: 789-807 doi: 10.1016/j.cma.2019.04.001
[22] Aronson RM, Evans JA. Divergence-conforming isogeometric collocation methods for the incompressible Navier-Stokes equations. Computer Methods in Applied Mechanics and Engineering, 2023, 410: 115990 doi: 10.1016/j.cma.2023.115990
[23] Sederberg TW, Zheng JM, Bakenov A, et al. T-splines and T-NURCCs. ACM Transactions on Graphics, 2003, 22(3): 477-484 doi: 10.1145/882262.882295
[24] 王悦, 崔雅琦, 於祖庆等. 基于T样条的变网格等几何薄板动力学分析. 力学学报, 2021, 53(8): 2323-2335 (Wang Yue, Cui Yaqi, Yu Zuqing, et al. Dynamic analysis of variable mesh isogeometric thin plate based on T-spline. Chinese Journal of Theoretical and Applied Mechanics, 2021, 53(8): 2323-2335 (in Chinese) Wang Yue, Cui Yaqi, Yu Zuqing, et al. Dynamic analysis of variable mesh isogeometric thin plate based on T-spline. Chinese Journal of Theoretical and Applied Mechanics, 2021, 53(8): 2323-2335 (in Chinese)
[25] Bazilevs Y, Calo VM, Cottrell JA, et al. Isogeometric analysis using T-splines. Computer Methods in Applied Mechanics and Engineering, 2010, 199: 229-263 doi: 10.1016/j.cma.2009.02.036
[26] Scott MA, Li X, Sederberg TW, et al. Local refinement of analysis-suitable T-splines. Computer Methods in Applied Mechanics and Engineering, 2012, 213: 206-222
[27] Chen L, De Borst R. Adaptive refinement of hierarchical T-splines. Computer Methods in Applied Mechanics and Engineering, 2018, 337: 220-245 doi: 10.1016/j.cma.2018.03.032
[28] Casquero H, Wei X, Toshniwal D, et al. Seamless integration of design and Kirchhoff-Love shell analysis using analysis-suitable unstructured T-splines. Computer Methods in Applied Mechanics and Engineering, 2020, 360: 112765 doi: 10.1016/j.cma.2019.112765
[29] Garau EM, Vázquez R. Algorithms for the implementation of adaptive isogeometric methods using hierarchical B-splines. Applied Numerical Mathematics, 2018, 123: 58-87 doi: 10.1016/j.apnum.2017.08.006
[30] Jia Y, Anitescu C, Zhang YJ, et al. An adaptive isogeometric analysis collocation method with a recovery-based error estimator. Computer Methods in Applied Mechanics and Engineering, 2019, 345: 52-74 doi: 10.1016/j.cma.2018.10.039
[31] Fang WH, Yu TT, Yang Y. Adaptive extended isogeometric analysis for steady-state heat transfer in heterogeneous media. Computer Modeling in Engineering and Sciences, 2021, 126(3): 1315-1332 doi: 10.32604/cmes.2021.014575
[32] Xie X, Yang A, Jiang N, et al. Adaptive topology optimization under suitably graded THB-spline refinement and coarsening. International Journal for Numerical Methods in Engineering, 2021, 122(20): 5971-5998 doi: 10.1002/nme.6780
[33] Zhang HJ, Wang DD. Reproducing kernel formulation of B-spline and NURBS basis functions: A meshfree local refinement strategy for isogeometric analysis. Computer Methods in Applied Mechanics and Engineering, 2017, 320: 474-508 doi: 10.1016/j.cma.2017.03.034
[34] Li W, Nguyen-Thanh N, Zhou K. An isogeometric-meshfree collocation approach for two-dimensional elastic fracture problems with contact loading. Engineering Fracture Mechanics, 2020, 223: 106779 doi: 10.1016/j.engfracmech.2019.106779
[35] Huang J, Nguyen-Thanh N, Gao J, et al. Static, free vibration, and buckling analyses of laminated composite plates via an isogeometric meshfree collocation approach. Composite Structures, 2022, 285: 115011 doi: 10.1016/j.compstruct.2021.115011
[36] Chi SW, Chen JS, Hu HY, et al. A gradient reproducing kernel collocation method for boundary value problems. International Journal for Numerical Methods in Engineering, 2013, 93(13): 1381-1402 doi: 10.1002/nme.4432
[37] Wang DD, Zhang HJ. A consistently coupled isogeometric–meshfree method. Computer Methods in Applied Mechanics and Engineering, 2014, 268: 843-870 doi: 10.1016/j.cma.2013.10.014
[38] Gould PL. Introduction to Linear Elasticity. Springer-Verlag, 1999
下载:
计量
- 文章访问数: 141
- HTML全文浏览量: 22
- PDF下载量: 36
