return 0;
wmake

CFD理论入门

如果打不开图像,请右键在新标签页打开图像后刷新几次


序言和预习资料准备

本页面主要为参加OpenFOAM Knowledge Share(OKS)课程的老师和同学准备的预习资料,为【保证学习效果,所有人员务必对本页内容进行预习】。内容可能看起来较长,但只是CFD的冰山一角。如果对本页内容均已经掌握的老师和同学可略过。【未参加课程的人员也可用来学习,本文算是非常浅显的CFD理论入门资料】。本页面内容将在报名截止前更新完毕。

OKS主要讲授OpenFOAM的使用。OpenFOAM是一个采用C++编写的CFD求解器,其主要的特点在于开源并且免费。开源的特性方便各类研究人员修改代码并植入自己的模型。国内CFD用户目前使用最常用的软件是ANSYS Fluent,ANSYS Fluent在使用的过程中,用户可以通过UDF来自定义函数,但ANSYS Fluent提供的接口有限。OpenFOAM用户则可以任意的修改代码进行任意的函数定义。但是OpenFOAM的学习过程非常缓慢并且资料非常少。在OKS课程中,PCFD主要讲解必要的CFD理论。OKSS1课程主要讲解OpenFOAM的基本使用。在参加OKSS1之前必须参加PCFD,才能保证学习效果。同时,参加OKSS1的人员主要需要预习一些Linux的基本操作。OKSS2课程主要讲解CFD方程在OpenFOAM中的实现以及OpenFOAM代码定制,需要老师和同学具备基本的C++基础。下面是一些相关书籍推荐。

一些人经常问ANSYS Fluent和OpenFOAM的区别。实际上,这个问题可以理解为商软和开源的区别。这个问题,对于精准度、客户支持、模型方面大家都比较了解。但开源软件,最重要的特点是可以近似理解是万能的。商业软件,你无法改动它的代码,因此就不能说是万能的。比如采用C++编写的OpenFOAM,其只不过是庞大的C++库,C++可以做的,OpenFOAM都可以做。如果用户有C++写的LBM代码,完全可以植入到OpenFOAM。如果用户有了一个GUI框架,也可以把OpenFOAM作为后台相结合。有什么是C++编程做不到的呢?除了让它帮你拿拖鞋之外。

高等数学

【第2章】(导数和微分):CFD控制方程即为各种导数和微分的结合;【第3章第3节】(泰勒方程):大量的CFD方程基于泰勒方程导出;【第8章】(向量):了解CFD中的矢量操作;【第9章第7节】中(梯度):了解CFD中的梯度操作;【第11章第6节】:了解有限体积法中的高斯积分以及通量。

C++ Primer Plus

【第7-8章】(函数):学习OpenFOAM中大量的类使用,需要首先了解函数的概念;【第10章】(类):务必要详细学习C++中的类,当初OpenFOAM创始人Henry Weller就是因为不能忍受FORTRAN的特性,将代码全部用C++改写,可以看这个有关Henry Weller的采访;【第11-12章】(类的高级用法):在OpenFOAM中类的使用不仅多,而且还很高级,因此类的高级使用也需要了解;【第13章】(类继承):在OpenFOAM中所有的类均使用了继承用法;再次强调:课堂上没有充足的时间系统讲解C++,参加OKSS2的老师和同学务必在课程之前预习好C++。

OpenFOAM-5.x虚拟机

本节标题指向的链接内含仅供老师和同学内部交流学习使用的、适用于win7和win10的vmware虚拟机以及注册码。首先解压,后安装vmware软件;安装好后,打开Ubuntu文件夹下的dyfluid.vmx即可运行已经安装好的OpenFOAM-5.x。本虚拟机支持中文输入法输入;root密码压缩包内已经提供;现在课程必要的介绍已经完毕。


预习节点

张量标识

学习CFD理论首先遇到的就是各种各样的算符。【能看懂、拆分CFD方程是研究算法的最基本步骤】。实际上,给出一个CFD方程,比如动量方程,最基本最基本的,同学可以不知道是怎么推出来的,但要能看懂。【本节介绍纯粹的张量标识。不涉及到任何的CFD理论】。

本网站【标量采用斜体】,如压力$p$。【矢量采用正体加粗】,如速度矢量$\bfU$,其具有三个分量$u_1$,$u_2$和$u_3$或$u$,$v$,$w$。【二阶张量也采用正体加粗】,比如应力张量$\boldsymbol{\tau}$,其分量为:

\begin{equation} \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}

首先来看下面这个方程(实际上其为粘度为$1$的不可压缩流体动量方程)。再次强调一下,本节不涉及到任何CFD算法,只是介绍如何看懂数学方程,因此方程直接给出:

\begin{equation} \frac{\partial \mathbf{U}}{\partial t}+\nabla \cdot (\mathbf{U}\mathbf{U})=-\nabla \frac{p}{\rho}+\nabla \cdot(\nabla \mathbf{U}) \label{mom} \end{equation}

其中的$\mathbf{U}$为速度矢量,$p$为压力,$\rho$为密度。【这个方程实际上表示3个方程】。一般来讲,CFD文献中通常采用方程\eqref{mom}的形式(紧凑形式)而并不进行展开。下面介绍如何将其展开为3个方程:

结合方程\eqref{time}、\eqref{conv}、\eqref{p}、\eqref{diff}这四项,有三个方程。下面仅列出$x$方向(取$[ ]$中的第一行):

\begin{equation}\label{result} \frac{\partial u_1}{\partial t} + \frac{\partial u_1 u_1}{\partial x}+\frac{\partial u_2 u_1}{\partial y}+\frac{\partial u_3 u_1}{\partial z} = -\frac{1}{\rho}\frac{\partial p}{ \partial x} + \frac{\partial}{\partial x}\left(\frac{\partial u_1}{\partial x}\right)+\frac{\partial}{\partial y}\left(\frac{\partial u_1}{\partial y}\right)+\frac{\partial}{\partial z}\left(\frac{\partial u_1}{\partial z}\right) \end{equation}

方程\eqref{result}即为$u_1$的方程。同样的,各位老师和同学们应该仔细观察方程\eqref{result},其只不过是各种导数的加和。除了$u_1$之外,如果所有的变量均为已知,那么就可以求出$u_1$关于时间的步进(看那个第一项速度关于时间的导数)。很遗憾,方程\eqref{result}中的压力$p$和$u_1$是耦合在一起的(也就是二者相互影响),对这个方程的求解非常复杂,这也是CFD中SIMPLE算法等要处理的问题。本文中,理解到方程\eqref{result}的形式即可,顺便要知道,方程\eqref{result}解析解是否存在都无从考证(改变数学历史的百万美元大奖)。

