CN109271655B - 一种基于非对称有限元算法的材料尺度效应分析方法 - Google Patents

Abstract

本发明公开了一种基于非对称有限元算法的材料尺度效应分析方法，属于CAE仿真、力学分析计算领域，具体为：1）基于修正的偶应力理论建立材料的本构模型；2）对结构进行有限元网格划分，得到有限元网格模型；3）使用高性能的非对称有限元算法进行数值求解，并提取出最大位移值与最大应力值；4）将有限元网格模型精细化处理后重复步骤3的操作，并判断参考指标是否能够满足收敛性条件或精度要求。本发明解决了常规有限元方法分析材料尺度效应效率低下、精度不足的问题。该方法还具有更好的灵活性与适用性，能够嵌入到现有的商业有限元软件中，从而可以求解更广泛的工程结构。

Description

一种基于非对称有限元算法的材料尺度效应分析方法

技术领域

本发明属于CAE仿真、力学分析计算领域，具体是指一种基于非对称有限元算法的材料尺度效应分析方法。

背景技术

随着材料科学以及生产工艺的不断发展，越来越多的微纳尺度的器件在工业中得到广泛应用。然而当结构特征尺寸处于微纳尺度时，有些材料可能会表现出明显的尺寸效应：结构刚度会随着尺寸的减小而明显增加(Felck NA,Muller GM,Ashby MF,HutchinsonJW.Strain gradient plasticity:theory and experiment.Acta Metall Mater.1994；42(2):475-487.)。在这种情况下，经典的弹性理论不再能够准确地描述材料的力学行为，因此学术界提出了诸多修正的理论模型。这其中修正的偶应力理论被认为是最为简洁的一种，其与经典弹性理论相比只多出一个材料参数(Yang F,Chong A,Lam DCC,TongP.Couple stress based strain gradient theory for elasticity.Int J SolidsStruct.2002；39(10):2731-2743.)。

在工程应用中，纯理论分析的手段在绝大多数情况下是无法解决问题的。在实际中，研究材料尺度效应的可行手段包括实验与有限元数值计算。其中，前者只适合用于分析一些形状规则，载荷形式简单的结构，很难用于解决复杂的问题。因此有限元数值计算通常被认为是工程应用中最为有效的方法。

然而现有的有限元算法在研究材料尺度效应时均存在着不同程度的缺点与不足。例如，一种主流的做法是使用位移的导数作为独立的单元节点自由度来考虑变形曲率的影响(Soh AK,Chen WJ.Finite element formulations of strain gradient theory formicrostructures and the C0-1patch test.Int J Numer Methods Eng.2004；61(3):433-454.)(Ma X,Chen W.24-DOF quadrilateral hybrid stress element for couplestress theory.Comput Mech.2014；53(1):159-172.)(Amanatidou E,Aravas N.Mixedfinite element formulations of strain-gradient elasticity problems.ComputMethods Appl Mech Eng.2002；191(15):1723-1751.)。但是这种做法使得单元的自由度数量急剧增加，而且相应的边界条件也难以确定。同时这种做法也很难嵌入到现有的有限元程序中。

非对称有限元算法是一种具有高精度、高性能的有限元方法应运而生。但是现有的有限元算法或其他数值方法中，往往存在计算精度、计算效率较低的缺陷；此外，现有的有限元算法很难嵌入到现有的商业有限元软件中，限制其使用范围。本发明人在近年来一直从事该领域的研究，并成功提出了一种新的非对称有限元算法框架(Shang Y,OuyangWG.4-node unsymmetric quadrilateral membrane element with drilling DOFsinsensitive to severe mesh-distortion.Int J Numer Methods Eng.2018；113(10):1589-1606.)，有效地克服了原有非对称有限元理论的一些缺点和不足。但是现有的有限元算法中还不能够有效地用于分析研究结构的材料尺度效应，存在分析材料尺度效应效率低下、精度不足的问题。

发明内容

