一种基于p型有限元法和围线积分法求解复合型应力强度因
子的计算方法
技术领域
背景技术
[0002] 应力强度因子是判断含裂纹构件断裂及计算裂纹扩展速度的重要参数，在裂纹体分析中占据着重要地位。目前已有众多理论和数值解法求解应力强度因子，如
解析法、混合法、权函数法、线
弹簧法、有限元法、边界元法等。其中有限元法由于其强大的建模能力，可充分利用计算机的计算能力，能够适用于各种复杂的几何情况，能够在各种工程问题中获得较高的
精度，成为求解应力强度因子较为有效的一种方法。
[0003] 有限元解按其结构可分为三类：h型；p型；hp型有限元法，即分别通过加密网格、提高插值多项式阶次或两者同时进行来提高有限元解的精度。相比传统的有限元法，p型有限元方法具有网格划分少、收敛速度快、计算精度高、前处理少等优点。
[0004] 围线积分法基于Betti功互等定理，只需知道远场积分路径上的位移和应力应变，就能计算出混合型裂纹的应力强度因子，导出应力强度因子时具有超收敛性。
[0005] 使用有限元法对断裂力学问题进行研究的大多是传统的h型有限元法或者基于传统有限元的发展(如扩展有限元法等)，因而存在如下问题：
[0006] (1)传统的h型有限元缺乏有效的误差估计模式，计算精度的控制较为依赖研究人员的经验。为获得较好的计算精度，需要多次重新划分网格，来判断有限元解是否收敛；
[0007] (2)传统的h型有限元基于低阶的插值函数，在模拟裂纹尖端的高梯度应力应变场时存在先天不足，为获得较高精度，需要大量的加密网格，前处理及计算成本较大；
[0008] (3)由于裂尖应力场的奇异性，基于传统有限元导出应力强度因子时，由于裂尖应力应变场的误差较大，导出应力强度因子时精度较低；或者要获得较高精度应力强度因子时需要增大前处理和计算成本。
[0009] 本发明基于国家自然科学基金(资助编号：51769011)，提出了将p型有限元法应用于断裂力学领域，并结合围线积分法来求解复合型应力强度因子。
发明内容
[0010] 针对上述
现有技术存在的问题及不足，本发明提供一种基于p型有限元法和围线积分法求解复合型应力强度因子的计算方法。本发明将p型有限元法应用于断裂力学领域，并结合围线积分法，得到一种求解复合型应力强度因子的自适应计算方法。本发明[0011] 一种基于p型有限元法和围线积分法求解复合型应力强度因子的计算方法，按以下步骤进行：
[0012] 步骤1、构建含裂纹构件有限元模型，该步骤包括：
[0013] 根据裂纹构件的实际情况，采集构件的几何参数、材料参数和荷载参数。
[0014] 根据几何参数建立
几何模型，进行几何模型的建立。
[0015] 进行网格划分：依据几何模型建立有限元网格。由于裂纹尖端的奇异性，裂纹尖端的网格应适当加密。裂纹外围应力应变场的梯度较小，网格可以较为稀疏。
[0016] 进行材料设定：根据采集的材料参数建立有限元的材料模型，将对应的材料模型赋予相应的单元。
[0017] 进行
载荷与边界条件的施加：根据采集的荷载参数(位移参数和荷载参数)，在相应边界上施加相应的边界条件；
[0018] 步骤2、将步骤1构建含裂纹构件有限元模型采用p型有限元法计算得到含裂纹构件的位移场、应力场和应变场，
[0019] 具体步骤包括：
[0020] 步骤2.1、求解得到结构位移列阵a
[0021] 根据方程：Ka＝F (1)，其中K＝∑eGTKeG，结构整体
刚度矩阵，Ke为单元刚度矩阵；F＝∑eGTfe，结构结点载荷列阵；G为单位
节点自由度和结构节点自由度的转换矩阵，GT为G的转置，上标T表示矩阵的转置；e表示单个单元；a为结构位移列阵；
[0022] 再根据公式单元刚度矩阵Ke＝∫ΩBTDBdΩ (2)，
[0023] 单元等效节点载荷列阵
[0024] 单元内部节点力
[0025] 外部节点力
[0026] 上述公式(2)至(5)中，Ω表示在单元内部；B＝LNI，L为微分算子，NI为插值函数矩阵或形函数矩阵，I为单元的标号，Γt为单元的外部边界；D为应力矩阵，b为体力，为荷载边界条件，分别由步骤1中的材料属性以及荷载和位移边界条件确定；
[0027] 将公式(2)至(5)以及NI为p型有限元法的插值多项式代入到(1)中，通过求解线性方程组(1)得到结构位移列阵a；
[0028] NI采用p型有限元法的插值多项式，典型的二维p型有限元形函数基于勒让德
正交多项式，以二维四边形单元为例(
附图1所示)，在
坐标系(η,ξ)下，点p1、p2、p3和p4分为单元的四个
顶点，Γ1、Γ2、Γ3和Γ4分为单元的四条边，形函数的构造形式如下所示：
[0029] p≥1，点模式基函数与常规的拉格朗日型基函数时一致的，以四边形四节点单元为例按如下展开：
[0030]
[0031]
[0032]
[0033]
[0034] p≥2，边模式基函数：
[0035] 边Γ1(ξ＝-1,-1≤η≤1)上的形函数可按如下公式给出( 中1为边号，i为插值多项式的阶次)：
[0036]
[0037] 式中：
[0038] 这里，Pn(t)是阶数为n≥0的勒让德多项式
[0039] 类似的，可以如下定义与边Γk(2≤k≤4)相关联的形状函数：
[0040]
[0041]
[0042]
[0043] p≥4时内部模式基函数：
[0044]
[0045] 本发明采用的p型有限元方法，其插值多项式的阶次p＝1,…,∞，提升插值多项式时低阶的刚度矩阵可继续沿用，只需计算高阶部分，避免了刚度矩阵低阶部分的重复计算，具有良好的承袭性，节约了前处理的成本。
[0046] 步骤2.2、根据公式aI＝Ga求解得到单元位移列阵aI；
[0047] 步骤2.3、根据位移场公式u＝NIaI，应变场公式ε＝Lu，应力场公式σ＝Dε＝DBaI，求解得到位移场u，应变场ε，应力场σ。
[0048] 步骤3、采用围线积分法选择包含裂纹尖端的一条围线，根据步骤2得到的该围线上的位移值和应力值导出裂纹处的应力强度因子；
[0049] 本发明采用围线积分法导出应力强度因子。考虑如附图2所示的弹性裂纹场Ω，有局部笛卡尔坐标系(x1,x2,x3)和局部极坐标系(r,θ,z)，z轴与x3重合，用u(x1,x2,x3)和u(r,θ,z)分别表示局部笛卡尔和局部极坐标系下导出域Ωs中裂纹尖端的应变场。其中KⅠ、KⅡ、KⅢ分别为Ⅰ型、Ⅱ型、Ⅲ型应力强度因子(即分别为张开型、滑开型和撕开型裂纹应力强度因子)，位移场u下计算得到围线上
牵引力向量用T(u)表示：
[0050] 在围线上(6)
[0051] 其中公式(6)中， 为位移场u中计算得到的应力张量，nj为围线上外方向法向量分量。
[0052] 在围线积分法中，导出域Ωs退化到域Ωs的一个平面片 上，如附图3所示。应力强度因子KⅠ、KⅡ和KⅢ便可通过围线 和 上的线积分得到：
[0053]
[0054]
[0055]
[0056] 式(7)至(9)中： 和 为导出函数，是裂尖附近弹性解渐近展开的表达式在相应模式下负的特征值，由以下公式给出：
[0057]
[0058]
[0059]
[0060] 其中： G为剪切强度模量，κ为克罗斯夫常数，
平面应力状态下，κ＝(3-ν)/(1+ν)，平面应变状态下，κ＝3-4ν，ν为泊松比；
[0061] 和 分别为虚位移场 和 下计算得到的牵引力向量，同公式(6)；p3和p4为裂纹面上预先给定的牵引力荷载，由步骤1结构的模型参数确定；这里u为步骤2中计算得到的位移场。
[0062] 选择一条包含裂纹尖端的围线，由上述公式导出应力强度因子KⅠ、KⅡ和KⅢ。围线Γ2的选择是很灵活的，不必接近裂纹前端。由于裂纹前端最内层单元中数值解的精度较低，此时导出的应力强度因子误差较大，因而选择围线Γ2时，一般包含裂尖附近的最内层单元。
[0063] 步骤4、判断步骤3所得的应力强度因子是否满足精度要求，如未满足精度要求，提高插值多项式的阶次，返回步骤2：
[0064] 对标准模型，将数值计算结构与理论解或实验结果进行对比，判断数值计算结果是否满足精度要求。
[0065] 对非标准模型，没有理论解或者实验结果进行对比的情况下，依次提高插值多项式阶次，再由围线积分法导出应力强度因子。观察后得到的结果是否与前面计算的结果偏差是否在可接受范围内。
[0066] 精度不满足要求，或所得结果直接偏差不在可接受范围内。再次返回步骤二，重新计算位移场和应力应变场，导出应力强度因子，直至数值计算结果满足精度要求[0067] 所述步骤3中的围线为一条远离裂尖且包含裂尖的围线，且由于裂尖外层单元受奇异性的影响，也应当包含裂纹外层的一层单元。
[0068] 上述同一个标号代表同一个参数值，上标T表示相应该矩阵的转置。
[0069] 本发明的有益效果是：
[0070] (1)本方法能减少计算成本，提高计算效率、收敛速率和精度。
[0071] (2)本方法计算方法能直接应用在实际工程问题中，对含裂纹构件断裂进行计算，且计算精度高。
附图说明
[0072] 图1是本发明标准四边形母单元的示意图；
[0073] 图2是本发明裂纹尖端邻域及全局坐标、局部坐标系的示意图；
[0074] 图3为裂纹尖端退化后的导出域及局部笛卡尔坐标系、极坐标系的示意图；
[0075] 图4是本发明
实施例1边缘裂纹构件的结构示意图；
[0076] 图5为实施例1边缘裂纹网格示意图；
[0077] 图6为实施例1采用p型有限元法计算有限元模型得到的位移场
云图；
[0078] 图7为实施例1采用p型有限元法计算有限元模型得到的应力场云图；
具体实施方式
[0079] 下面结合附图和具体实施方式，对本发明作进一步说明。
[0080] 实施例1
[0081] 如图4所示一边缘裂纹板，其中a＝5为裂缝宽度，板宽b＝10，板长h＝20，板顶端受竖向单位应力σ＝1.0作用，底部左端固定，底部受竖直固定约束。板的材料参数：
杨氏模量E＝107Pa，泊松比ν＝0.3，结构处于平面应变状态。
[0082] 步骤1、构建含裂纹构件有限元模型：
[0083] 首先根据结构的几何参数构建几何模型、定义裂纹；
[0084] 根据结构的几何特征划分网格，模型的网格划分如附图5所示，外围采用稀疏的网格划分，裂纹尖端处于板中央，采用了双层网格加密，最内层网格尺寸是其外层单元的0.15倍。
[0085] 裂纹尖端奇异点的存在，应力梯度较大，这里按0.15的几何因子对裂纹尖端进行了两个网格加密，即内层网格尺寸为外层网格尺寸的0.15倍。本发明中，除了裂纹尖端网格尺寸比较密集之外，其他地方网格网格非常稀疏，p型有限元方法提升阶次的过程中，提高了稀疏网格处的计算精度。
[0086] 根据采集的材料参数建立有限元的材料模型，将相应的材料模型赋予对应的单元；
[0087] 施加结构的位移边界条件以及荷载边界条件：在底部施加竖向的固定约束，左端点施加横向固定约束(位移边界)；顶部施加竖直向上的单位应力(σ＝1.0)(荷载边界)。
[0088] 步骤2、将步骤1构建含裂纹构件有限元模型采用p型有限元法计算得到含裂纹构件的位移场、应力场和应变场，设定有限元计算的插值多项式阶次；这里直接计算了p＝1～8时的应力场和位移场。
[0089] 通过本步骤得到的位移场云图如附图6所示，应力场云图如附图7所示(p≥6时，有限元解的
能量范数误差较低(下表中给出)，云图变化很小，仅给出了p＝1～6时的位移云图和应力云图)。
[0090] 步骤3、采用围线积分法选择包含裂纹尖端的一条围线，根据步骤2得到的该围线上的位移值和应力值导出裂纹处的应力强度因子。
[0091] 选择以裂纹尖端为圆心，以r为半径的一个圆作为围线，利用前面提高的应力强度因子导出公式计算应力强度因子。这里依次导出了p＝1～8时应力强度因子的值。
[0092] r的选择没有特殊的要求，由于裂纹尖端奇异性的存在，当围线处于最内层单元中时，导出的应力强度因子振荡较大，因而围线以包含最内层单元为佳，继续增大围线的半径r，应力强度因子的误差变化较小。
[0093] 本模型具有经验解，由经验公式 计算此时应力强度因子的经验解。
[0094] 上式中： 为经验公式。a/b≤0.6时，函数F：
[0095]
[0096] 当a/b＝0.5时，KI＝11.2018。
[0097] 基于p型有限元求解该构件的位移场、应力应变场通过围线积分法导出应力强度因子，下表给出了p＝1,2,…,8时有限元解的能力范数误差和对应的情形下导出的应力强度因子，及相对误差
[0098]
[0099] 基于p型有限元法，随插值多项式阶次p的增加，有限元解的能力范数误差逐渐收敛。基于有限元解，采用围线积分法导出应力强度因子，应力强度因子的解也随着应力场精度提高而提高。
[0100] 本次计算的模型，结构简单，网格划分较少，所需的计算成本较低，因而直接计算了p＝1～8时的所有情形。
[0101] 当计算的模型结构复杂，网格划分较多时，可依次提高插值多少式的阶次，并依次导出计算的应力强度因子，继续计算至应力强度因子的变化较小时，可认为是比较理想的结果。
[0102] 依次提高插值多项式时，低阶插值多项式的部分不用重复计算，只需计算高阶插值多项式的部分，重新组成高阶的刚度矩阵。
[0103] 基于上表可看出，本发明计算应力强度因子时，具有自适应性，前处理少，计算收敛快，采用较低的计算成本，便可获得较高精度的数值解。
[0104] 实例2
[0105] 模型如实例1所示，改变裂纹尺寸a。
[0106] 实例1中，当插值多项式阶次p＝5时，应力强度因子的精度低于1％。
[0107] 在实例2中，在相似网格条件下，选择插值多项式p＝6，确保应力强度因子精度满足要求。
[0108] 建立模型、划分网格、材料设定、施加边界条件后，直接计算p＝6时的位移场、应力场、应变场，采用围线积分法应力强度因子。改变裂纹长度，继续计算当a/b＝0.1,0.2,…,0.6时的应力强度因子，具体数据如下表：
[0109]
[0110] 表中KΙ为使用本方法计算的结果，KⅠ解析为经验公式计算的结果。
[0111] 由上表可看出改变裂纹长度，网格条件略有不同，对比经验解，此时应力强度因子的精度较高，且具有良好的数值
稳定性。
[0112] 以上结合附图对本发明的具体实施方式作了详细说明，但是本发明并不限于上述实施方式，在本领域普通技术人员所具备的知识范围内，还可以在不脱离本发明宗旨的前提下作出各种变化。