除了本网站这种表示法外,张量还可以用张量标识法表示,其在SCI论文中也非常常见。在张量标识法中,符号的下标表示张量的阶数。如标量压力$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_i$表示速度矢量,$u_1$表示$x$方向速度分量】。这和不使用张量标志法是不同的,例如本文中,$\bfU$表示速度矢量,【$u_i$表示各个方向的分量而不是速度矢量】,其中$i=1,2,3$。因此,遇到$u_i$,采用不同的数学表示方法,其含义可能不同。

张量标志法通常结合爱因斯坦操作符使用变成非常强大的武器。举例说明:当在做乘积操作的时候,如果【某一项中下标重复,那么则需要进行加和】如:

\begin{equation} a_i b_i=\sum_{i=1}^3 a_ib_i=a_1b_1+a_2b_2+a_3b_3=\mathbf{a}\cdot\mathbf{b} \end{equation}
\begin{equation} \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}=\boldsymbol{\tau}:\boldsymbol{\tau} \end{equation}
\begin{equation} \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}=\boldsymbol{\tau}\cdot\mathbf{U} \label{tautan} \end{equation}

方程\eqref{tautan}中$\tau_{ij}u_j=$的$j$出现两次,需要加和。$i$出现一次,表示矢量。

对于二阶张量,将在一项中最少出现两个下标如$i,j$。例如下文方程中$D_{ik}E_{kj}$中的$i,j$表示二阶张量,$k$表示加和:

\begin{equation} \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}=\boldsymbol{D}\boldsymbol{E} \end{equation}
\begin{equation} \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}=\boldsymbol{D}^\mathbf{T}\boldsymbol{E} \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}

同样的,对于偏微分方程组也可以进行类似的表示。例如对于下面的方程(即连续性方程):

\begin{equation} \frac{\partial \rho}{\partial t} + \frac{\partial \rho u_1}{\partial x} + \frac{\partial \rho u_2}{\partial y} + \frac{\partial \rho u_3}{\partial z}\equiv \frac{\partial \rho}{\partial t} + \frac{\partial \rho u_i}{\partial x_i} =0 \end{equation}

其中$\frac{\partial \rho u_i}{\partial x_i}$中出现两次$i$,则需要加和。对于下面的方程(即动量方程),其可以表示为:

\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(\nu\frac{\partial u_i}{\partial x_j}\right) \end{equation}

其中$\frac{\partial u_i u_j}{\partial x_j}$中出现两次$j$,则需要加和。出现一次$i$,则表示分量。因此其为一个矢量的方程。现在随意从文献中抽取一个方程,如《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}

到此为止,各位老师和同学应该知道对于文献中的一个CFD方程,怎样去做拆分。


本次预习预计耗时30分钟,今天预习内容结束

流体与计算流体力学(CFD)

自本节起,进入到计算流体力学章节。

流体通常泛指气体和液体。不同于固体,流体在施加剪切力的时候,会发生形变。且流体没有固定的形态。流体在快速变形的时候,也会存在一种抵抗的力,这种抵抗的力一旦流动停止,也即消失。流体这种抗拒本身发生变形的力称之为【粘性力】。不同的流体存在不同的粘性力。有些流体的粘性很大,有些流体的粘性很小。大部分日常生活中接触的为【牛顿流体】,其形变率和剪切力呈现线性的关系。部分流体为【非牛顿流体】,其剪切率和形变率呈现非线性的关系。在这里同学们可能不理解剪切力和形变率的概念,没关系,这里只是初步的概念讨论,下文有详述。

液体不同于气体的典型特征为其对体积变化的巨大抵抗力。例如我们不能把一瓶矿泉水压缩成一瓶盖的大小。然而气体对体积变化的抵抗力是相对小的。例如我们可以轻易地用手对一个封闭的空的矿泉水瓶进行压缩。基于流体的密度是否为常数,进一步可以分为【不可压缩流体】和【可压缩流体】。液体通常考虑为不可压缩流体。虽然气体比液体的可压缩性高很多,但是是否考虑气体的可压缩性主要取决于具体情况。比如在空气动力学中,如果我们的气流速度很小,气体的可压缩性对汽车的设计影响较小,但是对于飞机,其速度有可能和声速接近或者大于声速,在这种情况下气体的可压缩性对飞机的设计非常重要。有些高压容器中,即时是液体,也应该考虑可压缩性。因为只有考虑了可压缩性,才能解释某些奇怪的物理现象,比如水锤。

计算流体力学(CFD)是采用计算机求解控制流体流动的偏微分方程组的学科。以日常生活中的电风扇吹风为例。电风扇吹出来的空气的流动看不见、摸不着,但是能明显的感觉出来吹得风不一样。风扇转的越快,风越大越凉快,相对应的能耗也越多。【工业上的设计目标就是用最小的能耗,吹出最大的风让一群人都很凉快】。这对于土豪往往无所谓,土豪通常对任何事情都无所谓。毕竟有钱,随便用电。在这物理现象的背后,风扇吹出来的风的速度、压力、温度等均遵循相关的物理定律,这个物理定律可以写成数学中偏微分方程的形式,其主要有两个方程:一个是【连续性方程】(表明流动着的空气质量是守恒的),另一个是【动量方程】(表明空气的流动满足牛顿第二定律)。暂且可以这样理解:求解这两个关于时间的方程,就可以求出风扇出吹来的风在某一时刻的流速。有了流速之后,同时通过风扇叶片表面的受力计算一下风扇的功率,就可以算出耗电量进行相关的设计。如果通过实验,工程师可以设计100个叶片的形状,检测风力的大小和能耗并做分析,取最优的叶片形状。通过CFD,工程师可以生成100个叶片的几何模型,用计算机求解这100个几何模型,就可以计算出风力的大小和能耗进行分析。在工业应用中,一个机翼的实验工作可能要花费上百万元,一个化工反应器的中试实验工作可能需要上亿。但是通过CFD计算模拟,可能需要几万元就可以了。可见,CFD的目标是让工程师从实验中解脱,并节省实验资金,进而进行工业设计。在有些情况下,可能会听到一些新词汇,如数值风洞、数值水池等,这实际上就是用CFD来模拟风洞内的流场、用CFD来模拟水池内水的流动。

问题在于,用于描述CFD的方程非常难解,甚至目前【解析解】(方程的准确解,比如$x=-1$是$x+1=0$的解析解)的存在与否都无从考证(改变数学历史的百万美元大奖)。目前CFD在学术上的研究,即如何获得这些方程的解,这些解通常被称之为【数值解】。在理想的情况下,数值解无限趋近于解析解。然而,理想很丰满,现实很骨感。用户获得的数值解,总是和解析解存在或多或少的差距。