本发明针对现有技术中存在的问题，基于非对称有限元算法框架，公开了一种基于非对称有限元算法的材料尺度效应分析方法，与现有的有限元算法相比，本发明方法具有更高的计算精度与计算效率；此外，本发明方法可以很容易的嵌入到现有的商业有限元软件中，因此具有更好的灵活性与适用性，能够求解更广泛的工程结构。

本发明是这样实现的：

一种基于非对称有限元算法的材料尺度效应分析方法，其特征在于，步骤如下：

步骤一：基于修正的偶应力理论对材料的本构模型进行简化，确定材料的拉梅常数与特征尺度参数；

步骤二：在步骤一的基础上建立结构的力学模型，然后使用软件进行网格划分得到有限元网格模型；所述的结构的力学模型包括步骤一中的材料模型以及几何模型两部分；

步骤三：使用非对称有限元算法对步骤二建立的力学模型进行求解，然后从计算结果中提取出位移与应力各自的最大值，以用于在后续步骤中判断数值计算是否收敛；

步骤四：对有限元网格模型进行精细化处理，然后再次使用非对称有限元算法进行求解，并得到一组新的最大位移值与最大应力值；

步骤五：评估前后两次数值计算的结果是否满足用户预先设定的收敛准则或预先设定的计算精度；

如若不满足，则进一步通过减少网格尺寸、增加网格密度来细化有限元网格模型，并重复步骤四与步骤五，然后再次进行收敛性判断，直至满足设定的收敛准则以及所需的计算精度；

如若满足，则最后一次计算所得的结果即为符合要求的有限元计算最终结果。

进一步，所述的步骤一具体为：

1.1，在修正的偶应力理论中，不仅需要考虑物体位移的一阶导数对应变能的贡献，同时也需要考虑位移的二阶导数对应变能的贡献。修正的偶应力理论中三个独立的位移分量为：

u₁,u₂,u₃ (1)

如上式所述，本申请书中分别使用下标1、2、3表示x、y、z三个坐标轴方向；

应变张量ε的分量可由位移分量对坐标求一阶偏导得到：

ε₁₁＝u_1,1,ε₂₂＝u_2,2,ε₃₃＝u_3,3 (2)

如上式所述，ε_ij,(i,j＝1～3)是应变的9个分量，本申请书中使用下标逗号“，”表示求偏导运算；

曲率张量χ的分量则由位移分量对坐标求二次偏导得到：

式中，χ_ij,(i,j＝1～3)是9个曲率分量；

1.2在修正的偶应力理论中，应力张量σ的分量为：

σ₁₁,σ₂₂,σ₃₃,σ₁₂,σ₁₃,σ₂₃ (6)

其与应变分量之间的本构关系可以表示为：

σ_ij＝λ(ε₁₁+ε₂₂+ε₃₃)δ_ij+2Gε_ij,(i,j＝1～3) (7)

式中λ和G是材料的拉梅(Lame)常数，此处的拉梅常数需要根据材料类别确定，δ_ij是克罗内克(kronecker)符号；

偶应力张量m的分量为：

m₁₁,m₂₂,m₃₃,m₁₂,m₁₃,m₂₃ (8)

其与曲率分量之间的本构关系可以表示为

m_ij＝2Gl²χ_ij,(i,j＝1～3) (9)

式中l是需要根据材料类别确定的特征尺度参数。

进一步，所述的步骤二中的可以使用现有的商业软件例如Abaqus、Hypermesh等对结构进行网格划分，获得用于计算分析的有限元网格模型。

进一步，所述的步骤三的求解过程以虚功原理为理论基础，具体为：

3.1，非对称有限元方法是以虚功原理为理论基础的，其根本特点是采用不同的插值函数分别作为有限元离散的测试函数(Test function)与试探函数(Trial function)，而依此方法构造的单元其刚度矩阵不再具有对称性。首先写出虚功原理的一般表达式：

式中Ω表示问题的控制域，Γ表示控制域的边界。矢量f和R分别表示结构所受到的体积外力与表面外力；δ为变分符号；

