二维弹性板无网格EFG数值模拟程序

资 源 简 介

本项目是一个基于MATLAB环境开发的二维弹性板力学分析系统，专门用于演示和实现无网格伽辽金法（Element-Free Galerkin Method, EFG）。作为无网格法学习的入门级必学程序，该项目详细展示了EFG方法的核心算法流程。主要功能包括：1. 域内节点生成与管理，支持规则或不规则节点分布；2. 采用移动最小二乘法（Moving Least Squares, MLS）构造形函数，实现场变量的高精度逼近；3. 建立背景积分网格（Background Integration Mesh），利用高斯积分方案计算刚度矩阵和力向量，主要针对二维平面应力或平面应变问题；4. 实现了这一类无网格方法中至关重要的边界条件处理算法，通常采用罚函数法（Penalty Method）或拉格朗日乘子法来强制施加本质边界条件，解决了MLS形函数不具备Kronecker delta性质的问题；5. 求解线性方程组得到节点位移，并进一步通过形函数导数计算全场的应变和应力分布；6. 提供完整的结果后处理模块，能够绘制位移云图、应力云图及变形后的网格形态。代码编写规范，结构模块化，旨在帮助用户深入理解无网格法的理论基础及其在计算力学中的应用。

详 情 说 明

二维弹性板无网格EFG数值模拟程序

项目简介

本项目是一个基于MATLAB环境开发的计算力学仿真系统，专门用于演示和实现无网格伽辽金法（Element-Free Galerkin Method, EFG）。该程序针对二维弹性力学问题（平面应力状态），详细展示了无网格方法的核心算法流程，包括节点生成、移动最小二乘（MLS）形函数构造、背景网格积分以及基于罚函数法的边界条件处理。作为无网格法的入门级学习资源，代码结构清晰，算法逻辑严密，非常适合用于理解EFG方法的理论基础及其数值实现。

功能特性

• 物理模型：模拟二维悬臂梁结构的弹性变形，基于平面应力（Plane Stress）假设。
• 节点离散：支持生成规则分布的场节点（Field Nodes），作为近似函数的基点。
• 支持域管理：采用圆形支持域，通过比例因子动态控制节点的影响范围。
• 无网格近似：利用移动最小二乘法（MLS）构建高精度的形函数场。
• 数值积分：构建背景积分网格（Background Integration Mesh），均布高斯积分点以计算系统刚度矩阵。
• 边界条件处理：实现了罚函数法（Penalty Method），有效解决了无网格法中形函数缺乏Kronecker delta性质导致的本质边界条件施加难题。
• 后处理可视化：包含全场位移重构、应变计算、应力计算（含Von Mises等效应力）及云图数据通过。

系统要求

• MATLAB R2016a 或更高版本
• 无需额外的工具箱（Toolbox），主要依赖MATLAB基础矩阵运算功能。

算法实现与逻辑详解 (基于 main.m)

该程序的核心逻辑 main.m 按照标准的有限元/无网格分析流程编写，具体实现细节如下：

1. 前处理与参数定义

程序首先定义了物理几何参数（板长、板高、弹性模量、泊松比）以及数值计算参数（节点密度、背景网格密度、支持域半径比例、罚函数因子）。在此阶段，程序根据平面应力假设，计算了各向同性线弹性材料的本构矩阵（D矩阵），为后续刚度计算奠定基础。

2. 节点生成与背景网格构建

• 场节点生成：代码通过双重循环生成了规则分布的节点坐标矩阵，这些节点承载位移自由度。
• 支持域定义：计算节点间距，并根据设定的比例因子（alpha）确定MLS近似的支持域半径（dmax）。
• 背景积分网格：不同于有限元的单元网格，程序划分了规则的矩形背景网格仅用于数值积分。代码记录了每个背景单元的几何中心和长宽尺寸，用于后续的高斯积分映射。

3. 刚度矩阵组装

这是程序的核心计算模块，采用了基于背景网格的全局积分方案：

• 稀疏矩阵初始化：根据节点总数和每个节点的自由度（X, Y方向），初始化全局刚度矩阵和力向量。
• 高斯积分循环：程序遍历每一个背景积分单元，并在每个单元内进行 4x4 的高斯点循环。
• MLS形函数计算：在每个高斯积分点处，调用MLS算法计算该点处的形函数值及其对坐标的偏导数（dphi/dx, dphi/dy）。程序会自动识别落入支持域内的邻近节点。
• 应变-位移矩阵(B矩阵)构造：利用形函数导数构建B矩阵，并不再依赖单元的拓扑连接，而是基于支持域内的节点动态构建。
• 刚度累加：计算局部切线刚度（B' * D * B * 权重），并将其利用索引映射累加至全局稀疏矩阵中。

4. 载荷施加

程序模拟了右端悬臂端的剪切载荷：

• 识别位于右边界（x=L）的所有节点。
• 将总剪切力均分或按设定分布施加到这些节点的Y方向自由度上，直接修改全局力向量。

5. 边界条件处理（罚函数法）

由于MLS形函数不具备插值特性，程序采用罚函数法强制施加本质边界条件（左端固定）：

• 边界积分：沿着左边界（x=0）设置专门的积分点（而非仅仅在节点上）。
• 罚刚度矩阵：在每个边界积分点处重新计算MLS形函数，构建罚项刚度矩阵。该矩阵的形式为大数（Penalty Factor）乘以形函数向量的外积。
• 系统修正：将罚刚度矩阵叠加到全局刚度矩阵的对应自由度上，从而在数值上强制约束左端位移为零。

6. 方程求解

使用MATLAB内置的高效稀疏矩阵求解器（反斜杠运算符），求解线性方程组 K * U = F，获得所有场节点的位移向量。

7. 后处理与结果计算

求解完成后，程序进入可视化数据准备阶段：

• 绘图网格生成：生成高密度的规则网格用于绘制平滑云图。
• 物理量重构：遍历绘图网格的每一个点，再次调用MLS算法。利用求得的节点位移，通过形函数插值计算该点的位移。
• 应力应变计算：利用形函数的导数计算该点的应变分量，进而通过本构关系计算应力分量（Sigma_xx, Sigma_yy, Sigma_xy）。
• Von Mises应力：根据平面应力状态下的应力分量，计算Von Mises等效应力，用于评估材料的屈服状态。

使用方法

1. 确保MATLAB当前工作路径包含 main.m 及相关的MLS计算函数（代码中调用的 calculateMLS 和 getGaussPoints 需存在于同一路径下）。
2. 直接运行 main.m 脚本。
3. 程序将由此在命令行输出执行进度（网格初始化、刚度组装、求解过程）。
4. 运行结束后，程序将自动弹出图形窗口，按代码逻辑绘制节点分布、变形后的网格形态及应力云图（代码主要部分已呈现此逻辑）。