各种算法犹如各路神仙过河一样,总会通过各种奇葩的方法来获得CFD方程的数值解。比如CFD中有一种方法为光滑粒子法,其将流体看做是粒子,然后对这些粒子进行跟踪进行求解并获得数值解。最近,还有研究学者通过深度学习来假拟CFD的方程并求解。目前CFD求解方法中最广泛的,莫过于【有限体积法】。目前主流的CFD软件,均使用的有限体积法进行求解。

有限体积法的思想也很简单。举例,如果要分析整个河流的流动状态、整个房间的热流密度、血管内的血流状态,有限体积法是通过从局部到整体的方法进行分析。例如我们把整个房间分成一个个的小格子,如果要获得每个小格子内的流动速度,那么这些格子构成的整体就是整个房间的流动速度。格子的数量越多,格子的速度场越能代表整个房间的速度场。本文主要介绍有限体积法。

流动模型

在谈具体的方程之前,首先需要介绍流动模型的概念。有限体积法是通过对每一个具体的流动模型进行分析,进而获得整体解的过程。在这里来理解流动模型可能会觉得比较抽象。但是在阅读完本文之后,回过头再来看流动模型,会有柳暗花明的感觉。

有限控制体模型

有限控制体其进一步可以分为【空间位置固定的有限控制体】,和【随流线运动的有限控制体】。空间位置固定的有限控制体可以理解为一个空旷的房间(当然控制体的体积要小的多),其会流入、流出某些物理量(如质量能量等),这些量的产生和消失符合物理学基本定律。举例,这个房间内部充满了人(不能再容纳更多的人),如果挤进来10个人,那么必然会挤出去10个人。即为【质量守恒】。随流线运动的有限控制体可以理解为随风飘扬的一个气球,气球内部存在大量的空气,且空气和外界不发生交换(质量不变)。这个气球的运动当然满足物理的基本定律。在后续的讨论中,可以看出【从有限控制体推导的方程为积分形式】,且若采用空间位置固定的控制体,方程为【守恒型方程】。若采用随流线运动的控制,方程为【非守恒型】。要透彻的理解CFD中守恒和非守恒的概念,目前来讲对部分老师和同学来说比较不容易。在此只需要知道有不同的方程的概念就可以。

如果打不开图像,请右键在新标签页打开图像后刷新几次
图1. 左图:空间位置固定的有限控制体,流体流入控制体后继续流出。右图:空间位置移动的有限控制体,其沿着流体的流线进行移动。

无穷小微团模型

类似有限控制体,无穷小微团也可以进一步可以分为【空间位置固定的无穷小微团】,和【随流线运动的无穷小微团】。无穷小微团和有限控制体的重要区别在于无穷小微团的体积、质量等变量可以用数学中的$\rd V$、$\rd m$来表示。【有限控制体的体积要比无穷小微团要大】,也即有限控制体中包含了大量的$\rd V$。从数学的角度出发,$\rd V$表示这个无穷小微团的体积足够小。但是需要注意的是,为了使得CFD方程适用于连续介质,这个$\rd V$还要包含一定量的流体分子。正因为无穷小微团的各个属性可以用数学中的 $\rd...$ 来表示,同时有限控制体包含了大量的流体微团,这也就体现了从有限控制体推导出来的方程为积分形式,【从无穷小微团推导出来的方程为微分形式】。

微分形式积分形式的方程,这些概念性的东西,跟具体的CFD求解有什么关系么?

这个问题是个好问题。就像大一的时候,完全不知道学习高等数学有什么用一样。在这里,不要求同学们了解这些概念性的东西。对其有印象就可以了。具体来讲,就是说,目前老师和同学们知道【CFD中的方程有微分形式、有积分形式】俩种形式就行了。在透彻的学习CFD之后回过头来理解,比现在理解要简单的多。就像各位老师和同学现在做CFD,才知道高等数学有什么用一样。

如果打不开图像,请右键在新标签页打开图像后刷新几次
图2. 左图:有限控制体(大圈)和无穷小微团(小正方体)。中图:空间位置固定的无穷小微团,流体流入微团后继续流出。右图:空间位置移动的无穷小微团,其沿着流体的流线进行移动。

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

CFD要做的,就是通过计算机,求解这些划分的小块上的物理量。


本次预习预计耗时20分钟,今天预习内容结束

物质导数

在CFD中,第一个莫名其妙的算符应该就是物质导数。然而物质导数其实非常好理解,如果较为熟悉高等数学,简单来讲,【物质导数就是高数中的全导数】。

透彻的理解物质导数需要借助上文中提及的流动模型。现选取随流线运动的无穷小微团上的温度$T(x,y,z,t)$进行分析。需要注意的是,在上文中指出无穷小微团上的变量可以微分算符$\rd$来表示,这里的温度并没有表示为$\rd T$,其原因在于温度并不是一个随物质的量变化而变化的量。不管存在多少的分子,温度均为$T$。类似的这种量,在CFD被称之为【intensive property】。反之,质量体积等为【extensive property】。

在$t_1$时间点,流体微团的位置在$(x_1,y_1,z_1)$,温度为$T_1(x_1,y_1,z_1,t_1)$。随着微团的移动,在$t_2$时刻,位置变为$(x_2,y_2,z_2)$,同时温度为$T_2(x_2,y_2,z_2,t_2)$。既然温度为时间和位置的函数,则可以对温度进行泰勒展开并省略高阶量(如果对泰勒方程不熟悉,请翻看一下上文中提及的高等数学第3章第3节):

\begin{equation}\label{T} T_2 = T_1 + \left.\frac{\partial T}{\partial x}\right|_{x=x_1}(x_2-x_1) +\left.\frac{\partial T}{\partial y}\right|_{y=y_1}(y_2-y_1) +\left.\frac{\partial T}{\partial z}\right|_{z=z_1}(z_2-z_1) +\left.\frac{\partial T}{\partial t}\right|_{t=t_1}(t_2-t_1) \end{equation}

对方程$\eqref{T}$左右两边除以$t_2-t_1$有:

\begin{equation}\label{T2} \frac{T_2-T_1}{t_2-t_1} =\left.\frac{\partial T}{\partial x}\right|_{x=x_1}\frac{x_2-x_1}{t_2-t_1} +\left.\frac{\partial T}{\partial y}\right|_{y=y_1}\frac{y_2-y_1}{t_2-t_1} +\left.\frac{\partial T}{\partial z}\right|_{z=z_1}\frac{z_2-z_1}{t_2-t_1} +\left.\frac{\partial T}{\partial t}\right|_{t=t_1} \end{equation}

