多尺度边界元方法
本帖最后由 hdy 于 2025-6-7 22:52 编辑弹性非均匀材料全局/局部力学分析的直接多尺度边界元方法
[*]南京航空航天大学能源与动力工程学院，南京 210016
[*]南京航空航天大学机械结构力学及控制国家重点实验室，南京 210016
摘要提出了一种用于弹性非均匀材料全局和局部力学分析的直接多尺度边界元方法。该方法比基于有限元的多尺度方法更准确、方便。通过结合均质化方法和边界元方法的基本方程推导了该方法的公式。该解给出了计算全局弹性常数和局部应力场的便捷公式。最后，给出了两个多孔材料的数值算例，以证明所提方法的准确性和效率。结果表明，该方法不需要迭代即可获得微观水平下的位移解。关键词：多尺度方法；边界元方法；微观结构；均质化方法；全局弹性性质
中图分类号：TB332
文献标识码：A
文章编号：1005-1120(2013)02-0145-10引言评估非均匀材料有效力学性能的多尺度力学方法在当今工程中正成为一种重要方法。自Bensoussan、Lions和Papanicolaou完成开创性工作以来，渐近均质化方法已在许多应用中得到发展，如复合材料拓扑优化、材料降解的多尺度计算建模等。到目前为止，上述领域的大多数多尺度模型都是在有限元方法（FEM）的框架内进行的。基于FEM的主要问题是需要离散化整个代表性体积单元（RVE），当材料微观结构不规则时（尤其是多孔材料中随机缺陷分布的情况），会导致建模过程复杂化。边界元法（BEM）作为一种与FEM截然不同的重要方法，为工程师在某些情况下简化分析提供了新途径。BEM的主要优势在于降低问题的维度，这对于计算成本高昂的大规模建模问题极具吸引力。此外，快速多极方法（FMM）的发展显著提高了BEM的效率。在渐近均匀化理论框架下，已开发出基于BEM的多尺度模型用于预测线性弹性纤维复合材料的有效弹性性能。Kazuhiro等人使用小波方法降低基于BEM的均匀化分析计算成本；Okada等人应用边界元均匀化方法获取橡胶改性环氧树脂和颗粒复合材料的弹性性能；Song Young Seok等人采用多尺度BEM预测聚合物基碳纳米管复合材料的有效弹性性能；近期FMM被引入多尺度BEM用于估计多孔材料的有效弹性性能。然而，这些基于渐近均匀化理论开发的多尺度BEM由于需要求解不同尺度下的平衡方程组而显得复杂。为获得微观位移解，必须使用迭代方案，这降低了方法效率。本文提出的多尺度BEM公式无需使用渐近均匀化理论，比先驱公式更简单直接，解决了微观位移求解的迭代问题。1 基本方程1.1 多尺度方法的理论框架采用均匀化方法构建代表材料微观结构的代表性体积单元（RVE）内位移场的表达式。设为RVE内一点的位置，位移场可分解为线性项和波动项：其中是线性系数，是波动位移。基于均匀化概念，定义均匀（有效或全局）应变和应力为局部场的RVE平均：当存在孔洞或裂纹等不连续时，均匀化应变和应力重新定义为沿RVE外边界的积分：式中，为位移，为弹性刚度，为RVE外边界的法向量。准确估计RVE的整体响应与所施加的边界条件类型直接相关。可使用四种边界条件：均匀牵引力、均匀位移、混合边界条件和周期性边界条件。
其中，上标L、R、T、B分别代表RVE边界的左、右、上、下部分；是通过柯西关系定义的牵引力：通过组合公式(1, 3-4)，可以证明公式(1)中出现的线性系数等同于均匀应变。因此，公式(1)被重新表述为：图2. 周期性边界条件下典型RVE示意图图中是一个带有圆形孔洞的正方形区域，表示代表性体积单元（RVE）。区域四周标注有L（左）、R（右）、T（上）、B（下），表示不同的边界部分。坐标轴和从区域左下角指向右上角。区域内有位移、、、分别表示左、右、上、下边界的位移，、、、表示相应边界的法向量方向。图中还标注了均匀应变与位移的关系。 目前，周期性边界条件（PBC）通常更受青睐，因为即使在微观结构非周期的情况下，它们也能为非均匀材料的力学性能提供最合理的估计。公式(7)连接了均匀应变和不同尺度下的位移场。最终，理论框架可以由一组方程表示：(1) 位移分解（公式(7)）(2) 周期性边界条件（公式(4-5)）(3) 均匀量的定义（公式(3)）(4) 平衡方程(5) 位移-应变关系(6) 微观尺度本构关系在本文提出的理论框架中，波动位移是基本未知变量。这些方程的解给出了均匀应变与波动位移之间的关系，进而通过公式(9-10)推导出微应变和应力的表达式。应用均匀方程(3)可获得全局应力和应变，最终这些均匀量之间的关系提供了材料有效弹性常数的计算公式。基于此方法，通过应用不同数值方法可发展多种多尺度力学方法。例如Suquet等学者开发了基于有限元的方法和快速傅里叶变换(FFT)方法。通过将该多尺度方法与加权残值法结合，还可推导出高保真度的广义胞元法。在下一节中，我们将通过理论框架与边界元法(BEM)的结合，推导出一种直接的基于边界元法的多尺度方法公式。1.2 多尺度边界元法的基本方程对于二维问题，边界积分方程可表示为：(11)其中， 和分别是边界上的位移和牵引力， 和是附录中给出的基本解， 是所谓的自由项， 为边界。 在RVE内的固体域，前面提到的符号是的子集。为了推导多尺度边界积分方程，上述方程中的位移由方程(3)替代为： 其中是源点的位置，是源点处的波动位移。为数值求解RVE积分方程，将边界离散为一系列单元。选择常量边界单元离散RVE边界，因此边界波动位移和牵引力保持连续。写出方程(12)对每个离散单元的表述： 进行积分后得到线性代数方程组：其中： 重复此过程直至方程(14)中所有冗余未知变量被消除，得到最终方程组：1.3 均匀化本构方程 均匀化本构关系可表示为与之间的关系。根据方程(3)-(7)，各向同性材料的应力公式可表示为： 写出方程(19)对每个单元的离散形式为： 方程(18)表明边界单元的波动位移是均匀应变的线性函数： 将方程(21)代入方程(20)可得与的关系： 简化方程(22)可得均匀化本构关系：其中是均匀化弹性常数。1.4 微场应力公式推导基于等式(22-23)，等效弹性模量的表达式为：该式仅包含RVE的几何和物理参数，在已知材料微观结构时可用于计算等效弹性常数。微应力估计在材料退化与断裂模拟中具有重要意义。根据边界元理论，RVE内任意点的位移表达式为：其离散形式可表示为：）2 数值结果分析2.1 周期性孔洞材料案例为验证方法的准确性，选取孔洞呈周期性分布的多孔材料进行测试。如图3所示，选取图中虚线框标示的微观结构最小单元作为RVE，正方形RVE边长为50m。图4对比展示了RVE的两种离散化方法：(a)基于边界元多尺度方法的60个常数边界单元；(b)基于有限元多尺度方法的400个线性四边形单元（880个节点）。材料固体部分弹性模量为70GPa，泊松比为0.2。
2.2 带孔平面应力分析图6展示了边界元法(BEM)和有限元法(FEM)获得的微观应力场对比，包含六个子图分别显示不同方向的应力分布：
图6 BEM与FEM获得的微观应力场对比图7 受分布力和作用的带孔平面示意图图7为带孔大平面在分布力GPa和GPa作用下的受力示意图。根据Kirsch公式，孔附近应力集中表达式为：其中 GPa， GPa为远场应力值。BEM多尺度方法在孔洞附近显示出比FEM更高的计算精度。图8 不同方法预测应力的比较
[*](a) 沿方向应力分布
[*](b) 沿方向应力分布
对比显示BEM结果（红线）与分析解（黑线）吻合度最高，FEM方法（蓝线）在孔洞附近偏差较大，验证了BEM在应力集中问题中的优势。图9 样本和网格示意图
[*](a) 随机孔洞RVE样本
[*](b) BEM方法边界网格
[*](c) FEM方法体域网格
BEM仅需边界离散，网格生成效率显著高于需要全域划分的FEM方法，尤其适用于复杂结构。
对100个样本的计算表明：
[*]BEM得到的分布与FEM结果接近
[*]BEM总耗时 min，显著低于FEM的 min
图展示了所有样本的值分布，横轴为样本编号（），纵轴为值（）。边界元法（）与有限元法（）数据波动趋势相似，计算时间分别为秒和秒。结论数值实验证明，基于边界元法的多尺度方法能高效估算均匀弹性常数并预测微观场变量（如微应变、微应力）。其优势在于：
[*]解析特性：边界元法的解析特性便于离散化
[*]适用性：适用于随机微观结构的多相材料模拟及复合材料拓扑优化
[*]可扩展性：当前二维公式可通过Kelvin解扩展至三维（作者将在后续论文中详述）
页:
[1]
