• 返回主页
  • OKS+NOM预习资料

    2017.07.04:经国旭邮件勘误(已修订):公式9和10最后一项遗漏1/2
    2017.06.27:更新完毕
    2017.06.21:补充部分内容
    2017.06.18:补充部分内容
    2017.05.23:创立页面


    OKS课程工具介绍:OpenFOAM是一个采用C++编写的CFD求解器,其主要的特点在于开源并且免费。开源的特性方便各类研究人员修改代码。国内CFD用户目前使用最常用的软件是ANSYS Fluent,OpenFOAM相对于ANSYS Fluent的优势在于可以方便的进行代码定制以及开发。ANSYS Fluent在使用的过程中,用户可以通过UDF来自定义函数,但ANSYS Fluent提供的接口有限。OpenFOAM用户则可以任意的修改代码进行任意的函数定义。但是OpenFOAM的学习过程非常缓慢并且资料非常少。
    NOM课程工具介绍:ANSYS ICEM是目前非常受欢迎的网格生成程序,但在使用ANSYS ICEM制作全六面体网格的时候需要非常烧脑的block划分技巧,市面上虽然存在一些ANSYS ICEM教程,但是这种block划分思想主要通过各类工程实践才能获得。

    本预习资料介绍
    为了让参加OpenFOAM Knowledge Share以及Not Only Mesh课程的学员在课堂上能跟上思路来获得最佳的学习效果,并鉴于CFD以及计算机语言较难上手的学习特性,所有学员(尤其是参加OKSS2的学员)务必在开课前详细充分预习本文内容。
    OKSS1:OKSS1课程主要偏向OpenFOAM应用,只需要非常基本的数学和编程技巧。参加OKSS1的学员主要需要预习Linux的基本操作,其他概念可简单理解;
    OKSS2:OKSS2课程主要讲解CFD方程在OpenFOAM中的实现以及OpenFOAM代码定制。需要扎实的高等数学基础和C++基础。本预习资料简单介绍CFD比较基本的概念且在大多数教科书中均有说明。但是不同的教科书针对的读者层次不同导致读起来较难。本预习资料尽可能的浅出。在预习本文内容后,课堂上将直接从课程列表中OKSS2中的有限体积中的高斯积分进行。若学员对此部分内容已经充分掌握,可略过。另外需要注意,OKS课堂偏向介绍CFD在C++中的实现,学员务必提前学习好C++中有关类的概念,不要着急运行OpenFOAM算理,望厚积薄发。
    NOMS1:NOMS1主要介绍一些网格的数学特性,对编程没有要求,请预习好本网页中的CFD理论即可;
    NOMS2:NOMS2直接使用ANSYS ICEM划分网格,不需要做预习;

    OKS以及NOM除了预习资料外,还会发放一些预先打印好的文献资料在课堂上分发,除此之外无PPT等任何资料。但是会每人发一个记事本,请充分做好笔记。课程费用约每小时150元,若已报名,一定要努力预习并做好笔记获得大于费用的成果!

    开始上学 :-)

    下面的内容已修改为默认展现(部分浏览器不兼容),可点击按钮收起方便浏览

    在本节,链接指向京东的书籍建议直接购买来进行透彻详细的学习,指向的连接可直接下载。

    在linux中有众多发行版,OKS课程选择比较接近windows使用习惯的Ubuntu系统,可将Ubuntu系统安装在虚拟机中,也可以选择双系统的方式。已经预装好OpenFOAM-4.x的Ubuntu虚拟机可点击下载。

    在Linux系统中,大部分操作在终端中进行,如下图中右半部分的窗口所示,且OpenFOAM的运行主要在终端中操作,熟悉终端中的操作会大大增加工作的效率。

    下述代码切换到用户文件夹下:

    cd ~

    用户如果不清楚什么是用户文件夹,键入下述代码可以列出当前的文件路径:

    pwd

    返回上一个文件夹:

    cd ..

    返回上面多个文件夹:

    cd ../../../

    创建一个空的名字为folder的文件夹:

    mkdir folder

    创建一个空的名字为file的文件:

    touch file

    复制file文件为file1:

    cp file file1

    将file1重命名为file2:

    mv file1 file2

    将folder下的file文件移动到folder1下:

    mv folder/file folder1/

    拷贝folder1文件夹并重新命名为folder2:

    cp -r folder1 folder2

    删除文件file:

    rm file

    删除文件夹folder:

    rm -r folder

    显示文件file的内容:

    less file

    其中可以用上箭头和下箭头进行导航,也可以用j和k导航,点击q退出。也可以用more来显示文件file的内容:

    more file

    使用more的时候,使用回车来换行。 有时用户并不想查看全部文件的内容,指向查看log文件中计算最后的时间步到哪了,可以用tail来显示文件最后的内容。比如如果打算查看file文件最后的10行:

    tail -10 file

    对于正在运行的算例,log文件可能一直在被写入,如果打算查看动态的log文件,可以使用:

    tail -f file

    这样会在终端持续的输出log文件。

    查找src/目录下header.H文件的路径:

    find src/ -name header.H

    查找src/目录下哪些文件具有关键词FOAM:

    grep -R "FOAM"

    查找log文件夹下是否存在FOAM关键词:

    grep -l FOAM log

    有的话输出log,没有的话不显示。

    上述内容是部分常用的Linux基本操作。对于入门的学员掌握一些基本操作就可以,若对Linux操作有更深的兴趣,可学习一些功能更强大的一些操作如awksed,以及一些命令行下的编辑器如VIM等。

    符号说明

    在所有东岳流体的文档(包括本文)中,将采用统一的符号标记。标量为斜体如压力,温度,湍流动能为: \begin{equation} p,T,k \end{equation} 矢量为加粗正体如速度,面矢量为: \begin{equation} \bfU,\bfS \end{equation} 因此,速度的分量为三个标量$u$,$v$,$w$。二阶张量和矢量相同。

    计算流体力学中的“流体”

    流体在通常泛指气体液体。不同于固体,流体在施加剪切力的时候,会发生形变。并且,流体没有固定的形态。流体在快速变形的时候,也会存在一种抵抗的力,这种抵抗的力一旦流动停止,也即消失。流体这种抗拒本身发生变形的力称之为粘性力。不同的流体存在不同的粘性力。有些流体的粘性很大,有些流体的粘性很小。进阶概念:在CFD中,气体液体固体可以用不同的相来称呼,虽然流体通常不包含固体,但是在CFD的某些模型中可以把大量固体颗粒构成的集(如沙子)认为是一种近似的流体,例如在双流体模型中,其中的固体集也被看做为流体。

    液体的典型特征为其对体积变化的巨大抵抗力。例如我们不能把一瓶矿泉水压缩成一瓶盖的大小。然而气体的典型特征为其对体积变化的抵抗力是相对小的。例如我们可以轻易地用手对一个封闭的空的矿泉水瓶进行压缩。气体这种压强和体积的关系可以用玻意尔定律来表示: \begin{equation} p_1 \cdot V_1 = p_2 \cdot V_2 \end{equation} 其中$p_1$和$p_2$表示$V_1$体积下以及$V_2$体积下的压强。

    气体除了压强和体积有关之外,还和温度有关。普适性来讲,有理想气体状态方程来描述这些变量之间的关系: \begin{equation} pV=nRT \end{equation} 其中$R$为理想气体常数,$T$为温度,$n$为物质的量。除了温度以外,在计算流体力学中更重要的变量为密度。基于气体的密度是否为常数,进一步可以分为不可压缩流体和可压缩流体。液体通常考虑为不可压缩流体。 虽然气体比液体的可压缩性高很多,但是是否考虑气体的可压缩性主要取决于具体情况。比如在空气动力学中,如果我们的气流速度很小,气体的可压缩性对飞机的设计影响较小,但是如果飞机的速度和声速接近或者大于声速,气体的可压缩性对飞机的设计非常重要。

    总结:计算流体力学中的流体通常指气体和液体(在某些情况下也可能为固体),气体和液体的区分主要为可压缩性。液体和低速气体通常被认为是不可压缩的,高速气体通常被认为是可压缩的。

    有限体积法中的有限控制体

    CFD的根本就是使用计算机对流体的运动进行预测来指导工业生产。求解CFD的方法有很多种,本文仅仅讨论比较流行的有限体积法。

    如果要分析一大片河流的流动状态、整个房间的热流密度、血管内的血流状态,CFD是通过从局部到整体的方法进行分析。举例:在有限体积法中,其采用的方法是分析局部(非常小一块区域中)的流体的运动来获得整片河流的流动。在这里首先需要介绍流动模型的概念。流通模型通常分为两大类(四种):

    无穷小微团
    图1. 有限控制体和无穷小微团示意图

    以空间位置固定的有限控制体为例,我们要了解一个房间内流体的流动状况。如果房间内每一点的流动均已知,那么我们就可以获得这个房间完整的流体流动图像。举例:我们把房间分为9个位置并用9个坐标来表示,如果我们采用实验的方法对这9个坐标检测速度,那么我们就可以认为这9个速度的场即为整个房间的速度场。进一步的,如果我们把这个房间分为$N$块,并用实验的方法检测这$N$个位置的速度,我们会得到一个更精确的速度场。实际上,房间内的$N$块,即可以近似的理解为CFD中空间位置固定的有限控制体的概念。

    无穷小微团
    图2. 一个网格单元示意图

    图2中可以理解为CFD中的一个网格单元,其为边长为$\delta x$,$\delta y$,$\delta z$的空间位置固定的无穷小微团或有限控制体。在这个空间位置固定的无穷小微团中,其体心坐标为$(x,y,z)$。其具有6个面($N$,$S$等)。流体的相关物理量均为空间位置固定的无穷小微团的函数,如压力表示为$p(x,y,z,t)$,速度表示为$u(x,y,z,t)$。注意的是:如果考虑此网格单元为无穷小微团,相关推导出来的方程为守恒形式(因为位置固定)的微分形式(因为是无穷小微团),如果考虑此网格单元为有限控制体,相关推导出来的方程为守恒形式(因为位置固定)的积分形式(因为是有限控制体)。

    总结:推导CFD方程之前需要定义流动模型,其主要分为空间位置固定的有限控制体和随流体流动的有限控制体,空间位置固定的无穷小微团和随流体流动的无穷小微团。不同的流动模型推导出来的方程形式不同,但可以彼此转化。CFD中的网格可以被看做是位置固定的有限控制体。

    本节内容主要引自Computational fluid dynamics basics with applications,可阅读原书进行拓展理解。

    连续性方程

    本节以连续性方程为例子,一步一步进行推导,旨在让学员了解CFD中基本控制方程的来源,让学员在CFD模拟的时候有根有据,有理论基础。复杂的方程离散、对流格式等将在OKSS2课程中讲述。

    CFD的控制方程无论具有什么形式,都是建立在流体力学基本控制方程——连续性方程、动量方程、能量方程的基础之上。他们主要来源于著名的质量守恒定律、牛顿第二定律、能量守恒定律。CFD中最基本的假设为:把流体看做是连续介质。当我们从宏观尺度来看的时候(比如一微米),流体的分子动力学可以忽略$^2$。这样,我们可以描述流体的一些宏观物理量如速度、压力、密度和温度等。因为这些物理量可以看做大量分子的平均而不具有分子级别的运动。

    回到我们图2所示的网格单元,从这个网格单元出发,可以推导出著名的连续性方程。上一节中我们提到过,其可以看做为空间位置固定的有限控制体或者无穷小微团,我们首先以空间位置固定的无穷小微团为例,推导连续性方程。采用这种观点推导的连续性方程为守恒的微分形式。质量守恒意味着:

    无穷小微团质量的变化 = 流入无穷小微团质量$-$流出无穷小微团质量

    也即为:

    无穷小微团质量的变化率 = $($流入无穷小微团质量$-$流出无穷小微团质量$)$的变化率

    首先,我们用$m$表示无穷小微团的质量,其为关于空间位置$\bfx$和时间$t$的函数(东岳流体所有矢量用直体加粗表示如速度$\bfU$),因此在数学上可以表示为$m(\bfx,t)$。那么其关于时间$t$的变化率(单位:kg/s)可以表示为: \begin{equation} \frac{\partial m(\bfx,t)}{\partial t} \end{equation} 其中我们在这里需要使用偏导数$\frac{\p}{\p}$,这是因为在这里的质量为空间和时间的函数。如果我们的质量仅仅为时间的函数即$m(t)$,那么则可以用: \begin{equation} \frac{\rd m(t)}{\rd t} \end{equation} 来表示质量的变化率。CFD中通常所有变量均为关于时间和位置的函数,有的情况下还会涉及到更多的变量,因此在表示某某变量仅仅关于某某变量的变化率的时候,需要使用偏导数。同时,在下文中我们还会了解到全导数的概念。学员一定要透彻的了解偏导数的概念,在OKS课程中,会出现大量的偏导数。

    进一步的,我们的质量$m$可以表示为密度$\rho$和体积$V$的乘积。再次需要提醒的是,这里的密度$\rho$也是关于时间和位置的函数即$\rho(\bfx,t)$,体积$V$则为网格单元的体积,在本例中是固定的。进阶概念:在动网格中,网格单元的体积具有速度且和空间位置有关。 若我们假定此网格单元为正方体如图3所示,我们的体积$V$可以表示为三个边长的乘积。如我们假定此网格单元为无穷小微团,则三个边长可以用数学符号$\delta x$,$\delta y$,$\delta z$来表示(也可以表示为$\rd x$,$\rd y$,$\rd z$)。因此,我们的质量对于时间的变化率可以表示为: \begin{equation}\label{mass11} \frac{\partial m}{\partial t}=\frac{\partial (\rho V)}{\partial t}=\frac{\partial (\rho \delta x \delta y \delta z)}{\partial t}=\frac{\partial \rho}{\partial t}\delta x \delta y \delta z \end{equation} 在方程$\eqref{mass11}$中,因为$\delta x$,$\delta y$,$\delta z$和时间无关,因此可以提出。
    有限控制体
    图3. 流体中的相关量在无穷小微团中的表示


    现在看图3中的无穷小微团,定义微团中心处的位置矢量为$(0,0,0)$,即位置矢量$\bfx$的各个分量$x=y=z=0$,此点的速度分量分别为$u$,$v$,$w$。在这里需要简单介绍一下泰勒公式,在泰勒公式中,如果对于一个函数,知道$x_0$点的值为$f(x_0)$,那么$x_0+h$的值可以通过下述公式求得: \begin{equation} f(x_0+h)\approx f(x_0)+f'(x_0)h \end{equation} 同样的,对于图三中的微团中心处的速度为$\bfU$(对应分量为$u$,$v$,$w$),那么我们在单位时间内、单位面积上流出$x$方向的$E$面的质量为: \begin{equation} \rho u(0+\frac{1}{2}\delta x) \approx \rho u(0) + \frac{\p \rho u}{\p x}|_{x=0}\frac{1}{2}\delta x = \rho u + \frac{1}{2}\frac{\p \rho u}{\p x}\delta x \label{Emass} \end{equation} 同理,单位时间内、单位面积上流入$W$面的$x$方向的质量为: \begin{equation} \rho u(0) - \frac{\p \rho u}{\p x}\frac{1}{2}|_{x=0}\delta x = \rho u - \frac{1}{2}\frac{\p \rho u}{\p x}\delta x \label{Wmass} \end{equation} 需要注意的是在讨论流入流出的质量的时候,需要附加单位时间内,但是在讨论速度的时候,由于其本身就表示单位时间内移动的距离,因此不需要附加“单位时间内”来进行描述。

    这样,$x$方向上单位时间内质量的变化(质量变化率)就可以通过单位面积上流进来的质量$-$流出去的质量最后乘以面积来求得。将方程\eqref{Wmass}减去方程\eqref{Emass}并乘以$x$方向的面积有: \begin{equation} -\frac{\p \rho u}{\p t}\delta x \delta y \delta z \label{xMass} \end{equation} 其表示从$x$方向上计算的质量变化率。 同理,从$y$方向上计算的质量变化率为: \begin{equation} -\frac{\p \rho v}{\p t}\delta x \delta y \delta z \label{yMass} \end{equation} 从$z$方向上计算的质量变化率为: \begin{equation} -\frac{\p \rho w}{\p t}\delta x \delta y \delta z \label{zMass} \end{equation} 将\eqref{xMass},\eqref{yMass},\eqref{zMass}相加,则表示微团内的总质量变化率: \begin{equation} -\left( \frac{\partial \rho u}{\partial x} + \frac{\partial \rho v}{\partial y} + \frac{\partial \rho w}{\partial z} \right) \delta x \delta y \delta z \label{flux} \end{equation} 结合方程($\ref{flux}$)和方程($\ref{mass11}$)有: \begin{equation}\label{cont} -\left( \frac{\partial \rho u}{\partial x} + \frac{\partial \rho v}{\partial y} + \frac{\partial \rho w}{\partial z} \right) \delta x \delta y \delta z = \frac{\partial \rho}{\partial t}\delta x \delta y \delta z \end{equation} 也即: \begin{equation}\label{contF} \frac{\partial \rho}{\partial t} + \frac{\partial \rho u}{\partial x} + \frac{\partial \rho v}{\partial y} + \frac{\partial \rho w}{\partial z}=0 \end{equation} 方程($\ref{contF}$)即为连续性方程

    在这里或许可以思考一下,连续性方程是几阶的?

    本节讨论笛卡尔坐标系下的张量操作。连续性方程\eqref{contF}在各种教材中具有不同的表现形式,但内容相同。例如方程\eqref{contF}和下面不同的数学表示是相同的: \begin{equation}\label{nabla} \frac{\partial \rho}{\partial t} + \nabla \cdot \left( \rho \mathbf{U} \right)=0 \end{equation} \begin{equation}\label{div} \frac{\partial \rho}{\partial t} + \mathrm{div} \left( \rho \mathbf{U} \right)=0 \end{equation} \begin{equation}\label{ein} \frac{\partial \rho}{\partial t} + \frac{\partial \rho u_i}{\partial x_i}=0 \end{equation} 在这里首先我们看方程\eqref{nabla}中的$\nabla\cdot$操作符,其数学意义表示如下: \begin{equation} \nabla\cdot\mathbf{U}=\frac{\p u}{\p x}+\frac{\p v}{\p y}+\frac{\p w}{\p z} \end{equation} 可见,对一个矢量做$\nabla\cdot$运算,将变成一个标量。$\nabla\cdot\bfU$也通常被称之为求速度的散度。

    同时,这个倒三角还可以用来表示梯度,其和散度的区别在于没有后面的小圆点,求一个标量,如CFD中的压力$p$的梯度如: \begin{equation} \nabla p = \left[\begin{matrix} \frac{\partial p}{\partial x} \\ \frac{\partial p}{\partial y} \\ \frac{\partial p}{\partial z} \end{matrix} \right] \end{equation} 可见,对一个标量做$\nabla$运算,将变成一个矢量。对于压力$p$,我们也可以先求梯度,再求散度: \begin{equation} \nabla \cdot(\nabla p)=\frac{\partial^2p}{\partial x^2}+\frac{\partial^2p}{\partial y^2}+\frac{\partial^2p}{\partial z^2} \end{equation} 其中$\nabla\cdot\nabla$也被称之为拉普拉斯算符,其也可以写成为$\nabla ^2$。在讨论完矢量的散度(变为标量)和标量的梯度(变为矢量)之后,我们还可以对矢量进行梯度操作如: \begin{equation} \nabla \mathbf{U} = \left[ \begin{matrix} \frac{\partial u}{\partial x} & \frac{\partial v}{\partial x} & \frac{\partial w}{\partial x}\\ \frac{\partial u}{\partial y} & \frac{\partial v}{\partial y} & \frac{\partial w}{\partial y} \\ \frac{\partial u}{\partial z} & \frac{\partial v}{\partial z} & \frac{\partial w}{\partial z}\\ \end{matrix} \right] \end{equation} 可见,在对矢量做梯度之后,变为一个二阶张量。如果我们对矢量做梯度之后再做散度,应该是标量还是矢量还是二阶张量?我们可以用一个例子表示:对速度先求梯度再求散度: \begin{equation} \nabla \cdot(\nabla \mathbf{U})= \left[ \begin{matrix} \frac{\partial}{\partial x}\left(\frac{\partial u}{\partial x}\right)+\frac{\partial}{\partial y}\left(\frac{\partial u}{\partial y}\right)+\frac{\partial}{\partial z}\left(\frac{\partial u}{\partial z}\right)\\ \frac{\partial}{\partial x}\left(\frac{\partial v}{\partial x}\right)+\frac{\partial}{\partial y}\left(\frac{\partial v}{\partial y}\right)+\frac{\partial}{\partial z}\left(\frac{\partial v}{\partial z}\right)\\ \frac{\partial}{\partial x}\left(\frac{\partial w}{\partial x}\right)+\frac{\partial}{\partial y}\left(\frac{\partial w}{\partial y}\right)+\frac{\partial}{\partial z}\left(\frac{\partial w}{\partial z}\right)\\ \end{matrix} \right] \end{equation} 以上主要是这个倒三角符号的用法。在这里我们在谈一下$\mathrm{div}$,很简单,$\mathrm{div}$操作即为$\nabla\cdot$,同理$\mathrm{grad}$操作即为$\nabla$。方程\eqref{ein}则在后续中介绍。

    除了对速度和压力进行各种操作外,本节引入二阶张量的概念。在CFD中,最重要的二阶张量主要为应力张量以及形变率。

    粘性应力张量(Viscous Stress)

    单位面积上的应力张量我们通常用字符$\mathbf{\tau}$来表示。带后缀的$\tau_{ij}$表示应力张量的某个分量。第一个下标$i$表示此应力分量作用于和$i$方向垂直的平面,第二个下标$j$表示应力分量的方向。例如,$\tau_{xy}$表示作用于与$x$轴垂直的平面,$y$方向的应力分量。$\mathbf{\tau}$具有9个分量,如下图所示:
    应力张量分量
    图4. 应力张量分量示意图

    应力张量在数学上通常这样表示: \begin{equation}\label{tau} \mathbf{\tau} = \left[ \begin{matrix} \tau_{xx}& \tau_{xy} & \tau_{xz}\\ \tau_{yx}& \tau_{yy} & \tau_{yz}\\ \tau_{zx}& \tau_{zy} & \tau_{zz} \end{matrix} \right] \end{equation} 通常情况下,所有的气体以及大多数液体都是各项同性的。部分高分子液体可能呈现的为各向异性应力张量。如果我们考虑应力张量为各项同性的。则有:$\tau_{xy}=\tau_{yx}$,$\tau_{xz}=\tau_{zx}$,$\tau_{zy}=\tau_{yz}$。也即,我们的应力张量为对称二阶张量

    我们更进一步的可以把应力张量的分量分为正应力切应力。切应力使得有限控制体形状发生改变,但是体积不变。正应力使得有限控制体体积发生变化。图5形象的表示了2维情况下正应力和切应力对有限控制体作用引起的变化。
    应力张量分量
    图5. 正应力和切应力作用于有限控制体示意图

    在图5中,切应力$\tau_{yx}$使得控制体变形,但是体积不变。正应力$\tau_{yy}$使得控制体体积发生变化。在大多数流体中,正应力要比切应力小得多。

    除了正应力外,有限控制体的受力还受压力的作用,即静压。单位面积上的静压其可以这样表示: \begin{equation}\label{sp} -p\bfI=\begin{bmatrix} 1 & 0 & 0\\ 0 & 1 & 0\\ 0 & 0 & 1 \end{bmatrix}=\begin{bmatrix} -p & 0 & 0\\ 0 & -p & 0\\ 0 & 0 & -p \end{bmatrix} \end{equation} 压力和正应力的作用方向正好相反,压力在有限控制体的各个面上对有限控制体进行压缩,因此方程\eqref{sp}的符号为负号,正应力则倾向于使得有限控制体变得膨胀。压力和正应力同样使得有限控制体体积发生变化,但不会发生侧向形变。

    形变率(Strain Rate)

    类似于粘性应力张量,控制体的形变率通常用$\mathbf{S}$来表示。各项同性流体的形变率也为对称二阶张量。其如下表示: \begin{equation}\label{S} \mathbf{S}=\left[ \begin{matrix} S_{xx}& S_{xy} & S_{xz}\\ S_{yx}& S_{yy} & S_{yz}\\ S_{zx}& S_{zy} & S_{zz} \end{matrix} \right] = \left[ \begin{matrix} \frac{\partial u}{\partial x}& \frac{1}{2}\left( \frac{\partial u}{\partial y}+\frac{\partial v}{\partial x} \right) & \frac{1}{2}\left( \frac{\partial u}{\partial z}+\frac{\partial w}{\partial x} \right)\\ \frac{1}{2}\left( \frac{\partial u}{\partial y}+\frac{\partial v}{\partial x} \right) & \frac{\partial u}{\partial y} & \frac{1}{2}\left( \frac{\partial v}{\partial z}+\frac{\partial w}{\partial y} \right) \\ \frac{1}{2}\left( \frac{\partial u}{\partial z}+\frac{\partial w}{\partial x} \right)& \frac{1}{2}\left( \frac{\partial v}{\partial z}+\frac{\partial w}{\partial y} \right) & \frac{\partial w}{\partial z} \end{matrix} \right] \end{equation} 不难发现,方程($\ref{S}$)中定义的形变率也即$S=\frac{1}{2}\left(\nabla \mathbf{U} + \nabla \mathbf{U}^T \right)$。

    牛顿流体

    牛顿流体假设粘性应力和形变率成线性的比例。牛顿流体的粘性应力主要和两个形变率有关,第一个即为上文讨论的形变率$\bfS$,第二个为体形变率$\lambda\nabla\cdot\bfU$。牛顿流体的粘性应力张量和形变率的关系用数学公式可以表示为: \begin{equation} \tau_{xx}=2\mu\frac{\p u}{\p x}+\lambda\nabla\cdot\bfU \end{equation} \begin{equation} \tau_{yy}=2\mu\frac{\p v}{\p y}+\lambda\nabla\cdot\bfU \end{equation} \begin{equation} \tau_{zz}=2\mu\frac{\p w}{\p z}+\lambda\nabla\cdot\bfU \end{equation} \begin{equation} \tau_{xy}=\tau_{yx}=\mu\left(\frac{\p u}{\p y}+\frac{\p v}{\p x}\right) \end{equation} \begin{equation} \tau_{xz}=\tau_{zx}=\mu\left(\frac{\p u}{\p z}+\frac{\p w}{\p x}\right) \end{equation} \begin{equation} \tau_{yz}=\tau_{zy}=\mu\left(\frac{\p v}{\p z}+\frac{\p w}{\p y}\right) \end{equation} 也即为: \begin{equation} \tau=2\mu\bfS+\lambda(\nabla\cdot\bfU)\bfI=\mu \left(\nabla \mathbf{U} + \nabla \mathbf{U}^T \right)+\lambda(\nabla\cdot\bfU)\bfI \end{equation} 其中$\bfI$为单位二阶张量: \begin{equation} \bfI=\begin{bmatrix} 1 & 0 & 0\\ 0 & 1 & 0\\ 0 & 0 & 1 \end{bmatrix} \end{equation} 斯托克斯进一步假设$\lambda=-\frac{2}{3}\mu$,最终有: \begin{equation} \tau=\mu \left(\nabla \mathbf{U} + \nabla \mathbf{U}^T \right)-\frac{2}{3}\mu(\nabla\cdot\bfU)\bfI \end{equation} 对于不可压缩牛顿流体,$\nabla\cdot\bfU=0$,粘性应力的形式更加简单:$\tau=\mu \left(\nabla \mathbf{U} + \nabla \mathbf{U}^T \right)$。

    下面我们看著名的动量方程,在这里对动量方程的推导不做介绍,不可压缩动量方程可以表示为: \begin{equation}\label{mom} \frac{\p \bfU}{\p t}+\nabla\cdot(\bfU\otimes\bfU)=-\nabla p+\nabla\cdot\tau \end{equation} 如何把将方程\eqref{mom}变成分量的形式?

    张量标识法

    本网站速度矢量统一表示为$\bfU$,除了这种表示法外,张量还可以用张量标识法(Tensor Notation)表示。在张量标志法中,符号的下标表示张量的阶数。如标量压力$p$可以记为(二者无区别): \begin{equation} p\equiv p \end{equation} 矢量$\mathbf{U}$可以写为$u_i$: \begin{equation} u_i\equiv\mathbf{U}= \begin{bmatrix} u_1\\ u_2\\ u_3 \end{bmatrix} \end{equation} 二阶张量$\boldsymbol{\tau}$定义为$\tau_{ij}$: \begin{equation} \tau_{ij}\equiv \boldsymbol{\tau}= \begin{bmatrix} \tau_{11} & \tau_{12} & \tau_{13} \\ \tau_{21} & \tau_{22} & \tau_{23} \\ \tau_{31} & \tau_{32} & \tau_{33} \end{bmatrix} \end{equation} 需要注意的是,在使用张量标识法的时候,$u_i$的分量用数字表示如:$u_1$,$u_2$,$u_3$。举例说明:如果某文章采用张量标志法,则$u_i$表示速度矢量,$u_1$表示$x$方向速度分量。对于本文中,$\bfU$表示速度矢量,$u_i$表示各个方向的分量,其中$i=1,2,3$。

    爱因斯坦操作符

    张量标志法通常结合爱因斯坦操作符(Einstein Notation)使用变成非常强大的武器。我们举例说明:当在做乘积操作的时候,如果某个下标重复,那么则需要进行加和。如对于标量$s$: \begin{equation} s=\mathbf{a}\cdot\mathbf{b}\quad\equiv\quad s=a_i b_i=\sum_{i=1}^3 a_ib_i=a_1b_1+a_2b_2+a_3b_3 \end{equation} \begin{equation} s=\boldsymbol{\tau}:\boldsymbol{\tau}\quad\equiv\quad s=\tau_{ij}\tau_{ij}=\sum_{i=1}^3\sum_{j=1}^3\tau_{ij}\tau_{ij}=\tau_{11}\tau_{11}+\tau_{12}\tau_{12}+\tau_{13}\tau_{13}+\tau_{21}\tau_{21}+\tau_{22}\tau_{22}+\tau_{23}\tau_{23}+\tau_{31}\tau_{31}+\tau_{32}\tau_{32}+\tau_{33}\tau_{33} \end{equation} 如对于矢量$c_i$: \begin{equation} \mathbf{c}=\boldsymbol{\tau}\cdot\mathbf{U} \quad\equiv\quad c_i=\tau_{ij}u_j=\sum_{j=1}^3\tau_{ij}u_j=\begin{bmatrix} \tau_{11}u_1+\tau_{12}u_2+\tau_{13}u_3\\ \tau_{21}u_1+\tau_{22}u_2+\tau_{23}u_3\\ \tau_{31}u_1+\tau_{32}u_2+\tau_{33}u_3 \end{bmatrix} \end{equation} 对于二阶张量$\tau_{ij}$: \begin{equation} \boldsymbol{\tau}=\boldsymbol{D}\boldsymbol{E}\quad\equiv\quad \tau_{ij}=D_{ik}E_{kj}=\sum_{k=1}^3D_{ik}E_{kj}=D_{i1}E_{1j}+D_{i2}E_{2j}+D_{i3}E_{3j} \end{equation} \begin{equation} \boldsymbol{\tau}=\boldsymbol{D}^\mathbf{T}\boldsymbol{E}\quad\equiv\quad \tau_{ij}=D_{ki}E_{kj}=\sum_{k=1}^3D_{ki}E_{kj}=D_{1i}E_{1j}+D_{2i}E_{2j}+D_{3i}E_{3j} \end{equation} 克罗内克函数$\delta_{ij}$: \begin{equation} \delta_{ij}=\left\{\begin{matrix} 1 & i=j\\ 0 & i \neq j \end{matrix}\right. \end{equation} 很明显: \begin{equation} \delta_{11}=\delta_{22}=\delta_{33}=1 \end{equation} \begin{equation} \delta_{12}=\delta_{21}=\delta_{13}=\delta_{31}=\delta_{23}=\delta_{32}=0 \end{equation} 同样的,对于偏微分方程组也可以进行类似的表示。例如对于连续性方程($\ref{contF}$): \begin{equation} \frac{\partial \rho}{\partial t} + \frac{\partial \rho u}{\partial x} + \frac{\partial \rho v}{\partial y} + \frac{\partial \rho w}{\partial z}=0 \equiv \frac{\partial \rho}{\partial t} + \frac{\partial \rho u_i}{\partial x_i} =0 \end{equation} 对于动量方程($\ref{mom}$),其可以表示为: \begin{equation} \frac{\partial u_i}{\partial t}+\frac{\partial u_i u_j}{\partial x_j}=-\frac{\partial p}{\partial x_i}+\frac{\partial}{\partial x_j}\left(\frac{\partial u_i}{\partial x_j}\right) \end{equation} 现在随意从文献中抽取一个方程,如《Turbulence Modeling for CFD》中$k-\varepsilon$模型的湍流动能方程,即方程(4.9)为: \begin{equation} \rho\frac{\partial k}{\partial t}+\rho u_j\frac{\partial k}{\partial x_j}=\tau_{ij}\frac{\partial u_i}{\partial x_j}-\rho\varepsilon+\frac{\partial}{\partial x_j}\left[(\mu+\mu_t/\sigma_k)\frac{\partial k}{\partial x_j}\right] \end{equation} 对其分析可知首先其为一个方程而非方程组,且其展开形式为: \begin{equation} \rho\frac{\partial k}{\partial t}+\sum_{j=1}^3\rho u_j\frac{\partial k}{\partial x_j}=\sum_{i=1}^3\sum_{j=1}^3\tau_{ij}\frac{\partial u_i}{\partial x_j}-\rho\varepsilon+\sum_{j=1}^3\frac{\partial}{\partial x_j}\left[(\mu+\mu_t/\sigma_k)\frac{\partial k}{\partial x_j}\right] \end{equation}

    学员可对本网页详细预习后,进行下面的测试验证预习结果。
    学员若对上述材料有不理解的地方,请联系li.dy@dyfluid.com做内容增添。



    北京见 :-)


    东岳流体dyfluid.com
    勘误、讨论、补充内容请前往CFD中国