接下来定义【温度$T$在$(x_1,y_1,z_1)$点的物质导数为移动的无穷小微团通过$(x_1,y_1,z_1)$点的时候,无穷小微团针对时间的温度瞬时变化率】,即:

\begin{equation}\label{D} \left.\frac{\mathrm{D} T}{\mathrm{D} t}\right|_{x=x_1,y=y_1,z=z_1}=\lim_{t_2\rightarrow t_1}\frac{T_2 - T_1}{t_2-t_1} \end{equation}

同时,依据速度的定义,有$\lim_{t_2\rightarrow t_1}\frac{x_2-x_1}{t_2-t_1}=u$,$\lim_{t_2\rightarrow t_1}\frac{y_2-y_1}{t_2-t_1}=v$,$\lim_{t_2\rightarrow t_1}\frac{z_2-z_1}{t_2-t_1}=w$,代入到方程$\eqref{T2}$有$(x_1,y_1,z_1)$点的的物质导数:

\begin{equation}\label{Ddef0} \left.\frac{\mathrm{D} T}{\mathrm{D} t}\right|_{x=x_1,y=y_1,z=z_1} =\left.\frac{\partial T}{\partial x}\right|_{x=x_1}u +\left.\frac{\partial T}{\partial y}\right|_{y=y_1}v +\left.\frac{\partial T}{\partial z}\right|_{z=z_1}w +\left.\frac{\partial T}{\partial t}\right|_{t=t_1} \end{equation}

对于任意一点,有

\begin{equation}\label{Ddef} \frac{\mathrm{D} T}{\mathrm{D} t} =\frac{\partial T}{\partial x}u +\frac{\partial T}{\partial y}v +\frac{\partial T}{\partial z}w +\frac{\partial T}{\partial t} \end{equation}

方程$\eqref{Ddef}$即为【温度$T$在笛卡尔坐标系下的物质导数】定义。可见,【物质导数即为数学上物理量对时间的全导数】。参考内积的数学定义(这个链接里面可用于查询一些数学公式),方程$\eqref{Ddef}$可写为:

\begin{equation}\label{Ddef2} \frac{\mathrm{D} T}{\mathrm{D} t}= \frac{\partial T}{\partial t} +\mathbf{U} \cdot \nabla T \end{equation}

上述以温度$T$举例说明物质导数的由来。类似的,有速度$\bfU$的物质导数:

\begin{equation}\label{DUdef2} \frac{\mathrm{D} \bfU}{\mathrm{D} t}= \frac{\partial \bfU}{\partial t} +\mathbf{U} \cdot \nabla \bfU \end{equation}

【速度的物质导数表示流体微团的加速度】。除此之外,任何的流场数据都可以用物质导数来表示,如:

\begin{equation}\label{Ddef3} \frac{\mathrm{D}()}{\mathrm{D} t}=\frac{\partial ()}{\partial t}+u\frac{\partial ()}{\partial x}+v\frac{\partial ()}{\partial y}+w\frac{\partial ()}{\partial z}=\frac{\partial ()}{\partial t}+\mathbf{U}\cdot\nabla() \end{equation}

方程$\eqref{Ddef3}$中的$\frac{\partial ()}{\partial t}$又称为【局部导数】(local derivative),第二项$\mathbf{U}\cdot\nabla()$又称之为【对流导数】(convective derivative)。需要注意的是,如果求解器不调用拉格朗日粒子,方程\eqref{Ddef2}左边的物质导数项并不能够直接求出,其只能通过方程右边的方程计算而来。实际上,方程\eqref{Ddef2}可以理解为:

\begin{equation}\label{lagrangian} \underset{\mathrm{Lagrangian}}{\underbrace{\frac{\mathrm{D}}{\mathrm{D}t}}}=\underset{\mathrm{Euler}}{\underbrace{\frac{\partial }{\partial t}+\mathbf{U} \cdot \nabla }} \end{equation}

其中$\frac{\partial }{\partial t}$表示对时间的变化率,$\mathbf{U} \cdot \nabla $表示对流的变化率。


本次预习预计耗时30分钟,今天预习内容结束

连续性方程

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

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

从空间位置固定的无穷小微团推导连续性方程

连续性方程可以在上文中的4种流动模型基础上进行推导而来,本节从空间位置固定的无穷小微团角度出发,对连续性方程进行推导。这种推导形式是最容易理解的,只涉及到高等数学中的泰勒方程。从空间位置固定的角度来分析的方法,通常被称之为【欧拉观点】。【欧拉观点下的质量守恒意味着位置固定无穷小微团质量的变化$=$流入无穷小微团质量$-$流出无穷小微团质量】。也即为:位置固定无穷小微团质量的变化率$=($流入无穷小微团质量$-$流出无穷小微团质量$)$的变化率。图3为一个质量守恒的实例。

如果打不开图像,请右键在新标签页打开图像后刷新几次
图3. 质量守恒实例。在$t$时刻,有4单位的流体要流入无穷小微团(左图)。在$t+\rd t$时刻,有3单位的流体流出了无穷小微团(右图)。那么这个无穷小微团的变化,即为增加了$1$单位的流体(右图无穷小内增加$1$单位流体)。

下面首先看无穷小微团的质量的变化率。如上文所述,流体微团的体积可以表示为$\rd V=\rd x\rd y\rd z$。同时,有密度$\rho$,则对应的无穷小微团质量为

\begin{equation} \rd m=\rho\rd x\rd y\rd z \end{equation}

依据质量变化率的定义(质量的变化除以时间的变化),其可以表示为

\begin{equation} \frac{\p \rd m}{\p t}=\frac{\p \rho\rd x\rd y\rd z}{\p t}=\frac{\p \rho}{\p t}\rd x\rd y\rd z \label{left1} \end{equation}

下面考虑无穷小微团的流入和流出。在这里需要简单介绍一下泰勒方程,在泰勒方程中,如果对于一个函数,知道$x_0$点的值为$f(x_0)$,那么$x_0+h$的值可以通过下述方程求得:

\begin{equation} f(x_0+h)\approx f(x_0)+f'(x_0)h \end{equation}
有限控制体
图4. 无穷小微团$x$方向的通量。

现在看图4中的无穷小微团,定义立方体左侧【单位面积】的质量流量为$\rho u$,其表示【单位时间】内流入【单位面积】的质量(注意其中的单位时间和单位面积)。同时立方体左侧的面积为$\rd y\rd z$,因此有【单位时间】内流入的质量为

\begin{equation} \rho u\rd y \rd z \label{left} \end{equation}

