计算力学无网格学习报告

一、无网格法的概念

当前数值模拟欧拉方程和 N-S方程的空间离散方法主要有三种:有限差分法、有限体积法和有限元法。有限差分法的主要思想是将微分型控制方程应用到每个网格节点上，对某个节点的导数值用相关网格节点值进行差分求得。再利用各种时间推进格式进行迭代求得流场的稳定解。有限体积法的主要思想是把流体力学的积分型控制方程应用到每一个由网格剖分得到的控制体微元上，从而把积分方程的求解转化为代数方程或常微分方程的求解，再利用各种时间推进格式进行迭代，最后在每个控制体微元上得到流场的稳定解。有限元方法的数值基础是变分原理和加权余量法，其主要思想是在每个网格单元上选择若干合适的点作为解函数的插值节点，方程的近似解由各网格单元的近似函数逼近，而网格中的近似函数又由单元基函数的线性组合得到。以上三种空间离散方法都是在网格的基础上进行的。

无网格方法是一种只需要节点信息而不需将节点连接成网格单元的数值解法，这是无网格算法区别于网格算法之处。因此，无网格算法不能直接沿用网格算法中的有限差分和有限体积法等的离散格式。无网格算法的基本思想是在求解域内布置一系列的离散的点称之为“节点”，然后采用一种形函数(或称为核函数)近似拟合曲线，要求每个小域(称为“点云”)内节点的个数多于建立的方程的数目，在该点云上建立矛盾方程组，通过解矛盾方程组求得变量的梯度，进而可以求得问题的解。

典型无网格法

无网格法大致可分成两类：一类是以Lagrange方法为基础的粒子法（Particle method），如光滑粒子流体动力学（Smoothed particle hydrodynamics，简称SPH）法，和在其基础上发展的运动粒子半隐式（Moving－particle semi－implicit，简称MPS）法等；另一类是以Euler方法为基础的无格子法（Gridless methods），如无格子Euler/N—S算法（Gridless Euler/Navier-Stokes solution algorithm）和无单元Galerkin法（Element free Galerkin，简称EFG）等。

伽辽金型无网格法：

等效积分弱形式（虚位移原理）

特点：计算量大、精度高，稳定性好、需要背景网格进行积分、系数矩阵对称、不易施加本质边界条件处理。

背景网格积分：

加权余量法：

线弹性力学的控制方程

加权余量法不要求余量在各店均为零，而要求余量的加权积分为零 平均意义上满足方程：  等效积分形式

加权余量法的物理意义：选取合适的待定参数强迫余量在某种意义下为零

无网格方法可以方便地利用坐标点计算模拟复杂形状流场计算，但不足之处是在高雷诺数流动时提高数值计算精度较困难。

无网格方法中比较常见的还有径向基函数方法（Radious Basis Function），主要使用某径向基函数（如（MQ）f(r)=r^5）的组合，来逼近原函数。吴忠敏院士在这方面有比较突出的工作。以上方法中，无网格伽辽金法成为目前影响最大，应用最广的无网格计算方法，现有的 LS-dyna，Abaqus，Radioss等商业软件都加入了该方法的计算模块。

二、无网格法应用

悬臂梁在自由端受突加载荷作用后的动力学问题

%--------------------------------------------------------------------------

tic                         % 计时开始

clear;                      % 清屏，将要讨论的是悬臂梁的弹性动力学问题

% DEFINE BOUNDARIES/PARAMETERS   定义边界条件和几何参数

Lx = 100;                   % Lx指的是问题域（矩形域）x向的边长

Ly = 10;                    % Ly指的是矩形域y向的边长

young = 210e9;              % young是材料的弹性模量

nu=0.3;                     % nu是材料的泊松比

midu=10000;                 % midu是材料的弹性模量

%-------------------------------------------------

nx = 20;                    % x方向的划分网格

ny = 6;                     % y方向的划分网格

%-------------------------------------------------

ndivl=10;                   % x方向的积分背景划分网格                

ndivw=5;                    % y方向的积分背景划分网格

dmax=3.8;                   % 影响半径

%--------------------------------------------------

