有限体积法 订阅
有限体积法是计算流体力学中常用的一种数值算法,有限体积法基于的是积分形式的守恒方程而不是微分方程,该积分形式的守恒方程描述的是计算网格定义的每个控制体。有限体积法着重从物理观点来构造离散方程 ,每一个离散方程都是有限大小体积上某种物理量守恒的表示式 , 推导过程物理概念清晰 ,离散方程系数具有一定的物理意义 ,并可保证离散方程具有守恒特性。 [1] 展开全文
有限体积法是计算流体力学中常用的一种数值算法,有限体积法基于的是积分形式的守恒方程而不是微分方程,该积分形式的守恒方程描述的是计算网格定义的每个控制体。有限体积法着重从物理观点来构造离散方程 ,每一个离散方程都是有限大小体积上某种物理量守恒的表示式 , 推导过程物理概念清晰 ,离散方程系数具有一定的物理意义 ,并可保证离散方程具有守恒特性。 [1]
信息
外文名
finite volume method
别    称
有限容积法、控制体积法
归属学科
计算流体力学
种    类
属于加权剩余法中的子区域法
中文名
有限体积法
基本方法
子区域法
有限体积法定义
有限体积法(Finite Volume Method)又称为有限容积法、控制体积法。
收起全文
精华内容
下载资源
问答
  • 有限体积法

    2015-12-19 17:09:24
    有限体积法讲义,主要讲解了有限体积的理论和方法
  • 有限体积法基础

    2018-08-14 13:16:05
    第1章在比较了几种常用的流体流动数值计算方法特点的基础上着重介绍了有限体积法的基本思想和特点; 第2章介绍扩散问题的有限体积解法,从一维稳态扩散问题人手,简要介绍了区域离散方法、离散方程的推导和控制容积...
  • 有限体积法求解雷诺方程有有限元法、有限差分法、有限体积法,本程序就是采用有限体积法对雷诺方程进行求解。
  • 有限体积法MATLAB求解程序
  • 为了研究非稳态导热问题中的常用有限体积方案(FVM1)和新提出的有限体积方案(FVM2)的精度差异问题,采用理论推导和实例验证的方法,建立基于有限体积法原理的热平衡积分方程,利用矩形网格进行离散化,建立了平面非稳态热...
  • 介绍二维对流扩散方程的有限体积法的程序,通过进行离散化网格,最后计算出温度场。主要是c++程序。 运行环境:visual studio
  • 学过弹性力学的人应该都知道什么是有限元,而对学计算流体力学的来说,有限差分和有限体积法也是两种非常重要的方法。三者虽然目前形式各异,但是思想上有很多类似的地方。CFD(Computational Fluid Dynamics)中...

    学过弹性力学的人应该都知道什么是有限元,而对学计算流体力学的来说,有限差分和有限体积法也是两种非常重要的方法。三者虽然目前形式各异,但是思想上有很多类似的地方。CFD(Computational Fluid Dynamics)中主要的三种离散方法就是他们三个。

    而这篇文章主要目的是对三者进行比较,并给出三种方法计算同一个流体一维算例的过程。

    一维算例:

    流体在沿着一条流线传输过程中,关于流体的某个物理量的输运方程为:

            

    这个方程也可以看成对流-扩散方程的特例,其含义就是在控制体内部:

    由源产生的的量(因为全场没有源项,故等号右端为0) =

    由于流动导致的物理量的增加()+由于扩散导致物理量的减少(

    另外流场还应该满足连续条件:

             

    即:流入质量=流出质量

    那么控制方程组就是:

    下面,用三种方法分别来解这个流场的物理量

    1 有限差分法

    主要参考《计算流体力学入门》,John D.Andersion,JR .

    有限差分法的做法很简单:流体的控制方程组一般都是偏微分方程,解偏微分方程组很困难, N-S方程多数都是没有解析解的,但是解线性方程组是很简单的事情。所以,差分主要目的便是将微分近似为线性运算,进而把偏微分方程组转化为线性方程组。大多数用来代替偏微分形式是基于Taylor展开得到的,利用Taylor展开,对x一阶导数差分可以近似如下三种形式:

             

    对扩散项的二阶导数,近似为:

    将控制方程离散成差分方程,得到:

    这里我们把项转化为一个一阶导数和二阶导数,然后分别一阶和二阶格式进行差分。

    通过联立连续方程,可以求解出全场所有节点处流速值,代入对流扩散方程中,简化为:

    然后联立求解全场各个点有关的方程组,便可解得全场节点处数值。

    2 有限体积法

    主要参考《计算流体动力学分析:CFD软件原理与应用》,王福军

    (其实里面也就是流体力学控制方程简介+有限体积法,整个这个名字真让人迷惑QAQ)

    1.1离散

    有限体积法离散的核心和有限元法一样,使用有限个离散点来代替原来整个连续的空间。把计算区域分成不重叠的计算网格,然后确定每个节点位置和节点控制体体积(也就是节点所在的网格单元)。区域几何要素主要有以下几个:

    节点:需要解未知物理量的几何位置,一般在节点上定义所有的标量,下面图中的W、P、E三个点就是节点;

    控制体积:应用控制方程和守恒定律的最小单位,图中灰色区域;

    界面:规定了与各节点相对应的控制体积的分界面位置,图中w和e处;

    网格线:连接节点的曲线簇,图中就是x轴线。

    1.2建立离散方程

    有限体积最特殊的一步就是在控制体上积分控制方程,以便在控制体积的界面上产生离散方程。对一维算例在控制体积上积分有:

    然后控制体积分式等价于:

    其中 可以看成是垂直于 轴方向上断面的面积。

    这里我想详细的讨论一下,在书中作者说这是完全等价的转化,并没有取任何近似,但是我有不同的看法

    积分式里面的被积函数和A没有任何关系,所以我们可以先把做积分:

    然后根据分部积分方法:

    我们可以看出,等号左端正确结果应为:

    所以这一步其实也做了近似。否则应该加上的前提。

    1.3解方程

    我们已经把非线性的控制方程组通过离散,转化成为了线性方程组。

    但是在我们的解里只定义在了节点上,在界面处应该并没有的值,所以还需要某种差分方法,把在界面处的值用节点上的值来近似,这里就是所谓的使用差分格式。

    采用最煎蛋XD的方法……我们用界面两边的节点值取平均来近似,这应该是地球人都想得到的方法,方程继续化简为:

    记为:

    根据连续方程,我们可以求解整个流场,那么系数即为已知数,然后便可求解

    分别代表了相邻点处的物理量通过对流扩散作用对点影响的大小,我们假设流场是已知的,那么我们已经知道全场每个点的大小。这样把全场的节点离散方程联立求解,便可解得每个节点的值。

    展开全文
  • 有限体积法及其网格简介

    万次阅读 多人点赞 2018-07-08 11:26:58
    有限体积法及其网格简介有限体积法使目前CFD领域广泛使用的离散化方法,其特点不仅表现在对控制方程的离散结果上,还表现在所使用的网格上。1.有限体积法的基本思想有限体积法(Finite Volume Method)又称为控制...

    有限体积法及其网格简介

    有限体积法使目前CFD领域广泛使用的离散化方法,其特点不仅表现在对控制方程的离散结果上,还表现在所使用的网格上。

    1.有限体积法的基本思想

    有限体积法(Finite Volume Method)又称为控制体积法(Control Volume Method),其基本思路是:将计算区域划分为网格,并使每一个网格点周围有一个互不重复的控制体积;将待解微分方程(控制方程)对每一个控制体积积分,从而得出一组离散方程。其中的未知数是网格点上的因变量。为了求出控制体体积的积分,必须假定因变量的值在网格点之间的变化规律。从积分区域的选取方法看来,有限体积法属于加权余量法中的子域法,从未知解的近似方法来看,有限体积法属于采用局部近似的离散方法。简言之,子域法加离散,就是有限体积法的基本方法。

    有限体积法得出的离散方程,要求因变量的积分守恒对任意一组控制体积都得到满足,对整个计算区域,自然也得到满足,这就是有限体积法的优点。而有限差分法仅当网格极其细密时,离散方程才满足积分守恒;而有限体积法既使在粗网格情况下,也显示出准确的积分守恒。

    2.有限体积法所使用的网格

    有限体积法的核心体现在区域离散方式上,区域离散化的实质就是用有限个离散点来代替原来的连续空间。有限体积法的区域离散实施过程是:把所计算的区域划分为多个互不重叠的子区域,即计算网格(grid),然后确定每个子区域中的节点位置及该节点所代表的控制体积。区域离散化过程结束后,可以得到以下四种几何要素:

    l  节点(node):需要求解的未知物理量的几何位置。

    l  控制体积(control volume):应用控制方程或守恒定律的最小几何单位。

    l  界面(face):它规定了与各节点相对应的控制体积的分界面位置。

    l  网格线(grid line):联结相邻两节点而形成的曲线簇

    我们把节点看成是控制体积的代表,在离散过程中,将一个控制体积上的物理量定义并存储在该节点处。下图是一维和二维有限体积法计算网格。

    在上面两图中,节点排列有序,即当给出一个节点的编号后,立即可以得出其相邻节点的编号。这种网格称为结构网格(structured grid)。非结构网格(unstructured grid)的节点以一种不规则的方式布置在流程中,虽然生成较复杂,却有着极大的适应性。

    3.求解一维稳态问题的有限体积法

    无论是连续性方程、动量方程还是能量方程,都可以写成如下的通式。一维问题的控制方程:

    称之为一维模型方程,方程中包含对流项、扩散项及源项。方程中的Φ是广义变量,可以为速度、温度或浓度等一些待求的物理量,Г是相应于Φ的广义扩散系数,S是广义源项。变量Φ在端点A和B的边界值为已知。

    这里给出的是方程的守恒形式,这是因为采用有限体积法建立离散方程时,必须使用守恒形式。应用有限体积法求解方程所对应的对流-扩散问题,主要步骤如下:

    1.        在计算区域内生成计算网格,包括节点及其控制体积。

    2.        将守恒型的控制方程在每个控制体积上做积分(积分时要用到界面处未知量Φ及其导数的插值计算公式,即离散格式),得到离散后的关于节点未知量的代数方程组。

    3.        求解代数方程式,得到个计算节点的Φ值。

    生成计算网格

    有限体积法的第一步是将计算域划分为离散的控制体积,在点A和点B之间的空间域上放置一系列的节点,将控制体积的边界(面)取在两个节点中间得位置,这样,每个节点由一个控制体积所包围。

    建立离散方程

    有限体积法关键一步是在控制体积上积分控制方程,以在控制体积节点上产生离散的方程,对一维模型方程,在图中所示的控制体积上坐积分,有:

    其中△V是控制体积的体积值,当控制体很微小时,△V可以表示为△x•A,这里A是控制体积界面的面积,从而有:

    上式中对流项和扩散项均已转化为控制体积界面上的值,有限体积法最显著的特点之一就是老师方程中具有明确的物理插值,即界面的物理量要通过插值的方式由节点的物理量表示。

    为了建立所需要形式的离散方程,我们需要找出如何表示界面e和w处的ρ、u、Γ、Φ和,有限体积法规定,ρ、u、Γ、Φ和等物理量均是在节点处定义和计算的。因此,为了计算界面上的这些物理参数,需要有一个物理参数在节点间的近似分布,可以想象,线性近似是可用来计算界面特性值的最直接、也是最简单的方式。这种分布叫做中心差分。如果网格是均匀的,则单个物理参数的线性插值结果是:

    与梯度项相关的扩散通量的线性插值结果是:

    对于源项S,它通常是时间和物理量Φ的函数,为了简化处理,经常将S转化为如下线性方式:

    其中Sc是常数,Sp是随时间和物理量变化的项,将以上线性插值的结果代回方程可得:

    整理后得:

    记作:

    对于一维问题,控制体积界面e和w处的面积A均为1,即单位面积,△V=△X.

    离散方程的求解

    为了求解所给出的流体流动问题,必须在整个计算域的每个节点上建立离散方程,从而每个节点有一个对应的方程,这些方程组成一个含有节点未知量的线性代数方程组,求解这个方程组就可以得到物理量在各节点处的值。理论上,任何可用于求解代数方程组的方法:Gauss消去法都可以完成上述任务。

    展开全文
  • 有限体积法源代码

    热门讨论 2013-02-28 08:46:01
    有限体积法 源代码 包括算例 网格刨分 各种方法
  • 有限差分法、有限元法、有限体积法等离散方法介绍
  • 这一篇主要是有限体积法笔记的后续部分及第一个算例-Testing the Order of FV-Approximations,后附MATLAB代码。下一篇应该是剩下的两个算例。*顺便更正:上一篇中Simpson' rule分母误写为2,应该为6算例验证-比较...

    e00c71d41b46423b501bee023306c0e8.png

    这一篇主要是有限体积法笔记的后续部分及第一个算例-Testing the Order of FV-Approximations,后附MATLAB代码。下一篇应该是剩下的两个算例。

    *顺便更正:上一篇中Simpson' rule分母误写为2,应该为6

    算例验证-比较重要的结论

    FVM的精度取决于以下三个因素:

    1. 格心外其它点上的插值精度
    2. 积分近似方法的精度
    3. 微分方程中各项导数的近似精度
    4. (当然,除了方法对精度的影响,网格的影响是显然的)

    为了提高FVM的精度,仅仅提高其中某一步的精度(误差阶数)并不能提高结果的总体精度,这需要各近似步骤之间达到合适的平衡。

    笔记:

    61bda247ad154acd8f515248d407fa83.png

    97bcc57c8af8159f8939d55ec9b15090.png

    efea0010cd562b4af490c37e054f61ae.png

    4b8a1c02e0400fe3ad2536936d528342.png

    3afad68cf47512d7dc440f8717b9d80d.png

    404df13a0119e0852751796317e85cf4.png

    c74e4aae24622d72b9d32326bced1763.png

    85d9cdc1d031716da78b47192fb0a215.png

    70f9893be3200d02a7b7fb1390c35999.png

    算例 1:Testing the Order of FV-Approximations

    main_interpolation

    2阶CDS和4阶CDS

    close all
    clear all
    
    %% pre processing
    % objective function
    fun='cosine';
    
    % required point
    P=[0;0.5];
    switch fun
        case 'poly'
            Exact=poly(P);
        case 'cosine'
            Exact=cosine(P);
    end
    
    %% calculating
    % calculating the interpolation approximation
    CD2=zeros(1,7);
    CD4=zeros(1,7);
    for i=1:7
        CD2(i)=calculate(i,2,P,fun);
        CD4(i)=calculate(i,4,P,fun);   
    end
    
    %% post processing
    % value
    x=1./2.^(1:6);
    
    figure;
    loglog(x,CD2(2:end),'--k','LineWidth',1.5);hold on
    loglog(x,CD4(2:end),'-.k','LineWidth',1.5);hold on
    loglog(x,Exact*ones(1,6),'k','LineWidth',2);
        set(gca,'FontSize',12);
        legend({'CD2','CD4','Exact'},'Location','west','FontSize',14)
        xlabel('Delta x','FontSize',16)
        ylabel('Variable value','FontSize',16)
        grid on
        hold off
    
    % error
    err2=abs((CD2(2:end)-Exact)/Exact*100);
    err4=abs((CD4(2:end)-Exact)/Exact*100);
    
    figure;
    loglog(x,err2,'--k','LineWidth',1.5);hold on
    loglog(x,err4,'-.k','LineWidth',1.5);
        set(gca,'FontSize',12);
        legend({'CD2','CD4'},'Location','northwest','FontSize',14)
        xlabel('Delta x','FontSize',16)
        ylabel('Error(%)','FontSize',16)
        grid on
        hold off

    main_integral approximation_midpoint rule

    midpoint rule下的积分近似(2阶误差)

    clear all
    close all
    
    %% pre processing
    % objective function
    fun='poly';
    
    % exact solution
    exapoly=1;
    exacosine=1+sin(1);
    switch fun
        case 'poly'
            exa=exapoly;
        case 'cosine'
            exa=exacosine;
    end
    
    %% integral approximation
    CD2=zeros(1,7);
    CD4=zeros(1,7);
    Exact=zeros(1,7);
    for i=1:7
        dx=1/2^(i-1);
        dy=1/2^(i-1);
        num=2^(i-1);
        node=[zeros(1,num);((1:num)*2-1)/2^i];
        for j=1:num
            node2(j)=calculate(i,2,node(:,j),fun);
            node4(j)=calculate(i,4,node(:,j),fun);
        end 
        % approximate integral using midpoint rule
        CD2(i)=node2*(ones(num,1)*dy);
        CD4(i)=node4*(ones(num,1)*dy);
        switch fun
            case 'poly'
                Exact(i)=poly(node)*(ones(num,1)*dy);
            case 'cosine'
                Exact(i)=cosine(node)*(ones(num,1)*dy);
        end
        
    end
    
    %% post processing
    ...

    main_integral approximation_Simpson's rule

    Simpson's rule下的积分近似(4阶误差)

    clear all
    close all
    
    %% pre processing
    % objective function
    fun='cosine';
    
    % number of spacing levels
    N=6;
    
    % exact solution
    exapoly=1;
    exacosine=1+sin(1);
    switch fun
        case 'poly'
            exa=exapoly;
        case 'cosine'
            exa=exacosine;
    end
    
    %% integral approximation
    CD2=zeros(1,N);
    CD4=zeros(1,N);
    Exact=zeros(1,N);
    for i=1:N
        dx=1/2^(i-1);
        dy=1/2^(i-1);
        num=2^(i-1);
        node=[zeros(1,num);((1:num)*2-1)/2^i];
        vertex=[zeros(1,num+1);((0:num)*2)/2^i];
        for j=1:num
            node2(j)=calculate(i,2,node(:,j),fun);
            node4(j)=calculate(i,4,node(:,j),fun);
        end 
        for j=1:(num+1)
            vertex2(j)=calculate(i,2,vertex(:,j),fun);
            vertex4(j)=calculate(i,4,vertex(:,j),fun);
        end 
        switch fun
            case 'poly'
                Exactn=poly(node);
                Exactv=poly(vertex);
            case 'cosine'
                Exactn=cosine(node);
                Exactv=cosine(vertex);
        end
        % approximate integral using Simpson's rule
        for j=1:num
            CD2(i)=CD2(i)+dy/6*(vertex2(j)+4*node2(j)+vertex2(j+1));
            CD4(i)=CD4(i)+dy/6*(vertex4(j)+4*node4(j)+vertex4(j+1));
            Exact(i)=Exact(i)+dy/6*(Exactv(j)+4*Exactn(j)+Exactv(j+1));
        end
        
    end
    
    %% post processing
    ...

    functions

    %% 计算多项式函数
    function [value]=poly(position)
    x=position(1,:);
    y=position(2,:);
    value=-2*x+3*x.^2-7*x.^3+x.^4+5*y.^4;
    end
    
    %% 计算三角函数
    function [value]=cosine(position)
    x=position(1,:);
    y=position(2,:);
    value=cos(x)+cos(y);
    end
    
    %% 中心差分
    function [phi]=CDS(order,position,value)
    switch order
        case 2
            la=(position(2)-position(1))/(position(3)-position(1));
            phi=value(3)*la+value(1)*(1-la);
        case 4
            % only for uniform grid
            phi=(9*value(2)+9*value(4)-value(1)-value(5))/16;
    end
    end
    
    %% 主要计算过程(其实可以和CDS合到一起)
    function [appP]=calculate(level,order,P,fun)
    % spacing
    dx=1/2^(level-1);
    dy=1/2^(level-1);
    
    % needed nodes
    posE(1,:)=P(1)+(1:(order/2))*dx-dx/2;
    posE(2,:)=P(2);
    switch fun
        case 'poly'
            valE=poly(posE);
        case 'cosine'
            valE=cosine(posE);
    end
    
    posW(1,:)=P(1)-(1:(order/2))*dx+dx/2;
    posW(2,:)=P(2);
    switch fun
        case 'poly'
            valW=poly(posW);
        case 'cosine'
            valW=cosine(posW);
    end
    
    %% calculating approximation
    pos=[flip(posW(1,:)),P(1),posE(1,:)];
    val=[flip(valW),0,valE];
    
    appP=CDS(order,pos,val);
    
    end

    Solutions (figures)

    跟Ferziger书上结果一致

    078f968bc3aa4387cfd4f543eebec483.png
    多项式函数 插值

    bd5e89ed44253dd6d08babb7ceb71989.png
    多项式函数 插值误差

    d016ca07e6ba7d3a4cbabfc9ee1a08a3.png
    三角函数 插值

    a33ffd16b73eb884ac54a2acbecbfe68.png
    三角函数 插值误差

    a7b1c3be11ffa140dd5415693e988544.png
    多项式函数 中点法则积分近似

    2ed9f067a6ac6ba1e72a7ca8e55e0141.png
    多项式函数 中点法则积分近似误差

    09ebf9a5e2d97a07034f07f1fafd4a46.png
    三角函数 中点法则积分近似

    ead6fbd78ea30b0b49ff69b0b79f7c80.png
    三角函数 中点法则积分近似误差

    925605495f739dc08d3ebd3d9bca91d1.png
    多项式函数 Simpson法则积分近似

    d1a10b93b54b64dfd937ca2a7361baf7.png
    多项式函数 Simpson法则积分近似误差

    b682e46ecdd1be1873a85c68b8ffef6a.png
    三角函数 Simpson法则积分近似

    fc7e29d7ea11bb8f5c5126c63cb41c20.png
    三角函数 Simpson法则积分近似误差
    展开全文
  • 角度离散对有限体积法辐射计算的影响,吴明,刘立军,采用有限体积法计算四边形区域内吸收发射性介质的辐射传递问题,研究了区域内部角离散数目和边界处像素点数目对辐射计算结果精度
  • 根据能量守恒定律和傅立叶定律,采用有限体积法对二维稳态平面温度场进行求解,合理地确定控制体的大小,建立了平面稳态温度场问题的线性方程组及其离散模型。结果表现了平面温度场三角形单元变分计算公式的实际物理...
  • 有限体积法的粘性力矩问题及其改进,杨振,,本构方程的对称性是流体力矩平衡的必要条件,然而其粘性项的一部分合力为零,在Navier-Stokes方程中它们相互抵消,因而常规有限体积�
  • 为了研究弹性力学中三种典型的有限体积法方案(FVM1,FVM2和FVM3)的精度差异问题,采用理论推导和算例验证的方法,建立了基于有限体积法原理的力平衡积分方程,利用三角形网格进行离散化,在计算出内部单元和边界单元对力...
  • 有限体积法求解欧拉方程,翼型为naca0012
  • 有限体积法 求解 一维二维对流扩散问题 ,一维稳态问题,采用中心差分并与解析解比较。一维稳态 乘方格式
  • 基于无结构网格有限体积法的Roe隐格式探讨,杨彬,,本文采用非结构网格有限体积法,尝试建立适合河流、湖泊和近海浅水水流流动的Roe隐式离散格式的数学模型。对于隐式离散方程组采用
  • 基于有限体积法的二维溃坝水流模拟研究进展,刘铁锤,蔡华钦,溃坝计算可用于制定地区应急预案,加强区域防洪减灾,并用于后期的环境与生态评估。回顾和总结了国内外基于有限体积法对溃坝水流
  • 做了剩下有限体积法(Ferziger, Ch.4[1])另外两个Example(测试数值扩散的那一算例,书里似乎并没有给出CDS矩阵奇异下的求解方法,我直接按MATLAB里AB求解,结果并不完全一致):Scalar Transport in a Known ...

    做了剩下有限体积法(Ferziger, Ch.4[1])另外两个Example(测试数值扩散的那一算例,书里似乎并没有给出CDS矩阵奇异下的求解方法,我直接按MATLAB里AB求解,结果并不完全一致):

    1. Scalar Transport in a Known Velocity Field
    2. Testing the Numerical Diffusion

    果然FVM比FDM麻烦不少,一开始还是按面向过程写,中途感到不对劲于是写成类了。写到一半去搜了下Openfoam[2]发现自己写的类定义(思路)和它居然差不太多,挺巧。

    *网页端文章还可以把背景图片搞成这种花里胡哨的可还行

    问题

    Example 2 Scalar Transport in a Known Velocity Field

    考虑标量输运方程

    equation?tex=%5Cint_S%5Crho%5Cphi%5Cbm%7Bv%7D%5Ccdot%5Cbm%7Bn%7D+%5C+dS%3D%5Cint_S%5CGamma%5Cbm%7B%5Cnabla%7D%5Cphi%5Ccdot%5Cbm%7Bn%7D%5C+dS

    算一个给定速度场(驻点附近的无粘流动,

    equation?tex=u_x%3Dx%2Cu_y%3D-y ,流线如图)、边界条件等参数下,标量场
    equation?tex=%5Cphi 的分布情况。

    边界条件具体为:

    • 左侧(wall)和上侧(inlet)为第一类边界,左侧
      equation?tex=%5Cphi 为从0到1线性分布(如图),上侧
      equation?tex=%5Cphi%3D0
    • 右侧(outlet)和下侧(symmetric boundary)为第二类边界,
      equation?tex=%5Cphi_n%3D0

    要求扩散项采用CDS,对流项分别采用CDS和UDS计算。

    984f640a916dbb59b407d61cd8d488d8.png
    几何与边界条件

    Example 3 Testing the Numerical Diffusion

    考虑具有阶梯剖面(间断)斜向均匀流中的对流问题(无扩散项)

    4b67e042f84664c8d67bf3043de46941.png

    面对这一问题,CDS构造出的系数矩阵A主对角元均为0,矩阵奇异难以求解。总之就是验证CDS会产生震荡以及UDS通过数值扩散避免了震荡的发生。(TVD格式消除震荡暂时没考虑)

    基本思路

    设计了三个类:

    • solver
    • grid
    • cell

    solver类的目标在于针对设定的物理域、计算域、边界条件、各类参数及求解方法等,建立一个计算求解方案。它包含了grid类(网格),进一步的,grid类中每一个网格单元储存为cell对象。

    我是按照自底向上,即cell->grid->solver的过程编写的,毕竟只是解一两道题,而非编写一套CFD程序,不需要太宏观的架构。

    参数有两套存储形式,一是在solver类中同意存储所有网格单元上的各类信息,二是分别存储到对应的网格单元上。于是差分格式-积分近似这类计算在网格单元(cell类)中完成,而边界条件的处理及总体系数矩阵的组集则在solver类中计算。

    此外,写边界条件设置时需要脑子清醒(比如CDS改成单边近似这种事)。

    一些结果与讨论

    先分别给出40*40网格下,扩散率

    equation?tex=%5CGamma 分别取0.01和0.001时,CDS求解出的标量场等值线图(对照书上原图应该基本是一模一样的,一丝区别是书上还把边界值给画了进去,这里我就懒得整了)

    24d7990417bf013ed932e2e74ed7562f.png
    40*40 Gamma=0.01 CDS Isolines

    8b4321ce465f0fa6fc45715c27d773e3.png
    40*40 Gamma=0.001 CDS Isolines

    顺便看看40*40网格下的系数(也是稀疏)矩阵A,长得也符合预期:

    672fc9371691f2f63cf3c1fc709eb4e9.png
    40*40 Gamma=0.01 CDS sparse matrix A

    下面两张分别是320*320精细网格下的CDS和更直观的3d图示:

    3a328fd228d5425394850a10b46d5df3.png

    7cda006510e0f98c53a998b9cfbb0d40.png

    在粗糙网格(10*10)下,CDS出现了震荡解,而UDS的不震荡特性也得到了验证:

    3688e61e0b4c0dc5677774ab39d6e010.png
    10*10 Gamma=0.001 CDS oscillation (看颜色深浅深浅深浅)

    fc9bbddb3e5926ec0faebdc10e7037d5.png
    10*10 Gamma=0.001 UDS no oscillation

    然后是Example 2 左侧边界总通量(随网格量变化)的计算结果,与书里一致(书上通量值莫名其妙乘了个100?)

    337aad538b20d341b2600c75a095add4.png

    Example 3 没太多东西,改改参数就可以用这代码算了,一些结果:

    09b7b9e5ee140d053dfeec4f8923cb90.png
    UDS数值扩散避免震荡

    38f71129e1c6144e7763ff37a240fe8a.png
    UDS x=0.5处 phi分布 (及随网格量变化)

    f4c2e4eb9a175adc5fcf6c6e757f8356.png
    CDS x=0.5处 phi分布 (及随网格量变化)

    Codes

    class solver (求解器)

    classdef solver<handle
        %% Scalar Transprot solver 
        % using Finite Volume Method on uniform grid
        properties 
            grid            % 
            size            % Nx*Ny
            position_mp     % position of surface mid points
            velofield_mp    % velocity field (on surface mid points)
            massflux_mp     % mass flux through cell surface
            
            Q               % source term
            A               % coefficient matrix (sparse)
            phi             % solution
            
            rho             % constant density
            gamma           % constant diffusion coefficient
            
            drawx
            drawy
            drawp
            draw_on
            
            west_flux        % total flux of phi through west boundary
            
        end 
    	methods
            %% constructor
            function obj=solver(rho,gamma,size,domain,boundary,spacing)
                %% pre initialization %%
                switch nargin
                    case 2
                        obj.rho=rho;
                        obj.gamma=gamma;
                    case 4
                        obj.rho=rho;
                        obj.gamma=gamma;
                        obj.grid=uniformgrid(size,domain);
                    case 5
                        obj.rho=rho;
                        obj.gamma=gamma;
                        obj.grid=uniformgrid(size,domain,boundary);
                    case 6
                        obj.rho=rho;
                        obj.gamma=gamma;
                        obj.grid=uniformgrid(size,domain,boundary,spacing);
                end
                
                %% object initialization %%
                if nargin>=4
                    obj.size=obj.grid.Nx*obj.grid.Ny;
                    
                    % calculate information on surface mid points
                    obj.position_mp=cell(obj.size,1);
                    obj.velofield_mp=cell(obj.size,1);
                    obj.massflux_mp=zeros(obj.size,4);
                    n=[-1,0;0,-1;1,0;0,1]; % surface unit vector
                    s=[obj.grid.dy,obj.grid.dx,obj.grid.dy,obj.grid.dx];
                    for i=1:obj.size
                        obj.position_mp{i}=[...
                            obj.grid.cells{i}.position(1)-obj.grid.dx/2,...
                            obj.grid.cells{i}.position(1),...
                            obj.grid.cells{i}.position(1)+obj.grid.dx/2,...
                            obj.grid.cells{i}.position(1);...
                            obj.grid.cells{i}.position(2),...
                            obj.grid.cells{i}.position(2)-obj.grid.dy/2,...
                            obj.grid.cells{i}.position(2),...
                            obj.grid.cells{i}.position(2)+obj.grid.dy/2];
                        obj.velofield_mp{i}=velocity(obj.position_mp{i});
                        for j=1:4
                            obj.massflux_mp(i,j)=obj.rho*n(j,:)*obj.velofield_mp{i}(:,j)*s(j);
                            obj.grid.cells{i}.massflux(1,j)=obj.rho*n(j,:)*obj.velofield_mp{i}(:,j)*s(j);
                        end
                    end
                    
                    % clear useless information
                    obj.position_mp=[]; 
                    obj.velofield_mp=[];
                    obj.massflux_mp=[];
                    obj.Q=zeros(obj.size,1);
                    obj.draw_on=0;
                    
                end
                
            end
            
            %% calculate
            function new_calculate(obj)
                for i=1:obj.size
                    obj.grid.cells{i}.new_calculate();
                end
                obj.Q=zeros(obj.size,1);
                obj.A=[];
                obj.phi=[];
            end
            
            function calculate(obj,method)
                % diffusion without boundary condition
                % method{1} method for diffusion (default CDS)
                % method{2} method for convection (UDS or CDS)
                for i=1:obj.size
                    obj.grid.cells{i}.diffusion(method{1},obj);
                end
                
                % convection without boundary condition
                for i=1:obj.size
                    obj.grid.cells{i}.convection(method{2});
                end
                
                % setting boundary condition
                obj.boundary();
                
                % construct sparse matrix A
                %% delete it!!! %% 一开始偷懒直接构造矩阵,后来发现太耗空间了
    %             obj.A=zeros(obj.size);
    %             for i=1:obj.size
    %                 p=c2p(i,obj.grid.Ny);
    %                 x=p(1);y=p(2);
    %                 obj.A(i,i)=obj.grid.cells{i}.A(5);
    %                 ind=[p2c([x-1,y],obj.grid.Ny),...
    %                     p2c([x,y-1],obj.grid.Ny),...
    %                     p2c([x+1,y],obj.grid.Ny),...
    %                     p2c([x,y+1],obj.grid.Ny)];
    %                 for j=1:4
    %                     if obj.grid.cells{i}.A(j)~=0
    %                         obj.A(i,ind(j))=obj.grid.cells{i}.A(j);
    %                     end
    %                 end
    %             end
                %% sparse %%
    %             S=sparse(obj.A);
    %             obj.A=[];
                
                
                count=1;
                for i=1:obj.size
                    p=c2p(i,obj.grid.Ny);
                    x=p(1);y=p(2);
                    % diagonal element
                    sparsei(count)=i;
                    sparsej(count)=i;
                    sparseA(count)=obj.grid.cells{i}.A(5);
                    count=count+1;
                    % other element
                    ind=[p2c([x-1,y],obj.grid.Ny),...
                        p2c([x,y-1],obj.grid.Ny),...
                        p2c([x+1,y],obj.grid.Ny),...
                        p2c([x,y+1],obj.grid.Ny)];
                    for j=1:4
                        if obj.grid.cells{i}.A(j)~=0
                            sparsei(count)=i;
                            sparsej(count)=ind(j);
                            sparseA(count)=obj.grid.cells{i}.A(j);
                            count=count+1;
                        end
                    end
                end
                obj.A=sparse(sparsei,sparsej,sparseA);
                sparsei=[];
                sparsej=[];
                sparseA=[];
                
                obj.phi=obj.Aobj.Q;
            end
            
            %% boundary setting
            function boundary(obj)
                % setting boundary for each side
                for i=1:obj.size
                    boundtype=obj.grid.cells{i}.boundtype;
                    for j=1:4
                        switch boundtype{j} 
                            % boundary {1,2,3,4} for {w,s,e,n}
                            case 'first'
                                obj.grid.cells{i}.A(j)=...
                                    obj.grid.cells{i}.Ac(j)+...
                                    obj.grid.cells{i}.Ad(j)*2;
                                obj.grid.cells{i}.A(5)=...
                                    obj.grid.cells{i}.A(5)-...
                                    obj.grid.cells{i}.A(j);
                                obj.Q(i)=obj.Q(i)-...
                                    obj.grid.cells{i}.A(j)*...
                                    obj.grid.cells{i}.boundvalue{j};
                                obj.grid.cells{i}.A(j)=0;
                            case 'second'
                                obj.grid.cells{i}.A(j)=0;
                                obj.grid.cells{i}.A(5)=...
                                    obj.grid.cells{i}.A(5)-...
                                    obj.grid.cells{i}.A(j);
                            case '0'
                                obj.grid.cells{i}.A(j)=...
                                    obj.grid.cells{i}.Ac(j)+...
                                    obj.grid.cells{i}.Ad(j);
                                obj.grid.cells{i}.A(5)=...
                                    obj.grid.cells{i}.A(5)-...
                                    obj.grid.cells{i}.A(j);
                        end
                        
                    end
                end
            end
            
            %% west flux calculate
            function [wf]=westflux(obj)
                wf=0;
                for i=1:obj.grid.Ny
                    Fd=obj.gamma*obj.grid.dy*(obj.grid.boundary_value{1}(i)-...
                        obj.phi(i))/(-obj.grid.dx/2);
    %                 Fc=0.5*obj.massflux(i,1)*(obj.boundary_value{1}(i)+...
    %                     obj.cell{i}.value); 
                          % convection flux is zero on this boundary
                    wf=wf+(-Fd);%+Fc);
                end
                obj.west_flux=wf;
            end
            
            
            %% Post processing
            % draw figure
            function draw(obj,drawoption)
                % pre calculate
                if obj.draw_on~=1
                    obj.drawx=zeros(obj.grid.Nx,obj.grid.Ny);
                    obj.drawy=zeros(obj.grid.Nx,obj.grid.Ny);
                    obj.drawp=zeros(obj.grid.Nx,obj.grid.Ny);
                    for i=1:obj.grid.Nx
                        for j=1:obj.grid.Ny
                            obj.drawx(i,j)=obj.grid.cells{p2c([i,j],obj.grid.Ny)}.position(1);
                            obj.drawy(i,j)=obj.grid.cells{p2c([i,j],obj.grid.Ny)}.position(2);
                            obj.drawp(i,j)=obj.phi(p2c([i,j],obj.grid.Ny));
                        end
                    end
                    
                    obj.draw_on=1;
                end
                
                % draw
                figure;
                hold on
                for i=1:length(drawoption)
                    switch drawoption{i}
                        case 'contour'
                            levels=0.05:0.1:0.95;
                            contour(obj.drawx,obj.drawy,obj.drawp,levels,'-k',...
                                'ShowText','off','LabelSpacing',2000);
                                grid on
                                axis equal tight
                        case 'pcolor'
                            pcolor(obj.drawx,obj.drawy,obj.drawp);
                                shading flat
                                axis equal tight
                                
                    end
                end
    %             set(gca,'FontSize',12);
                hold off;
            end
            
    
        end
    end

    class uniformgrid (均匀正交网格)

    classdef uniformgrid<handle
        properties 
            cells       % store cell objects
            boundary_type
            boundary_value
            
            Nx          % number of nodes in x direction
            Ny          % number of nodes in y direction
            dx          % x direction spacing
            dy          % y direction spacing
            
            domain      % physical domain
            
        end 
    	methods
            %% constructor
            function obj=uniformgrid(size,domain,boundary,spacing)
                %% pre initialization %%
                switch nargin
                    case 2
                        obj.Nx=size(1);
                        obj.Ny=size(2);
                        obj.domain=domain;
                        obj.dx=(domain(2)-domain(1))/obj.Nx;
                        obj.dy=(domain(4)-domain(3))/obj.Ny;
                        obj.boundary_type=... % default
                            {'first','first','first','first'};
                        obj.boundary_value=... % default
                            {zeros(1,obj.Ny),zeros(1,obj.Nx),zeros(1,obj.Ny),zeros(1,obj.Nx)};
                        
                    case 3
                        obj.Nx=size(1);
                        obj.Ny=size(2);
                        obj.domain=domain;
                        obj.dx=(domain(2)-domain(1))/obj.Nx;
                        obj.dy=(domain(4)-domain(3))/obj.Ny;
                        obj.boundary_type=boundary{1};
                        obj.boundary_value=boundary{2};
                
                    case 4
                        obj.Nx=size(1);
                        obj.Ny=size(2);
                        obj.domain=domain;
                        obj.boundary_type=boundary{1};
                        obj.boundary_value=boundary{2};
                        obj.dx=spacing(1);
                        obj.dy=spacing(2);
                end
                
                %% Object Initialization %%
                % construct cells
                obj.cells=cell(1,obj.Nx*obj.Ny);
                for i=1:(obj.Nx*obj.Ny) % can be optimize
                    p=c2p(i,obj.Ny);
                    x=p(1);y=p(2);
                    position=[domain(1)+(x-1/2)*obj.dx,domain(3)+(y-1/2)*obj.dy];
                    neighbour=[p2c([x-1,y],obj.Ny),...
                        p2c([x,y-1],obj.Ny),...
                        p2c([x+1,y],obj.Ny),...
                        p2c([x,y+1],obj.Ny)];
                    boundtype={'0','0','0','0'};
                    boundvalue=cell(1,4);
                    if x==1
                        boundtype{1}=obj.boundary_type{1};
                        boundvalue{1}=obj.boundary_value{1}(y);
                    end
                    if x==obj.Nx
                        boundtype{3}=obj.boundary_type{3};
                        boundvalue{3}=obj.boundary_value{3}(y);
                    end
                    if y==1
                        boundtype{2}=obj.boundary_type{2};
                        boundvalue{2}=obj.boundary_value{2}(x);
                    end
                    if y==obj.Ny
                        boundtype{4}=obj.boundary_type{4};
                        boundvalue{4}=obj.boundary_value{4}(x);
                    end
                    boundary{1}=boundtype;
                    boundary{2}=boundvalue;
                    obj.cells{i}=rectcell(i,position,neighbour,boundary);
                end
            end
            
            
    	end
    
    end

    class rectcell (矩形网格单元)

    classdef rectcell<handle
        %% node on cell center
        properties 
            index       % index of the cell
            neighbour   % index of the neighbour cells ([W S E N])
            position    % position of the cell center
            value       % value at the cell center
            boundtype   % boundary type 
            boundvalue  % boundary value
            
            massflux    % mass flux
            
            Ac          % coefficient convection [AWc ASc AEc ANc APc]
            Ad          % coefficient diffusion [AWd ASd AEd ANd APd]
            A           % coefficient [AW AS AE AN AP]
            
        end 
    	methods
            %% contructor
            function obj=rectcell(index,position,neighbour,boundary)
                %% pre initialization %%
                switch nargin
                    case 1
                        obj.index=index;
                        obj.boundtype=cell(1,4); % default
                        obj.boundvalue=cell(1,4); % default
                        obj.massflux=zeros(1,4); % default
                        
                    case 2
                        obj.index=index;
                        obj.position=position;
                        obj.boundtype=cell(1,4); % default
                        obj.boundvalue=cell(1,4); % default
                        obj.massflux=zeros(1,4); % default
                        
                    case 3
                        obj.index=index;
                        obj.position=position;
                        obj.neighbour=neighbour;
                        obj.boundtype=cell(1,4); % default
                        obj.boundvalue=cell(1,4); % default
                        obj.massflux=zeros(1,4); % default
                        
                    case 4
                        obj.index=index;
                        obj.position=position;
                        obj.neighbour=neighbour;
                        obj.boundtype=boundary{1};
                        obj.boundvalue=boundary{2};
                        obj.massflux=zeros(1,4); % default
                        
                end
                
                obj.A=zeros(1,5);
                obj.Ac=zeros(1,5);
                obj.Ad=zeros(1,5);
            end
            
            %% new
            function new_calculate(obj)
                obj.A=zeros(1,5);
                obj.Ac=zeros(1,5);
                obj.Ad=zeros(1,5);
            end
            
            %% approximation of convection and diffusion
            %% convection %%
            function convection(obj,method)
                switch method
                    case 'UDS'
                        AcW=min(obj.massflux(1),0); 
                        AcS=min(obj.massflux(2),0);
                        AcE=min(obj.massflux(3),0);
                        AcN=min(obj.massflux(4),0);
                        AcP=-(AcW+AcE+AcS+AcN);
                        obj.Ac=[AcW,AcS,AcE,AcN,AcP];
                        
                    case 'CDS'
                        AcW=obj.massflux(1)/2; 
                        AcS=obj.massflux(2)/2;
                        AcE=obj.massflux(3)/2;
                        AcN=obj.massflux(4)/2;
                        AcP=-(AcW+AcE+AcS+AcN); % =0
                        obj.Ac=[AcW,AcS,AcE,AcN,AcP];
                        
                end
                
            end
            
            %% diffusion %%
            function diffusion(obj,method,solver) % default CDS
                dx=solver.grid.dx;
                dy=solver.grid.dy;
                gamma=solver.gamma;
                
                if method=='CDS'
                    AdE=-gamma*dy/dx;
                    AdW=-gamma*dy/dx;
                    AdN=-gamma*dx/dy;
                    AdS=-gamma*dx/dy;
                    AdP=-(AdW+AdE+AdS+AdN); 
                    obj.Ad=[AdW,AdS,AdE,AdN,AdP];
                end
    
            end
            
        end
    end

    test(Example 2)(Example 3 需要重设边界条件等)

    clear all;
    close all;
    
    %% settings
    % grid numbers
    Nx=40;Ny=40;
    % physical domain
    domain=[0,1,0,1];
    % parameters
    rho=1;
    gamma=0.001;
    
    % velocity field is defined by funciton velocity(...)
    
    % boundary type & boundary value
    type={'first','second','second','first'}; 
    value={1-((1:Ny)-1/2)*1/Ny,zeros(1,Nx),zeros(1,Ny),zeros(1,Nx)}; 
    bound={type,value}; 
    % method for diffusion and convection
    method={'CDS','CDS'}; 
    
    %% solve
    problem=solver(rho,gamma,[Nx,Ny],domain,bound);
    problem.calculate(method);
    
    %% solution 
    % scalar field
    problem.draw({'pcolor','contour'});
    
    % show the sparse matrix A
    % figure;
    % spy(problem.A(1:(2*min(Nx,Ny)),1:(2*min(Nx,Ny))));
    
    %% calculate west flux as grid spacing varying
    cal_wf=1;
    
    if cal_wf
        level=4;
        size=(2.^(0:(level-1))*10).^2;
        wf=zeros(2,level);
    
        for i=1:level
            Nx=sqrt(size(i));
            Ny=Nx;
            value={1-((1:Ny)-1/2)*1/Ny,zeros(1,Nx),zeros(1,Ny),zeros(1,Nx)};
            bound={type,value}; 
            problem=solver(rho,gamma,[Nx,Ny],domain,bound);
        
            problem.calculate({'CDS','CDS'});
            wf(1,i)=problem.westflux();
        
            problem.new_calculate(); % delete previous results
            problem.calculate({'CDS','UDS'});
            wf(2,i)=problem.westflux();
        end
        % show the convergence of west flux as grid spacing decreasing
        figure;
        semilogx(size(1:end),wf(1,1:end),'-k');hold on
        semilogx(size(1:end),wf(2,1:end),'--k');
            grid minor
            set(gca,'FontSize',14);
            axis([100,120000,0.01,0.014])
            legend(['CDS';'UDS'])
            xlabel('Number of CVs','FontSize',16)
            ylabel('Total flux (west)','FontSize',16)
            
        % calculate the error (via Richardson extrapolation)
        ...
        
    end
    

    functions

    function [v]=velocity(position) % 给定速度场计算
    %% calculate velocity at given point
    v(1,:)=position(1,:);
    v(2,:)=-position(2,:);
    % v(1,:)=ones(1,size(position,2));
    % v(2,:)=ones(1,size(position,2));
    end
    
    function [physical]=c2p(computational,Nj) 
    %% computational domain index to physical domain
    physical(1)=ceil(computational/Nj);
    physical(2)=computational-(physical(1)-1)*Nj;
    end
    
    function [computational]=p2c(physical,Nj)
    %% physical domain index to computational domain
    computational=(physical(1)-1)*Nj+physical(2);
    end

    参考

    1. ^Ferziger, J. H., et al. (2002). Computational methods for fluid dynamics, Springer.
    2. ^FVM in CFD 学习笔记_第7章_OpenFOAM和uFVM中的有限体积网格 https://blog.csdn.net/meiguanhua/article/details/104756024/
    展开全文
  • 有限体积法 求解 一维二维对流扩散问题 ,一维稳态问题,采用中心差分并与解析解比较。
  • 基于改进有限体积法的三维注塑成型充模过程数值模拟,严波,李阳,分别根据界面两侧相邻有限体积流动剪切应力相等的原则,考虑到相邻有限体积内粘度系数的差异,改进传统的有限体积法,推导了改进
  • 根据煤粒瓦斯解吸放散的特点,基于达西定律,建立了圆柱形煤粒瓦斯放散微分方程,并运用有限体积法对方程进行处理。引入无因次准数,对瓦斯放散微分方程及运用有限体积法进行离散处理后的方程进行无因次化。通过编写...
  • 基于有限体积法的近岸海域水环境二维数值模拟,孙佳佳,,本文根据海域的水动力特征,采用三角形网格划分计算区域,采用有限体积法离散控制方程,建立非稳态平面二维水动力数学模型。并将
  • 为研究梯形巷道围岩散热的特点,依据能量守恒定律,建立二维非稳态导热积分方程,运用有限体积法分别建立内部节点和边界节点的热平衡方程,采用VB语言编写围岩温度场及散热量数值计算程序,并通过一个实例进行数值分析....
  • 基于有限体积法的水柱分离-弥合瞬变流建模研究,周领,王欢,水柱分离-弥合水锤的数值模拟一直是水力瞬变计算的一个难题。本文提出一种基于有限体积法的时空均为二阶精度的计算模型。该模型�
  • 基于有限体积法的三维注塑结晶过程模拟,严波,李阳,半结晶型塑料的注塑熔体流变行为改变了结晶形态的演化规律,基于剪切应力引起分子链取向并导致平衡熔点升高,建立了非等温剪切诱
  • 有限元法 有限差分法 有限体积法

    千次阅读 2020-03-22 19:58:46
    有限也叫有限单元(finite element method, FEM),是随着电子计算机的发展而迅速发展起来的一种弹性力学问题的数值求解方法。五十年代初,它首先应用于连续体力学领域—飞机结构静、动态特性分析中,用以求得...
  • 有限体积法求解二维双层半透明介质辐射-导热耦合传热,金子程,路义萍,采用基于有限体积法的开源C++程序库OpenFOAM,开发了求解二维双层半透明介质辐射-导热耦合传热问题的计算程序。在对计算模型验证的基

空空如也

空空如也

1 2 3 4 5 ... 15
收藏数 299
精华内容 119
关键字:

有限体积法