对于立方体右侧的面,在已知立方体左侧面$\rho u$的定义的时候,右侧可以通过泰勒方程求出,即【单位时间】内流出立方体右侧的【单位面积】的质量为$\rho u+\frac{\p \rho u}{\p x}\rd x$。同样的,立方体右侧的面的面积为$\rd y\rd z$,因此有【单位时间】内流出的质量为

\begin{equation} \left(\rho u+\frac{\p \rho u}{\p x}\rd x\right)\rd y \rd z \label{right} \end{equation}

结合无穷小微团单位时间的流出质量,即方程\eqref{right},和无穷小微团单位时间的流入质量,即方程\eqref{left},有单位时间的$x$方向的净质量变化率为

\begin{equation} \rho u\rd y \rd z-\left(\rho u+\frac{\p \rho u}{\p x}\rd x\right)\rd y \rd z=-\frac{\p \rho u}{\p x}\rd x\rd y \rd z \label{leftRight} \end{equation}

同理,有$y$方向的净质量变化率为

\begin{equation} -\frac{\p \rho v}{\p y}\rd x\rd y \rd z \label{y} \end{equation}

$z$方向的净质量变化率为

\begin{equation} -\frac{\p \rho w}{\p z}\rd x\rd y \rd z \label{z} \end{equation}

将各个方向的质量变化率加和有

\begin{equation} -\frac{\p \rho u}{\p x}\rd x\rd y \rd z-\frac{\p \rho v}{\p y}\rd x\rd y \rd z-\frac{\p \rho w}{\p z}\rd x\rd y \rd z=-\left(\frac{\p \rho u}{\p x}+\frac{\p \rho v}{\p y}+\frac{\p \rho w}{\p z}\right)\rd x\rd y \rd z \label{sum} \end{equation}

同时,结合方程\eqref{left1},有

\begin{equation} \frac{\p \rho}{\p t}\rd x\rd y\rd z=-\left(\frac{\p \rho u}{\p x}+\frac{\p \rho v}{\p y}+\frac{\p \rho w}{\p z}\right)\rd x\rd y \rd z \label{cont0} \end{equation}

也即为连续性方程:

\begin{equation} \frac{\p \rho}{\p t}+\frac{\p \rho u}{\p x}+\frac{\p \rho v}{\p y}+\frac{\p \rho w}{\p z}=0 \label{contF} \end{equation}

方程\eqref{contF}在CFD中至关重要,其直接关系到压力泊松方程的导出。在这里,同学们应该尝试去多次的理解连续性方程推导的本源,即连续性方程实际上表示的是质量守恒。具体的,如何从连续性方程至压力方程、其对CFD的算法有什么稳定性作用,将在PCFD中做讲授。

重申,各位老师和同学应该多次的理解本节内容,【理解连续性方程推导过程的物理意义】。


本次预习预计耗时40分钟,今天预习内容结束

空间位置固定有限控制体推导的连续性方程

上文中从空间位置固定的无穷小微团推导的连续性方程是最容易理解的。下面考虑从空间位置固定有限控制体推导的连续性方程。本节涉及到高等数学的散度定律。若难以理解,可跳过6.2节。本节只不过是另一种推导连续性方程的途径。

如果打不开图像,请右键在新标签页打开图像后刷新几次
图5. 有限控制体上的速度和面矢量。

要提前强调的是,有限控制体和无穷小微团的关系可以理解为有限控制体包含着很多的无穷小微团(如图2左图所示)。因此,有限控制体中的量都需要用积分来进行计算。对于空间位置固定的有限控制体,其质量可通过下面的方程进行计算:

\begin{equation} \int \rho \rd V \label{mass} \end{equation}

其实际上就是高等数学中,质量的计算方程。其对应的质量变化率为

\begin{equation} \frac{\p }{\p t}\int\rho \rd V \label{massD} \end{equation}

下面考虑空间位置固定的有限控制体由于流入流出引起的质量变化。图5中的通量可以定义为

\begin{equation} \rho\bfU\cdot\rd\bfS \label{flux} \end{equation}

那么整个有限控制体上的净通量可以表示为

\begin{equation} -\int\rho\bfU\cdot\rd\bfS=-\int\nabla\cdot\rho\bfU\rd V \label{flux2} \end{equation}

其中负号的引入是因为考虑的为流入减去流出的通量(而不是流出减去流入)。结合方程\eqref{flux2}和方程\eqref{massD},有

\begin{equation} \frac{\p }{\p t}\int\rho \rd V=-\int\nabla\cdot\rho\bfU\rd V \label{flux3} \end{equation}

也即

\begin{equation} \frac{\p \rho}{\p t}+\nabla\cdot\rho\bfU=0 \label{flux4} \end{equation}

方程\eqref{flux4}即为【通过空间位置固定的有限控制体推导的连续性方程】,其和\eqref{contF}是等价的。

从空间位置移动的无穷小微团推导的连续性方程

上文是从空间位置固定的有限控制体以及无穷小微团下推导的连续性方程的。下面讨论另一种推导形式:从空间位置移动的无穷小微团推导的连续性方程。这种方式偏向于拉格朗日思想。若难以理解,可跳过6.3节。本节只不过是另一种推导连续性方程的途径。

在拉格朗日框架下,随流线运动无穷小微团质量的变化$=0$。参考前文物质导数的定义,即

\begin{equation}\label{lagrangian2} \frac{\mathrm{D} }{\mathrm{D} t} \rd m= 0 \end{equation}

其中采用$\rd m$是因为跟踪的为无穷小微团。依据物质导数的定义,方程\eqref{lagrangian2}可以化简为

\begin{equation}\label{lagrangian3} \frac{\mathrm{D} }{\mathrm{D} t} \rd m= \frac{\mathrm{D} \left(\rho\rd V\right)}{\mathrm{D} t}=\rho\frac{\mathrm{D} \left(\rd V\right)}{\mathrm{D} t}+\rd V\frac{\mathrm{D} \rho}{\mathrm{D} t}=0 \end{equation}

\begin{equation}\label{lagrangian4} \frac{\rho}{\rd V}\frac{\mathrm{D} \left(\rd V\right)}{\mathrm{D} t}+\frac{\mathrm{D} \rho}{\mathrm{D} t}=0 \end{equation}

依据物质导数和速度散度的物理意义,$\frac{1}{\rd V}\frac{\mathrm{D} \left(\rd V\right)}{\mathrm{D} t}$实际上即表示流动着的单位体积的流体微团的体积相对于时间的变化率,且有

\begin{equation}\label{lagrangian5} \frac{1}{\rd V}\frac{\mathrm{D} \left(\rd V\right)}{\mathrm{D} t}=\nabla\cdot\bfU \end{equation}