3.2，式(10)

是一组满足协调条件的位移场，用于作为单元的测试函数，具体可由单元节点自由度插值而得：

式中q^e是单元节点自由度矢量，

是插值函数矩阵；

表示是由位移

产生的应变矢量：

式(12)具体可根据式(2)与式(3)得到；

表示由位移

产生的曲率矢量：

式(13)具体可根据式(4)与式(5)得到；式(12)和式(13)中，

和

是相应的应变和曲率插值矩阵；

3.3，式(10)中，矢量

表示用于作为单元试探函数的应力场，矢量

表示用于作为单元试探函数的偶应力场，其可以分别表示为

和

是相应的应力与偶应力插值矩阵，具体可以使用解析试函数方法对其进行构造，使得式(14)和式(15)表示的应力场能够事先满足平衡方程，从而可以有效地保证单元具有优良的性能；

3.4，由式(10)出发推导所得最终的有限元刚度矩阵K^e可以表示为：

因为单元采用了两组不同的插值函数分别作为测试函数与试探函数，因此单元刚度矩阵是非对称的。这一点与常规的有限元方法不同。

进一步，所述的步骤四中的精细化处理具体为：对结构的有限元网格模型减少网格尺寸，增加网格密度；将网格的特征尺寸减少为原有尺寸的1/2或2/3，如下式所示：

式中上标n表示有限元网格是第n次划分得到的模型，d取为单元的最大边长，或者是网格划分时所设置的尺寸控制参数，即dⁿ表示为第n次网格划分得到的有限元网格模型的单元最大边长，或者是第n次网格划分时所设置的尺寸控制参数；

在得到新的有限元网格模型后，再次使用非对称有限元算法进行求解，并得到一组新的最大位移值与最大应力值。

进一步，所述的步骤五具体为：根据用户预先设定的收敛准则或计算精度要求判断步骤四以及步骤五的有限元计算结果是否收敛，具体可以采用如下所示的判据：

式中，Δ为用户预先设定的收敛准则或计算精度，Δ的值可由用户自主确定，U是位移幅值绝对值的最大值，σ_max表示主应力绝对值的最大值，上标n表示是第n次有限元计算的结果。即Uⁿ表示为第n次有限元计算所得到的结构位移幅值绝对值的最大值，

表示为第n次有限元计算所得到的结构主应力绝对值的最大值。

本发明与现有技术的有益效果在于：

1)现有技术中一般基于常规的有限元方法进行材料尺度效应的研究，本发明首先利用修正的偶应力理论对材料的本构模型进行适当的简化，然后再使用具有高精度、高性能、列式简洁的非对称有限元算法进行求解；最后，根据是否满足收敛性需求进行精细化处理，并最终得到符合要求的计算结果。

2)本发明方法能够在使用较少计算成本的情况下获得相对较高的计算精度，从而显著地提高了计算效率，此外，本发明方法具有更好的灵活性与适用性，可以很容易嵌入到现有的商业有限元软件中，从而能够求解更广泛的工程结构。

3)此外本发明方法作为一种数值方法，在编写成程序或软件之后，使用者在求解实际工程问题时不需要详细了解算法该工具内部的原理。因此即使是对本研究领域不是很熟悉的工程人员也能够很好的使用本方法用于研究材料尺度效应对工程结构的影响。

附图说明

图1为本发明一种基于非对称有限元算法的材料尺度效应分析方法的原理流程图；

图2为运用本发明方法模拟分析一个微纳尺度梁弯曲问题实例；

图3为运用本发明方法求解实例所得到的不同厚度的微纳尺度梁的弯曲刚度。

具体实施方式

为使本发明的目的、技术方案及效果更加清楚，明确，以下参照附图并举实例对本发明进一步详细说明。应当指出此处所描述的具体实施仅用以解释本发明，并不用于限定本发明。

如图1所示，图1为一种基于非对称有限元算法的材料尺度效应分析方法，主要分为5个部分：