%建立弹性刚度阵(平面应力问题）

Dmat = (young/(1-nu^2))*[1 nu 0;nu 1 0;0 0 (1-nu)/2];

% 建立节点信息

[x,numnod,dm] = mesh1(Lx,Ly,nx,ny,dmax);  % 形成区域内的节点坐标

%作悬臂梁的节点分布图

figure

hold on

plot(x(1,1:(ny+1)),x(2,1:(ny+1)),'k-','linewidth',3);axis equal;                    % 左边界

plot(x(1,(ny+1):(ny+1):numnod),x(2,(ny+1):(ny+1):numnod),'k-','linewidth',2); % 下边界

plot(x(1,numnod:-1:(numnod-ny)),x(2,numnod:-1:(numnod-ny)),'k-','linewidth',2);     % 右边界

plot(x(1,1:(ny+1):(numnod-ny)),x(2,1:(ny+1):(numnod-ny)),'k-','linewidth',2); % 上边界

%plot(x(1,:),x(2,:),'k.');axis off;

plot(x(1,:),x(2,:),'k.');axis equal; axis off;     % 所有的点

hold off

% 建立域内积分单元和角点信息

[xc,conn,numcell,numq] = mesh2(Lx,Ly,ndivl,ndivw);

%建立T1和T2边界上的积分单元和角点信息

[nnu,nnt,numT1,numT2] = mesh3(numq,xc,Lx,Ly);

% 建立高斯点基本信息

quado = 4;

[gauss] = gauss2(quado);

%建立整个求解域内的高斯点信息

numq2 = numcell*quado^2;

gs = zeros(4,numq2);

[gs] = egauss(xc,conn,gauss,numcell);  % 形成所有高斯点的坐标、牙科比、权

%  建立整体质量M矩阵

[M]=Mjuzhen(numnod,gs,x,dm,dmax,midu);

%建立k矩阵

[k]=kjuzhen(numnod,gs,x,dm,dmax,Dmat);

%建立ka矩阵

rfa=200e12;

[ka]=kajuzhen(numnod,nnu,numT1,xc,gauss,x,dm,dmax,rfa);

% 建立K矩阵

K=k+ka;

%建立f矩阵

[f] = fjuzhen( numnod,nnt,numT2,xc,gauss,x,dm,dmax);

%建立fa矩阵

fa = zeros(2*numnod,1);

%建立F矩阵

F=f+fa;

%给出初始条件

tbu=1e-3;

a=1/4;

delta=1/2;

b1=1/(a*tbu^2);

b2=1/(a*tbu);

b3=1/(2*a)-1;

u0=zeros(2*numnod,1);

v0=zeros(2*numnod,1);

a0 = M\\F;

%记入板悬臂端中点轴向位移随时间的变化规律

for i=1:numnod

    xi=x(1,i);

    yi=x(2,i);

if ((abs(xi-100)<100*eps)&(abs(yi-5)<100*eps))

        dian=i;

    end

end

%进入时间循环

uAx=zeros(71,1);

stressxx=zeros(71,1);

uAy=zeros(71,1);

%取B点的坐标，后面好计算该点的应力

gpos=[0;0.05];

v=domain(gpos,x,dm,numnod);

L=length(v);

en=zeros(1,2*L);

[phi,dphix,dphiy] = shape(gpos,dmax,x,v,dm);

Bmat=zeros(3,2*L);

for j=1:L

    Bmat(1:3,(2*j-1):2*j) = [dphix(j) 0;0 dphiy(j);dphiy(j) dphix(j)];

end

for i=1:L

    en(2*i-1) = 2*v(i)-1;

    en(2*i) = 2*v(i);

end

t1=zeros(71,1);

%t2=zeros(101,1);

ind1=1;

%ind2=1;

xiang1=K+b1*M;

for nt=1:1400

    if (nt==1)

        u1=u0;

        v1=v0;

        a1=a0;

    else

        u1=u2;

        v1=v2;

        a1=a2;

    end

    xiang2=b1*u1+b2*v1+b3*a1;

    xiang3=F+M*xiang2;

    u2=xiang1\\xiang3;

    a2=b1*(u2-u1)-b2*v1-b3*a1;

    v2=v1+(1-delta)*tbu*a1+delta*tbu*a2;

 %*********************************************************

    if (nt==ind1*20)  % 隔50个时间步长挑一个点出来

        ind1=ind1+1;

        uAx(ind1)=u2(2*dian-1)*1000;

        t1(ind1)=nt*tbu;

        %计算B点的应力

        stress=zeros(3,1);

        stress = Dmat*Bmat*u2(en);  %stress存放的是B点上的计算应力

        stressxx(ind1)=stress(1,1);

        %

        uAy(ind1)=u2(2*dian)*1000;

    end

%     if (nt==ind1*20)

%         ind1=ind1+1;

%         uAy(ind1)=u2(2*dian)*1000;

%         t2(ind1)=nt*tbu;

%     end

 end

% 计算悬臂端中点的扰度变化解析解i

tx = sqrt(12*midu*Lx^4/(young*Ly^2));

T = 2*3.1415926/(1.875)^2*tx;

Wmax = 2*Ly*Lx^3/(3*young*Ly^3/12);

dt = 0:0.02:1.4;

L=length(dt);

for i=1:L

    W(i) = 1000*(1-cos(2*3.1415/T*dt(i)))*Wmax/2;

end

%plot(dt,-W,'r*');

%作图比较结果

figure          

plot(t1,uAx,'r-');   % a点的水平方向的位移

legend('Displacement of A point in the x direction');

xlabel('t/s','fontweight','bold');

ylabel('UAx ','fontweight','bold');

figure  

hold on 

plot(t1,uAy,'r.');

plot(dt,-W,'K-');

legend('EFG Solution','Exact Solution');

xlabel('t/s','fontweight','bold');

ylabel('UAy ','fontweight','bold');

hold off

figure          

plot(t1,stressxx,'b-');   % B点的水平方向的正应力

legend('Stressxx of A point in the x direction');

xlabel('t/s','fontweight','bold');

ylabel('stressxx','fontweight','bold');

toc

节点分布图