代入到方程\eqref{lagrangian3},有

\begin{equation}\label{lagrangian6} \frac{\mathrm{D} \rho}{\mathrm{D} t}+\rho\nabla\cdot\bfU=0 \end{equation}

6.2和6.3节为采用不同的流动模型推导的连续性方程。各位老师和同学若感觉理解吃力,可略过,但6.1节要掌握。

连续性方程小结

方程\eqref{contF}、\eqref{flux4}以及\eqref{lagrangian6}均表示连续性方程,但各自的推导采用了不同的流动模型。有的模型出发点和推导比较容易理解,有的方法比较复杂。这些方程可以相互的转换。且各种教材中书写并不统一。例如,下面的方程均为采用不同书写方式书写的连续性方程:

\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}

实际上,$\nabla\cdot$和$\mathrm{div}$等价,均表示散度操作。方程\eqref{ein}则是另一种表示方式。初学CFD可能对方程中的这些符号感到很生疏,在此,老师和同学们记住下面的连续性方程形式就可以:

\begin{equation} \frac{\p \rho}{\p t}+\frac{\p \rho u}{\p x}+\frac{\p \rho v}{\p y}+\frac{\p \rho w}{\p z}=0 \label{contF2} \end{equation}

其和上面的三个方程是等价的。在后续的讨论中,老师和同学们需要熟悉各种张量表示形式。


本次预习预计耗时30分钟,今天预习内容结束

动量方程

本节讨论动量方程。动量方程依然可以通过上面四种不同的流动模型进行推导。由于篇幅所限本节不再一一推导,仅介绍一种较为容易理解的推导方法,即从【空间位置移动的无穷小微团进行分析】。在推导之前,需要介绍一些基础的知识。

受力分析

对动量方程进行推导,就需要对无穷小微团进行受力分析。正是因为流体的受力,才会引致流体的流动。无穷小微团的受力可区分为【体积力】和【表面力】。【表面力】为作用在无穷小微团面上的力,如压力、表面张力,海洋表面上的风也可以认为是表面力。【体积力】作用在无穷小微团全部的体积上(不仅表面有,体积内也有),例如重力、电磁力,或者一切引起旋转的力(科氏力)。在这些力中,最重要的表面力为压力和应力(后续会对应力进行介绍),最重要的体积力为重力(若不考虑其他源项力)。

首先看压力,如图6所示,CFD中的压力用$p$来表示(其实际为压强),且这个压力表示静压。其作用主要导致流体的压缩和膨胀,而不是切应力引致的变形。另外,压力为一种【正应力】(不同于下文中要介绍的切应力),其永远【作用于无穷小微团的面的垂直的方向】。

如果打不开图像,请右键在新标签页打开图像后刷新几次
图6. 左图:无穷小微团由于压力的作用导致体积的压缩。右图:三角形网格和四边形网格收到的垂直的压力。

如果打不开图像,请右键在新标签页打开图像后刷新几次
图7. 左图:一块蜂蜜由于剪切应力的作用导致的变形。右图:流体由于上方剪切力作用导致的流速区别,壁面附近速度较小。红色的箭头表示剪切应力,黑色的箭头表示速度大小和方向。

同时,无穷小微团还收到【剪切应力】的作用。其主要体现为流体的粘性。可以很明显的注意到,蜂蜜要比水要粘。这就是因为要获得同样大的形变率,施加给蜂蜜的剪切应力需要更大。【剪切应力的作用导致流体产生形变】。这个应力不是和形变成正比,而是【和形变的速率成正比】。也就是说,【应力越大,形变的速度越快】。

如果打不开图像,请右键在新标签页打开图像后刷新几次
图8. 左图:一个四边形无穷小微团所受的剪切应力。右图:一个四边形有限控制体受到剪切应力导致形变。

如图8左图所示,一个二维四边形无穷小微团所受的剪切应力分量分别为

\begin{equation} \tau_{yy},\tau_{yx},\tau_{xx},\tau_{xy} \label{tau} \end{equation}

其中的$\tau$表示剪切应力,【第一个下标表示作用于与某方向垂直的平面】,【第二个下标表示力的方向】。比如$\tau_{yy}$表示作用于和$y$垂直的平面的剪切应力分量,其方向为$y$方向。$\tau_{yx}$表示作用于和$y$垂直的平面的剪切应力分量,其方向为$x$方向。图8右图表示了流体由于剪切应力作用发生的形变。图9则是三维的无穷小微团受到的剪切应力。在三维的情况下,剪切应力的分量可以表示为

\begin{equation} \tau_{xx},\tau_{xy},\tau_{xz},\tau_{yx},\tau_{yy},\tau_{yz},\tau_{zx},\tau_{zy},\tau_{zz} \label{tau3D} \end{equation}

其通常写成下面的形式:

\begin{equation} \begin{bmatrix} \tau_{xx}& \tau_{xy} & \tau_{xz}\\ \tau_{yx}& \tau_{yy}& \tau_{yz}\\ \tau_{zx}& \tau_{zy} & \tau_{zz} \end{bmatrix} \label{tau3DM} \end{equation}
如果打不开图像,请右键在新标签页打开图像后刷新几次
图9. 三维剪切应力。

在这里,老师和同学们应该知道方程\eqref{tau3DM}中的剪切应力分量的定义是什么?分别作用在哪个面?哪个方向?其会引起什么样的物理过程。初学CFD目前理解到这里就可以。下面来分析一个具体的无穷小微团收到的力。(注:下面的图10看起来是非常复杂的。在找了若干资料之后,均没有发现一个不需要动脑的介绍方式。经过自己组织语言之后,发现只是更加复杂化。在这里只能摘取《计算流体力学基础及其应用》的描述来概括)

如果打不开图像,请右键在新标签页打开图像后刷新几次
图10. $x$方向各个面收到的力。

