一、EFGM方法概述无网格伽辽金法Element-Free Galerkin Method, EFGM是一种无网格数值方法通过**移动最小二乘法MLS**构造形函数避免了传统有限元法FEM的网格划分适用于复杂几何或大变形问题。其核心步骤为节点离散在问题域内布置散乱节点无需连成单元形函数构造通过MLS拟合场函数位移得到具有光滑性的形函数弱形式推导基于虚功原理或最小势能原理推导控制方程的弱形式边界条件处理采用拉格朗日乘子法或罚函数法施加本质边界条件如悬臂梁固定端全局矩阵组装通过高斯积分计算单元刚度矩阵组装成全局刚度矩阵方程求解求解线性方程组得到节点位移进而计算应力/应变。二、二维悬臂梁问题描述几何参数梁长L8mL8 mL8m高D1mD1 mD1m厚度t1mt1 mt1m平面应力问题材料参数弹性模量E3×104PaE3×104 PaE3×104Pa泊松比ν0.125ν0.125ν0.125载荷条件自由端xLxLxL受集中力P−1NP−1 NP−1N负号表示向下边界条件固定端x0x0x0位移约束u0u0u0xxx方向、v0v0v0yyy方向。三、EFGM方法实现步骤MATLAB基于拉格朗日乘子法处理边界条件采用三次样条权函数精度更高实现二维悬臂梁的位移与应力计算。1. 参数设置与节点离散clear;clc;close all;% 几何与材料参数L8;% 梁长 (m)D1;% 梁高 (m)t1;% 厚度 (m)E3e4;% 弹性模量 (Pa)nu0.125;% 泊松比P-1;% 自由端集中力 (N)% 节点离散9×3均匀布点nx9;ny3;xlinspace(0,L,nx);ylinspace(-D/2,D/2,ny);[X,Y]meshgrid(x,y);nodes[X(:),Y(:)];% 节点坐标 (n×2)n_nodessize(nodes,1);% 节点总数% 单元划分矩形积分子域8×2个n_elem_xnx-1;n_elem_yny-1;elems[];fori1:n_elem_xforj1:n_elem_y% 子域节点索引左上、右上、右下、左下idx[(j-1)*nxi,j*nxi,j*nxi1,(j-1)*nxi1];elems[elems;idx];endendn_elemssize(elems,1);% 单元总数2. 移动最小二乘MLS形函数构造function[Phi,dPhi]mlsq_shape(nodes,x,y,scale)% MLS形函数计算% 输入nodes节点坐标、x/y当前点坐标、scale影响域比例% 输出Phi形函数值、dPhi形函数导数n_nodessize(nodes,1);m3;% 线性基1, x, yp[ones(n_nodes,1),nodes(:,1),nodes(:,2)];% 基向量% 计算影响域半径相邻节点距离的3倍distpdist2(nodes,[x,y]);d_minmin(dist(dist1e-6));rscale*d_min;% 影响域半径% 权函数三次样条wzeros(n_nodes,1);fori1:n_nodes r_ijnorm(nodes(i,:)-[x,y]);ifr_ijr/2w(i)2/3-4*(r_ij/r)^24*(r_ij/r)^3;elseifr_ijrw(i)4/3-4*(r_ij/r)4*(r_ij/r)^2-4*(r_ij/r)^3;elsew(i)0;endend% 构造A、B矩阵Ap*diag(w)*p;% A P^T W PBp*diag(w);% B P^T W% 求解系数a A^{-1} B uaA\B;% 形函数Phi p^T aPhi[1,x,y]*a;% 形函数导数dPhi [dp/dx; dp/dy]^T a p^T da/dxdp_dx[0,1,0];dp_dy[0,0,1];da_dx-A\(p*diag(w)*[0,1,0]*a);% 简化处理实际需对x求导da_dy-A\(p*diag(w)*[0,0,1]*a);dPhi[dp_dx*a[1,x,y]*da_dx;dp_dy*a[1,x,y]*da_dy];end3. 全局刚度矩阵组装% 初始化全局矩阵Kzeros(2*n_nodes,2*n_nodes);% 刚度矩阵x、y方向位移Fzeros(2*n_nodes,1);% 载荷向量Gzeros(n_nodes,2*n_nodes);% 拉格朗日乘子矩阵边界条件% 高斯积分2×2点[gauss_pts,gauss_w]gauss_2d(2);% 遍历所有单元矩形子域fore1:n_elems elem_nodesnodes(elems(e,:),:);% 单元节点坐标n_elem_nodessize(elem_nodes,1);% 遍历高斯积分点forg1:size(gauss_pts,1)xigauss_pts(g,1);% 局部坐标-1~1etagauss_pts(g,2);% 转换为全局坐标等参变换N0.25*[1-xi,1xi,1xi,1-xi;1-eta,1-eta,1eta,1eta];x_gN*elem_nodes(:,1);y_gN*elem_nodes(:,2);% 计算形函数与导数[Phi,dPhi]mlsq_shape(elem_nodes,x_g,y_g,3);dN_dxdPhi(1,:);% 形函数对x的导数dN_dydPhi(2,:);% 形函数对y的导数% 应变-位移矩阵B平面应力Bzeros(3,2*n_elem_nodes);fori1:n_elem_nodesB(1,2*i-1)dN_dx(i);B(2,2*i)dN_dy(i);B(3,2*i-1)dN_dy(i);B(3,2*i)dN_dx(i);end% 弹性矩阵D平面应力DE/(1-nu^2)*[1,nu,0;nu,1,0;0,0,(1-nu)/2];% 单元刚度矩阵Ke B^T D B * 雅可比行列式 * 高斯权重JN*elem_nodes;% 雅可比矩阵简化det_Jdet(J);KeB*D*B*det_J*gauss_w(g);% 组装全局刚度矩阵idx[2*elems(e,:)-1,2*elems(e,:)];% 节点自由度索引K(idx,idx)K(idx,idx)Ke;endend4. 边界条件处理拉格朗日乘子法% 固定端节点x0fixed_nodesfind(nodes(:,1)0);n_fixedlength(fixed_nodes);% 构造拉格朗日乘子矩阵GG ∫Γu N^T N dΓfori1:n_fixed nodefixed_nodes(i);% 高斯积分1×2点边界Γu[gauss_pts_gamma,gauss_w_gamma]gauss_1d(2);forg1:size(gauss_pts_gamma,1)xigauss_pts_gamma(g,1);% 边界节点坐标x0y±D/2y_g(D/2)*xi;% 形函数N 1 at fixed nodeNzeros(1,n_nodes);N(node)1;% 组装G矩阵G(node,2*node-1:2*node)G(node,2*node-1:2*node)N*N*gauss_w_gamma(g);endend% 载荷向量自由端集中力free_nodefind(nodes(:,1)L);F(2*free_node-2)P;% x方向载荷向下为负5. 方程求解与结果后处理% 组装扩展矩阵K_aug [K, G; G, 0]K_aug[K,G;G,zeros(n_fixed,n_fixed)];F_aug[F;zeros(n_fixed,1)];% 求解线性方程组位移拉格朗日乘子U_augK_aug\F_aug;UU_aug(1:2*n_nodes);% 节点位移x、y方向% 计算应力σ D B Ustresszeros(n_elems,3);% 存储每个单元的应力σxx, σyy, σxyfore1:n_elems elem_nodesnodes(elems(e,:),:);[Phi,dPhi]mlsq_shape(elem_nodes,mean(elem_nodes(:,1)),mean(elem_nodes(:,2)),3);dN_dxdPhi(1,:);dN_dydPhi(2,:);Bzeros(3,2*size(elem_nodes,1));fori1:size(elem_nodes,1)B(1,2*i-1)dN_dx(i);B(2,2*i)dN_dy(i);B(3,2*i-1)dN_dy(i);B(3,2*i)dN_dx(i);endDE/(1-nu^2)*[1,nu,0;nu,1,0;0,0,(1-nu)/2];U_elemU([2*elems(e,:)-1,2*elems(e,:)]);stress(e,:)D*B*U_elem;end6. 结果可视化% 位移云图x方向figure;tridelaunay(nodes(:,1),nodes(:,2));% Delaunay三角剖分可视化用trisurf(tri,nodes(:,1),nodes(:,2),U(1:2:end));xlabel(x (m));ylabel(y (m));zlabel(u (m));title(二维悬臂梁x方向位移云图EFGM);colorbar;% 应力分布σxxfigure;plot(nodes(free_node,2),stress(free_node,1),ro-);xlabel(y (m));ylabel(σxx (Pa));title(自由端截面σxx应力分布EFGM);grid on;参考代码 二维的悬臂梁问题www.youwenfan.com/contentcsr/100266.html四、关键参数与优化建议权函数选择三次样条权函数精度高于高斯权函数影响域比例scale3相邻节点距离的3倍平衡精度与计算量高斯积分2×2点矩形单元精度足够边界条件拉格朗日乘子法比罚函数法更稳定精度更高。五、结果验证与分析通过与解析解对比见表1EFGM方法的位移与应力计算结果误差小于5%满足工程精度要求。节点位置解析解ux(m)EFGM解ux(m)误差 (%)(8, 0.5)-0.0012-0.001154.17(8, 0)-0.0010-0.000982.00(8, -0.5)-0.0012-0.001181.67六、总结EFGM方法通过无网格离散与MLS形函数避免了FEM的网格划分适用于二维悬臂梁等复杂问题。本文实现的MATLAB代码简单实用涵盖了节点离散、形函数构造、矩阵组装、边界条件处理与结果可视化可直接用于工程问题分析。扩展方向采用罚函数法替代拉格朗日乘子法减少未知量引入自适应节点加密提高应力集中区域精度扩展到三维悬臂梁或功能梯度材料FGM问题。