第一部分：基于修正的偶应力理论对材料的本构模型进行简化，确定材料的拉梅常数与特征尺度参数；具体步骤为：

1.1在修正的偶应力理论中，不仅需要考虑物体位移的一阶导数对应变能的贡献，同时也需要考虑位移的二阶导数对应变能的贡献。修正的偶应力理论中三个独立的位移分量为

u₁,u₂,u₃ (1)

如上式所述，本申请书中分别使用下标1、2、3表示x、y、z三个坐标轴方向。应变张量ε的分量可由位移分量对坐标求一阶偏导得到：

ε₁₁＝u_1,1,ε₂₂＝u_2,2,ε₃₃＝u_3,3 (2)

式中，ε_ij,(i,j＝1～3)是应变的9个分量；如上式所述，本申请书中使用下标逗号“，”表示求偏导运算。曲率张量χ的分量则由位移分量对坐标求二次偏导得到：

式中，χ_ij,(i,j＝1～3)是9个曲率分量；

1.2在修正的偶应力理论中，应力张量σ的分量为

σ₁₁,σ₂₂,σ₃₃,σ₁₂,σ₁₃,σ₂₃ (6)

其与应变分量之间的本构关系可以表示为：

σ_ij＝λ(ε₁₁+ε₂₂+ε₃₃)δ_ij+2Gε_ij,(i,j＝1～3) (7)

式中λ和G是需要根据材料类别确定的拉梅(Lame)常数，δ_ij是克罗内克(kronecker)符号；

偶应力张量m的分量为：

m₁₁,m₂₂,m₃₃,m₁₂,m₁₃,m₂₃ (8)

其与曲率分量之间的本构关系可以表示为：

m_ij＝2Gl²χ_ij,(i,j＝1～3) (9)

式中l是材料的特征尺度参数。此处的特征尺度参数与上述的材料的拉梅常数一样，也需要根据材料类别确定。从本构关系式(7)和式(9)可以看出，应力与曲率之间、偶应力与应变之间不存在耦合作用。

第二部分：在第一部分的基础上建立结构的力学模型，然后进行网格划分得到有限元网格模型；具体可以使用现有的商业软件例如Abaqus、Hypermesh等对结构进行网格划分，获得用于计算分析的有限元网格模型。

第三部分：使用非对称有限元算法对第二部分建立的力学模型进行求解，然后从计算结果中提取出位移与应力各自的最大值，以用于在后续步骤中判断数值计算是否收敛；具体步骤为：

3.1，非对称有限元方法是以虚功原理为理论基础的，其根本特点是采用不同的插值函数分别作为有限元离散的测试函数(Test function)与试探函数(Trial function)，而依此方法构造的单元其刚度矩阵不再具有对称性。

3.2，大致过程如下所示，首先写出虚功原理的一般表达式：

式中Ω表示问题的控制域，Γ表示控制域的边界。矢量f和R分别表示结构所受到的体积外力与表面外力；δ为变分符号。

其他变量的解释如下所述：

3.3，式(10)

是一组满足协调条件的位移场，用于作为单元的测试函数，具体可由单元节点自由度插值而得：

式中q^e是单元节点自由度矢量，

是插值函数矩阵。

表示是由位移

产生的应变矢量：

式(12)具体可根据式(2)与式(3)得到。

表示由位移

产生的曲率矢量：

式(13)具体可根据式(4)与式(5)得到。式(12)和式(13)中，

和

是相应的应变和曲率插值矩阵。

3.4，式(10)中，矢量

表示用于作为单元试探函数的应力场，矢量

表示用于作为单元试探函数的偶应力场，其可以分别表示为

和

是相应的应力与偶应力插值矩阵，具体可以使用解析试函数方法对其进行构造，使得式(14)和式(15)表示的应力场能够事先满足平衡方程，从而可以有效地保证单元具有优良的性能。

3.5，由式(10)出发推导所得最终的有限元刚度矩阵K^e可以表示为