图10表示施加在无穷小微团$x$方向上的全部表面力。在面abcd上,仅仅存在由切应力引起的$x$方向的分量$\tau_{yx}\rd x\rd z$。面efgh和面abcd的距离为$\rd y$,所以efgh面上的$x$方向的切应力为$(\tau_{yx}+\p \tau_{yx}/\p y \rd y)\rd x\rd z$。对于面abcd和面efgh上的切应力,要注意他们的方向。在底面,$\tau_{yx}$是向左的(与$x$轴方向相反),在顶面,$\tau_{yx}+(\p \tau_{yx}/\p y )\rd y$是向右的(与$x$轴方向相同)。这与下面的约定是一致的:即速度的三个分量的正的增量与坐标轴的正向一致。例如,图10中的平面efgh,因为$u$沿着$y$轴正向是增加的,所以在稍微高于平面efgh的地方,速度$u$要比平面efgh上的$u$大。于是就形成了“拉”的动作,试图将流体微团向右拉向$x$轴的正向。相反,若考虑平面abcd,则在稍稍低于平面abcd的地方,速度$u$要比平面abcd上的$u$小。于是对流体形成了一个阻碍的作用,作用在$x$轴的负向。图10中其他剪切应力的方向,都可以用相同的方式进行判断。特别是在面dcgh上,$\tau_{zx}$指向$x$轴负向,而在面abfe上,$\tau_{zx}+(\p \tau_{zx}/\p z)\rd z$指向$x$轴正向。在垂直于$x$轴的面adhe上,$x$方向的力有压力$p\rd x \rd z$,指向流体微团的内部;还有沿$x$轴负向的应力$\tau_{xx}\rd y\rd z$。依据前面提到的速度增量方向的约定,我们可以解释为什么在图10中,面adhe上$\tau_{xx}$是指向左边的。依据规定,速度$u$的正增量和$x$轴的正向一致,所以稍微离开面adhe左面一点点,$u$的值比面adhe上的u要小。因此,正应力的粘性作用在面adhe上就好像是一个吸力,产生一个向左拉的作用,想要阻止流体微团的流动。与此相反,在面bcgf上,压力$(p+(\p p/\p x)\rd x)\rd y \rd z$指向流体微团的内部(沿着$x$轴负向)。而由于在稍微离开面bcfg右面一点点的地方,$u$的值比面bcfg上的$u$要大,就会产生一个由正应力引起的吸力,将流体微团向右拉,这个力的大小为$(\tau_{xx}+(\p \tau_{xx}/\p x))\rd y\rd z$,方向指向$x$轴正向。

如果上面的描述同学们感觉难以理解,下面的结论要知悉。通过对流体微团$x$方向的受力做分析,有$x$方向所受到的表面力(压力和剪切应力贡献)为:

\begin{equation} F_x=\left(p-\left(p+\frac{\p p}{\p x}\rd x\right)\right)\rd y\rd z + \left(\left(\tau_{xx}+\frac{\p \tau_{xx}}{\p x}\rd x\right) - \tau_{xx}\right)\rd y\rd z + \left(\left(\tau_{yx}+\frac{\p \tau_{yx}}{\p y}\rd y\right) - \tau_{yx}\right)\rd x\rd z + \left(\left(\tau_{zx}+\frac{\p \tau_{zx}}{\p z}\rd z\right) - \tau_{zx}\right)\rd x\rd y \label{x} \end{equation}

\begin{equation} F_x=\left(-\frac{\p p}{\p x}+\frac{\p \tau_{xx}}{\p x}+\frac{\p \tau_{yx}}{\p y}+\frac{\p \tau_{zx}}{\p z}\right)\rd x\rd y \rd z \label{x2} \end{equation}

同理,有$y$方向和$z$方向:

\begin{equation} F_y=\left(-\frac{\p p}{\p y}+\frac{\p \tau_{xy}}{\p x}+\frac{\p \tau_{yy}}{\p y}+\frac{\p \tau_{zy}}{\p z}\right)\rd x\rd y \rd z \label{y2} \end{equation}
\begin{equation} F_z=\left(-\frac{\p p}{\p z}+\frac{\p \tau_{xz}}{\p x}+\frac{\p \tau_{yz}}{\p y}+\frac{\p \tau_{zz}}{\p z}\right)\rd x\rd y \rd z \label{z2} \end{equation}

方程\eqref{x2}、\eqref{y2}和\eqref{z2}即为流体微团在三个方向上的受力(为一个矢量):

\begin{equation} \bfF=\begin{bmatrix} \left(-\frac{\p p}{\p x}+\frac{\p \tau_{xx}}{\p x}+\frac{\p \tau_{yx}}{\p y}+\frac{\p \tau_{zx}}{\p z}\right)\rd x\rd y \rd z\\ \left(-\frac{\p p}{\p y}+\frac{\p \tau_{xy}}{\p x}+\frac{\p \tau_{yy}}{\p y}+\frac{\p \tau_{zy}}{\p z}\right)\rd x\rd y \rd z\\ \left(-\frac{\p p}{\p z}+\frac{\p \tau_{xz}}{\p x}+\frac{\p \tau_{yz}}{\p y}+\frac{\p \tau_{zz}}{\p z}\right)\rd x\rd y \rd z \end{bmatrix} \label{forceEnd} \end{equation}

本次预习预计耗时30分钟,今天预习内容结束

动量守恒

在分析完受力之后,回到牛顿第二定律,即【力等于质量乘加速度】。力已经在方程\eqref{forceEnd}中给出,现在来看质量和加速度。在上文中,已经知道流体微团的质量可以表示为$\rd m=\rho\rd x\rd y\rd z$,同时,流体微团的加速度为$\frac{\mathrm{D}\bfU}{\mathrm{D} t}$(参考物质导数一节),其分量形式为

\begin{equation} \left[\frac{\mathrm{D}u}{\mathrm{D} t},\frac{\mathrm{D}v}{\mathrm{D} t},\frac{\mathrm{D}w}{\mathrm{D} t}\right]^{\mathrm{T}} \label{uc} \end{equation}

那么有$x$方向的动量方程为

\begin{equation} \rho\rd x\rd y\rd z\frac{\mathrm{D}u}{\mathrm{D} t}=\left(-\frac{\p p}{\p x}+\frac{\p \tau_{xx}}{\p x}+\frac{\p \tau_{yx}}{\p y}+\frac{\p \tau_{zx}}{\p z}\right)\rd x\rd y \rd z \label{xx} \end{equation}

同理,$y$方向和$z$方向可以类似导出。即

\begin{equation} \rho\frac{\mathrm{D}u}{\mathrm{D} t}=-\frac{\p p}{\p x}+\frac{\p \tau_{xx}}{\p x}+\frac{\p \tau_{yx}}{\p y}+\frac{\p \tau_{zx}}{\p z} \label{xx2} \end{equation}
\begin{equation} \rho\frac{\mathrm{D}v}{\mathrm{D} t}=-\frac{\p p}{\p y}+\frac{\p \tau_{xy}}{\p x}+\frac{\p \tau_{yy}}{\p y}+\frac{\p \tau_{zy}}{\p z} \label{yy2} \end{equation}
\begin{equation} \rho\frac{\mathrm{D}w}{\mathrm{D} t}=-\frac{\p p}{\p z}+\frac{\p \tau_{xz}}{\p x}+\frac{\p \tau_{yz}}{\p y}+\frac{\p \tau_{zz}}{\p z} \label{zz2} \end{equation}

方程\eqref{xx2}、\eqref{yy2}和\eqref{zz2}即最终的动量方程。现在,老师和同学们应该结合上文介绍的相关知识,将方程\eqref{xx2}、\eqref{yy2}和\eqref{zz2}左边项写成欧拉形式(参考物质导数一节)。然后,再将这三个分量形式的动量方程,写成紧凑形式(结合张量一节)。若有问题,可邮件联系。若没有问题,恭喜,可以继续往下看。

实际上,方程\eqref{xx2}、\eqref{yy2}和\eqref{zz2}还不能求解,因为其中的未知量为方程左边的速度,但方程右边还存在未知量$\tau$,这种【待求未知量之外还存在未知量的情况,被称之为封闭问题】。为了将其封闭,需要将$\tau$表示为$\bfU$的函数。剪切应力$\tau$和速度$\bfU$的关系,被称之为【本构关系】。他们的方程,被称之为【本构方程】。17世纪末牛顿指出,流体的剪切应力和速度的梯度成线性关系。这种流体因此也被称之为【牛顿流体】。还存在一些不符合牛顿流体定义的流体,被定义为【非牛顿流体】,会展现一些不常见的流动特性。

下面讨论一个新的CFD变量:形变率,其也被称为剪切速率。这和上述几个方程的封闭有什么关系呢?这是因为未知量$\tau$可以表示为形变率的函数,形变率进一步可以表示为速度的函数。这样,方程\eqref{xx2}、\eqref{yy2}和\eqref{zz2}就可以封闭了。形变率和剪切应力类似,是一个二阶张量,一般用$\bfS$来表示。但由于在CFD中,形变率本身并不会封闭方程(而是形变率和速度的关系封闭了方程),因此在CFD方程中往往不会出现形变率这一项(注意,$\bfS$这个符号在有限体积法中,更倾向于被定义为面矢量)。若$\bfS$符号表示形变率,即:

\begin{equation} \bfS= \begin{bmatrix} S_{xx} & S_{xy} & S_{xz} \\ S_{yx} & S_{yy} & S_{yz} \\ S_{zx} & S_{zy} & S_{zz} \end{bmatrix} \label{S} \end{equation}

对于【各向同性流体,形变率的9个分量有6个是独立的】:

\begin{equation} S_{xy}=S_{yx},S_{xz}=S_{zx},S_{yz}=S_{zy} \label{Sequal} \end{equation}

当然也存在一些各向异性流体,如一些聚合物。对这些流体的研究会导致问题更加复杂。在以后的学习中,各位老师和同学们也会发现,很多模型都是从各向同性假定出发的。形变率进一步的可以表示为速度的方程:

\begin{equation} \bfS= \begin{bmatrix} \frac{1}{2}\left(\frac{\p u}{\p x}+\frac{\p u}{\p x}\right) & \frac{1}{2}\left(\frac{\p u}{\p y}+\frac{\p v}{\p x}\right) & \frac{1}{2}\left(\frac{\p u}{\p z}+\frac{\p w}{\p x}\right) \\ \frac{1}{2}\left(\frac{\p v}{\p x}+\frac{\p u}{\p y}\right) & \frac{1}{2}\left(\frac{\p v}{\p y}+\frac{\p v}{\p y}\right) & \frac{1}{2}\left(\frac{\p v}{\p z}+\frac{\p w}{\p y}\right) \\ \frac{1}{2}\left(\frac{\p w}{\p x}+\frac{\p u}{\p z}\right) & \frac{1}{2}\left(\frac{\p w}{\p y}+\frac{\p v}{\p z}\right) & \frac{1}{2}\left(\frac{\p w}{\p z}+\frac{\p w}{\p z}\right) \end{bmatrix} \label{Su} \end{equation}

这样,方程\eqref{Su}可以写成$\bfS$和$\bfU$的关系,各位老师和同学应该尝试将方程\eqref{Su}重写为紧凑形式。回到牛顿流体,其定义【剪切应力和形变率为线性的关系】(考虑不可压缩流体、即忽略体膨张率的情况会比较简单):

\begin{equation} \tau=2\mu\bfS \label{Stau} \end{equation}

其中$\mu$是粘度,其为一个物理属性。如果将方程\eqref{Su}代入到方程\eqref{Stau},然后代入到方程\eqref{xx2}、\eqref{yy2}和\eqref{zz2},那么其中的$\tau$就被速度表示出来,进而被封闭了。这个封闭的方程若展开,为4个方程,且存在未知项速度$\bfU$和压力$p$共4个未知量,即可求解。

至此,连续性方程已经完全的推导出来。各位老师和同学应该在方程\eqref{xx2}、\eqref{yy2}和\eqref{zz2}的基础上,通过剪切应力、形变率、速度的关系,自行推导出封闭的动量方程(不包含剪切应力项)。若推导成功,本次预习已经成功了95%。剩下的5%,老师和同学们应该仔细观察这4个方程,将其写成紧凑形式(2个方程),再转换成分量形式。观察紧凑形式下的散度、梯度算法如何和分量形式的方程对应。做到对动量方程不陌生、不恐惧。其只不过是各种偏导数放在了一起。具体如何通过有限体积法离散、各项在数值上有什么意义、耦合算法、数值格式以及算法稳定性和对角占优将在PCFD上介绍。要在一天点通求解N-S方程,还是强调,各位老师和同学应充分的熟悉N-S方程。


预习时间不定,今天预习内容结束

Linux操作

本节很轻松,讨论一下Linux的基本操作。在linux中有众多发行版,OKS课程选择比较接近windows使用习惯的Ubuntu系统。可将Ubuntu系统安装在windows下的虚拟机中,也可以选择windows和Ubuntu双系统共存的方式。在Linux系统中,大部分操作在终端中进行,OpenFOAM的运行也主要在终端中进行操作。下面我们介绍linux的主要基本操作。

返回到用户目录:

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 log 

这样会在终端持续的输出log文件。查找src/目录下header.H文件的路径:

find src/ -name header.H 

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

grep -R "FOAM" 

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

grep -l FOAM log

有的话输出log,没有的话不显示。上述内容是部分常用的Linux基本操作。对于入门的老师和同学们,掌握一些基本操作就可以,若对Linux操作有更深的兴趣,可学习一些功能更强大的一些操作如awk,sed,以及一些命令行下的编辑器如VIM等。

更新历史
最近更新:2018.06.14:将"学员"称呼改为"老师和同学"

东岳流体®版权所有
勘误、讨论、补充内容请前往CFD中文网

');