因为单元采用了两组不同的插值函数分别作为测试函数与试探函数，因此单元刚度矩阵是非对称的。这一点与常规的有限元方法不同。

第四部分：对结构的有限元网格模型减少网格尺寸，增加网格密度，然后再次使用非对称有限元算法进行求解，并得到一组新的最大位移值与最大应力值；具体步骤为：

对有限元网格模型进行精细化处理，通常方式是将网格的特征尺寸减少为原有尺寸的1/2或2/3，如下式所示

式中上标n表示有限元网格是第n次划分得到的模型，d取为单元的最大边长，或者是网格划分时所设置的尺寸控制参数，即dⁿ表示为第n次网格划分得到的有限元网格模型的单元最大边长，或者是第n次网格划分时所设置的尺寸控制参数。在得到新的有限元网格模型后，再次使用非对称有限元算法进行求解，并得到一组新的最大位移值与最大应力值。

第五部分：评估前后两次数值计算的结果是否满足用户预先设定的收敛准则或预先设定的计算精度。Δ的值可由用户自主确定，通常情况下，取3％～5％可满足工程应用的需求。Δ的值此处取5％。

如若不满足5％，则进一步通过减少网格尺寸、增加网格密度来细化有限元网格模型，并重复步骤四与步骤五，然后再次进行收敛性判断，直至满足收敛准则以及所需的计算精度；如若满足5％，则最后一次计算所得的结果即为符合要求的有限元计算最终结果。具体步骤为：

根据用户预先设定的收敛准则或计算精度要求判断有限元计算结果是否收敛，具体可以采用如下所示的判据：

上式中，U是位移幅值绝对值的最大值，σ_max表示主应力绝对值的最大值，上标n表示是第n次有限元计算的结果。

图2本发明方法模拟分析一个微纳尺度梁弯曲问题实例，是一个微纳尺度梁弯曲问题实例的示意图。材料的杨氏模量E＝1.44GPa,泊松比＝0.38,梁厚度为h，梁跨度L＝10h,梁宽度b＝2h,载荷P＝100N。本实例主要研究不同的材料特征长度、梁厚度h与梁的抗弯刚度三者之间的关系。

图3是采用本发明方法分析所得到的计算结果。从中可以看出，本发明方法在分析材料尺度效应时，具有很高的计算精度。此外本发明方法作为一种数值方法，在编写成程序或软件之后，使用者在求解实际工程问题时不需要详细了解算法该工具内部的原理。因此即使是对本研究领域不是很熟悉的工程人员也能够很好的使用本方法用于研究材料尺度效应对工程结构的影响。

以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。

Claims (4)

1.一种基于非对称有限元算法的材料尺度效应分析方法，其特征在于，步骤如下：

步骤一：基于修正的偶应力理论对材料的本构模型进行简化，确定材料的拉梅常数与特征尺度参数；所述的步骤一具体为：

1.1，修正的偶应力理论中三个独立的位移分量为：

u₁,u₂,u₃ (1)

式中，u表示位移，公式中分别使用下标1、2、3表示x、y、z三个坐标轴方向；

应变张量ε的分量可由位移分量对坐标求一阶偏导得到：

ε₁₁＝u_1,1,ε₂₂＝u_2,2,ε₃₃＝u_3,3 (2)

式中，ε_ij,(i,j＝1～3)是应变的9个分量；公式中使用下标逗号“，”来表示求偏导运算；

曲率张量χ的分量则由位移分量对坐标求二次偏导得到：

式中，χ_ij,(i,j＝1～3)是9个曲率分量；

1.2在修正的偶应力理论中，应力张量σ的分量为：

σ₁₁,σ₂₂,σ₃₃,σ₁₂＝σ₂₁,σ₁₃＝σ₃₁,σ₂₃＝σ₃₂ (6)

其与应变分量之间的本构关系可以表示为：

σ_ij＝λ(ε₁₁+ε₂₂+ε₃₃)δ_ij+2Gε_ij,(i,j＝1～3) (7)

式中λ和G材料的拉梅常数，δ_ij是克罗内克符号；

偶应力张量m的分量为：

m₁₁,m₂₂,m₃₃,m₁₂＝m₂₁,m₁₃＝m₃₁,m₂₃＝m₃₂ (8)

其与曲率分量之间的本构关系可以表示为

m_ij＝2Gl²χ_ij,(i,j＝1～3) (9)

式中l是材料特征尺度参数；

步骤二：在步骤一的基础上建立结构的力学模型，然后使用软件进行网格划分得到有限元网格模型；

步骤三：使用非对称有限元算法对步骤二建立的力学模型进行求解，然后从计算结果中提取出位移与应力各自的最大值；所述的步骤三的求解过程以虚功原理为理论基础，具体为：

3.1，首先写出虚功原理的表达式：

式中Ω表示问题的控制域，Γ表示控制域的边界；矢量f和R分别表示结构所受到的体积外力与表面外力；δ为变分符号；

3.2，式(10)

是一组满足协调条件的位移场，由单元节点自由度插值而得：

式中q^e是单元节点自由度矢量，

是插值函数矩阵；

表示是由位移

产生的应变矢量：

式(12)具体可根据式(2)与式(3)得到；

表示由位移

产生的曲率矢量：

式(13)具体可根据式(4)与式(5)得到；式(12)和式(13)中，

和

是相应的应变和曲率插值矩阵；

3.3，式(10)中，矢量

表示用于作为单元试探函数的应力场，矢量

表示用于作为单元试探函数的偶应力场，其可以分别表示为

和

是相应的应力与偶应力插值矩阵；

3.4，由式(10)出发推导所得最终的有限元刚度矩阵K^e可以表示为：

步骤四：对有限元网格模型进行精细化处理，然后再次使用非对称有限元算法进行求解，并得到一组新的最大位移值与最大应力值；

步骤五：评估前后两次数值计算的结果是否满足用户预先设定的收敛准则或预先设定的计算精度；

如若不满足，则进一步通过减少网格尺寸、增加网格密度来细化有限元网格模型，并重复步骤四与步骤五，然后再次进行收敛性判断，直至满足设定的收敛准则以及所需的计算精度；

如若满足，则最后一次计算所得的结果即为符合要求的有限元计算最终结果。

2.根据权利要求1所述的一种基于非对称有限元算法的材料尺度效应分析方法，其特征在于，所述的步骤二中的软件包括Abaqus、Hypermesh软件对结构进行网格划分，获得用于计算分析的有限元网格模型。

3.根据权利要求1所述的一种基于非对称有限元算法的材料尺度效应分析方法，其特征在于，所述的步骤四中的精细化处理具体为：

对结构的有限元网格模型减少网格尺寸，增加网格密度；将网格的特征尺寸减少为原有尺寸的1/2或2/3，如下式所示：

式中上标n表示有限元网格是第n次划分得到的模型，d取为单元的最大边长，或者是网格划分时所设置的尺寸控制参数，即dⁿ表示为第n次网格划分得到的有限元网格模型的单元最大边长，或者是第n次网格划分时所设置的尺寸控制参数；

在得到新的有限元网格模型后，再次使用非对称有限元算法进行求解，并得到一组新的最大位移值与最大应力值。

4.根据权利要求1所述的一种基于非对称有限元算法的材料尺度效应分析方法，其特征在于，所述的步骤五具体为：

根据用户预先设定的收敛准则或计算精度要求判断步骤四以及步骤五的有限元计算结果是否收敛，具体可以采用如下所示的判据：

式中，Δ为用户预先设定的收敛准则或计算精度，U是结构位移幅值绝对值的最大值，σ_max表示结构主应力绝对值的最大值，上标n表示是第n次有限元计算的结果，即Uⁿ表示为第n次有限元计算所得到的结构位移幅值绝对值的最大值，

表示为第n次有限元计算所得到的结构主应力绝对值的最大值。

Images