热物理过程及数值模拟-计算传热学.doc

热物理过程的数值模拟Numerical Simulation of Thermophysics Process讲稿主讲：李隆键第一章 概 论11 流动与传热过程的予测方法及特点流动、传热、燃烧问题是热工类各专业和机械类动力机械专业所研究和解决的主要问题之一，燃烧问题实际上是有化学反应的流动与传热问题，推而广之，在所有热物理过程中，几乎都涉及到流动、传热问题。预测的重要性： 在规定设计参数的相应的结构下，热物理过程是否满足要求，达到预定的指标？要预测； 优化设计，不同方案的比较，要预测； 减少设计、生产、再设计和再生产的费用； 减少设计更改； 减少试验和测量次数。问题的核心：速度场、温度场（传热量）、浓度场等。一、热物理问题的予测方法：理论分析法、实验测定、数值模拟1、理论分析以数学分析为基础，求解描述热物理过程的定解问题，获得函数形式的解，表示求解区域内物理量连续分布的场（速度场、温度场、浓度场）。控制方程+单值条件（数学模型）理论解（分析解，解析解）根据解的准确程度，又可再分为：（1）精确分析解（严格解）特点：函数形式的解；它在求解区域精确地满足定解问题。具体解法：直接积分法、分离变量法、积分变换法、热源法、映射法。（2）近似分析解法特点：函数形式的解，在求解区域上近似地满足定解问题（但在总量上满足相应的守恒原理，动量守恒、动量守恒、能量守恒、质量守恒）。具体解法：积分法（从积分方程出发）变分近似解法摄动法（从微分方程出发）2、实验测定（1）纯实验法（2）相似理论实验法：同类相似，减少变量数目减少工作量，得到规律性结果，可直接应用。（3）实验类比法：异类相似物理现象不同，规律相同：微分方程形式相同，单值性条件类似电热类比，水热类比3、数值模拟以数值计算方法为基础，借助（利用）电子计算机求解物理过程的方法热物理过程的数值模拟，对传热过程称为传热的数值模拟、数值传热、计算传热。如前述，传热过程函盖了流动、燃烧，所以计算传热学实质上就代表了热物理过理过程的数值模拟。用电子计算机对热物理问题进行数值计算就象在实验室中对该现象进行实验测定一样，可称之为“数值实验”。随着高速、大容量电子计算机的发展，特别是微型计算机的普及和推广，这种数值实验的方法越来越被更多的科技人员掌握和应用，成为解决热物理过程的一种重要方法。二、予测方法的比较与选择1、分析解法（1）精确预测了数学模型所控制的的热物理过程；（2）函数形式的解使得可以确定区域中任意位置物理量的大小；（3）以显函数的形式，展示各有关参量对该热物理过程的影响；（4）由于是函数形式的解，便于进一步的运算、处理，例如求导、积分。缺点：（1）获得分析解的可能较小；（2）即使能求得分析解，也常常是无究级数，特殊函数以及涉及特征值问题的超越函数，要得到具体的数值结果，也需要繁复的计算；（3）数学模型的结果也需要有实验检验。2、实验方法（1）可以获得热物理过程可靠的数据资料；（2）全比例设备实验可予测由它完全复制的同类设备在相同条件下将如何运行和变化；（3）是研究一种新的基本现象的唯一方法；（4）是检验其它预测方法准确程度的标准。缺点：（1）全比例实验代价大（投资，物力，人力，周期）；（2）缩小比例模型实验结果的外推受准则数实验范围的限制，有些在全比例设备上才能出现的特征在缩小比例模型上并非总是能模拟（例如流动的涡），降低了模型试验的效果；（3）测试困难及测量误差；（4）有些过程无法预先进行试验（航天，气象预报）。3、数值模拟（1）成本低：在大多数实际应用中，计算机运算的成本要比相应的实验研究成本低好几个数量级，对象愈庞大，过程愈复杂，此优点愈突出；同时，与大多数物品价格不断上涨的趋势相反，计算成本还会降低；（2）速度快，周期短；不同方案的对比计算和优选，这对某些大型实验几乎是不可能的。（3）信息完整：能提供计算区域内所有各个位置上有关变量的值（速度、压力、温度、浓度等），而实验则不可能测出整个区域各点处所有变量的值。（4）具有模拟真实条件的能力：几何条件、边界条件、物性条件、初始条件很容易模拟真实条件，不需要采用缩小模型或冷态实验，无论大小、高位，低温、过程快、慢。（5）具有模拟理想条件的能力：对于研究物理现象而不是工程问题时，注意力集中几个基本参数而要设法消除所有无关的因素。几何条件（维数变化，尺寸）、物性（常密度），BC（绝热表面），ic（特定的初始温度分布）。缺点：（1）数值模拟的对象是数学模型简化处理，结果的准确性有特价检验；（2）对一些十分复杂的问题（几何形状复杂，强烈非线性、物性变化大），数值解可能很难获得，或者即便可以获得，代价也是相当昂贵的，例如，对湍流问题，要想通过求解非稳态N-S方程来算出它们的全部与时间相关的结构，则仍然是计算所不能及的；（3）对解的唯一判断力较弱为了进一步讨论数值模拟的缺点，可以把所有的实际问题分成两大类：A类：有完整数学模型的一类问题，如热传导、层流问题、简单的湍流边界层问题；B类：迄今无完整数学模型的一类问题，如复杂湍流、某些非牛顿流体、某些两相流动等，问题的分类还有一标准问题，即描述到什么样的程度可以认为是“够了”，“合适”的。A类缺点：对这类问题，用计算机求解的优越性远远大于实验研究。 有数值模拟缺点的（2）、（3）两条；在某些情况下也需要进行实验检验。对于此类问题，研究计算方法的目的在于使这些计算方法理加可靠、准确和有效，随着研究的进展，其缺点将被不断克服。B类缺点：A类的缺点B类全有，此外，必须进行实验检验。数学模型的研究不断地把B类问题转化为A类：试算与修正。先提出一个模型计算求解与实验结果进行比较修正模型，并不断完善、湍流模型的最新发展就是这种转换的一个典型例子，k-双方程模型最初建立在科尔莫戈洛夫（Kolmgorov）(1942)及普朗特（Prandtl）(1945)的工作基础上的，但并未，也不可能付诸实现，只有到了20世纪70年代，当计算机和计算方法变得更加强有力的是候，该模型才逐步趋于完整并付诸实际应用。4、方法的选择三种方法或三种手段相辅相成，互为补充。（1）分析解可以为检验数值模拟结果的准确度提供比较依据；常常用有分析解的简单问题检验方法的准确度；（2）简单的解极解可以为发展数值方法中的某些算法提供理论依据，调和平均；（3）物理规律、数学模型的正确建立必须通过对现象的充分观察和测定，等；（4）出现在数学模型中的物性参数只有通过实验测定才能获得。数值模拟的对象是热物理过程的数学模型，所以其结果的准确度首先取决于数学模型反映实际热物理过程的准确度（包括所用的特性参数），然后才是所采用的数值方法，计算机并不能创造信息，发现规律，它只能把人们所送入的信息，按照计算者安排的程式对信息进行加工、处理，从而得到相应的结果。但是，一旦确立了与实际物理过程相符合的物理模型、数学模型、数值模拟又可以发挥很大的作用，它可以减少实验工作量，拓宽实验研究的范围，实现对理想单值性条件的模拟；对那些耗资巨大，条件恶劣的实验，或者难于进行的实验来说，“数值实验”更是一种有吸引力的辅助或替代手段。理论分析，实验测定和数值模拟有机而协调地结合，是研究热物理过程理想而有效的方法。12 本课程的内容及安排一、内容两个组成部分1. 理论部分： 基础理论（数学物理、数值方法）热物理过程的数值模拟：通用性，并以热传导、对流换热、流场的计算为例，推广到通用控制方程所描述的其它现象。论述方式的特点：（1）强调物理的概念和方法，而不过分倚重纯数学的推导。（2）以一维为基础，推广（扩展）到二维、三维。2、实践性环节：一个课程设计，程序设计，视情况而定。二、授课方式改变注入式，实行启发式，培养自学能力。三、教材、参考书1、SV帕坦卡著、张政译（郭宽良译），“传热与流体流动的数值计算”。（Numerical Heat Transfer and Fouid Flow），科学出版社，1984。2，陶文铨编著，数值传热学，西安交通学大出版社，1988。第二章 物理现象的数学描述数值计算的对象过程的数学模型，核心是控制方程。（数学模型：控制方程+单值性条件，？单值性条件）21控制微分方程1、控制微分方程的意义控制微分方程是一定守恒原理的数学表达式，影响因变量的各因素之间必定存在某种联系。回忆导热方程：热力学第一定律，对任意控制容积V，导入控制容积的热流量+控制容积中的内发热量=控制容积中物质能量的增量。导入控制容积的净热流量，控制容积内的发热量=，控制容积中物质内能的增量=。 散度定理： 非稳态项 扩散项 源项比容，对常密度（不可压缩）物质运动（流动）时，增加：物体宏观运动带入的能量= 压缩机械动 = 粘性摩擦功 =流体的变形率张量第二动力粘度，体积粘度对简单可压缩系统= 非稳 对流 扩散对理想气体及恒密度（固液）， ，对固、液，忽略p变化，则各项则代表着各因素在单位容积时的作用效果. 单位容积焓的变化率单位时间、单位面积、传递的焓。对流流量密度JC 单位容积流出的净焓量扩散流量密度 源项2、化学组分方程的质量分量扩散流密：3、能量方程如果为常数，OK时h=0，则由4、动量方程：某方向上的动量变化率，单位体积质量x方向的动量，其它粘性力项，5、湍流的时间平均方程湍流：给定点处的物理量随时间而变化，随机性。工程上：关心的是运动状态的时间平均特性非稳定层流方程的时间平均方程，滞流流动的时间平均方程，假设湍流中存在相对于时均值的脉动，平均化（时均化）运算附加项，雷诺应力（湍流热流），湍流扩散流量密度湍流模型，用平均性质来表示的附加项。引入湍流粘度（湍流扩散系数）滞流应力，流量必度时均方程层流方程，相应的层流交换系数用有效系数代替。模型6、通用方程（1）向量形式对于不同的，有特定的、S与之间相对应。大多数的传递过程，扩散流量密度Jd因变量的梯度确定，也有扩散流量密度不由梯度支配的情形，这时可将其（Jd）并入S中。连续性方程：（2）直角张量形式和法则：一项某下标重复，则该下标依次取1、2、3然后做和运算。如何把特定的控制微分方程改写成通用形式。把相关因变量的非稳态项，对流项及扩散项转化成标准形式；把扩散项内梯度前的系数取为的表达式；把其余所有项之和置于方程右端，定义为源项。通用方程可以是有量纲形式，也可以是无量纲形式，、S也相应无量纲化。作用：通用方程通用数值方法公式通用计算机程度2-2 坐标的性质及控制方程的类型一、坐标的性质坐标自变量。坐标（系）或自变量的数目对问题的难易程度有很大影响，而对一定的问题而言，坐标（自变量）数目是可变的，既可以用这个坐标系来描述，又可以用别的坐标系来描述。1、自变量的作用：在数值计算中，将选择用来计算值的自变量值，所需计算值的位置的多少与自变量的数目有关，自变量数目所需计算的位置。因变量与自变量的相对性，适用于温度场是坐标的单调函数的情形。2、坐标的选择，原则：恰当，合适。自变量数目最少的坐标系网格节点数少。（1）选择最简单的坐标系圆柱（简）中的轴对称导热：圆管内的轴对称流动：（2）运动坐标系准稳态的概念，（3）利用充分发展的概念：存在这样一个坐标，当过程发展到一定深度后，因变量的无量纲分布与该坐标无关。，充分发展后平面自由射流；但（4）相似变换：减少自变量数目的变换统称相似变换。半无限大物体（x0）在1B.C.下的非稳态导热：3、坐标的单、双向性采用单向坐标和双向坐标的概念，可以形象地描绘不同类型控制方程的物理作用上的区别。单向坐标：在一个坐标轴上，如果扰动（或影响）只能向一个方向传递，则称此坐标为单向坐标。何谓向一个方向传递？坐标上任意给定位置处因变量之值只受该位置一侧条件变化的影响，且该点因变量之值也只对其一侧位置上的固变量值发生影响。（“抛物型”表示一种“单向作用”的概念）时间坐标是一个典型的单向坐标。高温固体的冷却，其在某瞬时的温度只受该瞬时以前的条件的影响。双向坐标：在一个坐标轴上扰动（或影响）可以向两侧传递，称为双向坐标。何谓向两个方向传递？坐标上任意给定位置处因变量之值要受该位置两侧条件变化的影响，且该处因变量之值也会对其两侧位置上的因变量值发生影响。空间坐标是典型的双向坐标。但在一定条件下，空间坐标也可以成为单向坐标；在流体流动时，如果在某一个坐标上有很强的单向流动，则重要的影响只能从上游传播到下游，某处状态也主要受其上游条件的影响，受下游条件影响很小。总结：对流是一种单向过程，而扩散是一种双向过程：在流体流动时，二者同时存在，仅当流作用很强（流量很大）时，扩散作用可忽略不计，空间坐标才近似成为单向坐标。二、控制方程的类型讨论控制方程的类型对数值解的影响1、能量方程导热：对流换热：（已假定：）若常数，OK时h=0，则组成：非稳态项，对流项，扩散项，源项控制方程：微分容积中守恒原理的数学描述。2、能量方程的守恒性质数值计算中的控制容积无论多么小，总是有限容积，而不是微分容积。问题：对任意有限容积，控制方程是否也一定满足守恒原理。守恒型控制方程，对任意大小的有限容积能使守恒原理得到满足的控制方程。非 不能可以证明，在直角坐标中，当对流项为散度形式时，控制方程是守恒的。将上面的控制方程对任一有限容积V作积分有限容积V上的能量平衡原理的表达式：单位时间内有限容积中物质内能的增加量对流带入的净能量率导入的净能量率生成的能量率若用非守恒型控制方程，有： 简证：取平行六面体 3、控制方程的分类（1）从数学上分类，限于二阶偏微分方程，对二元二阶线性偏微分方程，下标xy表示对该变量的导数，a、b、c、d、e、f、是x、y的函数，对求解区域R，方程的特性由系数之值决定，在区域内任意一点有两条实的特征线 双曲型，在区域内任意一点有一条实的特征线 抛物型，在区域内任意一点无实的特征线 椭圆型三类方程在数学上的主要区别：影响区域和依赖区区域不相同。影响区域、依赖区域；R中任一点P的依赖区域，为了唯一定之道，必须给函数值的点B依赖区域的集合；P的影响区域：变化时，函数值发生变化的点的集合。椭圆型方程（即求解区域上每一点都是椭圆型的）：对于能量方程无非稳态项，自变量与时间无关。特点：无特征线，任一点P的依赖区域是包围该点的区域封闭边界曲线，而P点的影响区域则是整个求解区域。边值问题，对应于物理学上的一类平衡问题，或稳态问题。求解区域内各点处的因变量值是相互影响的，与双向坐标相联系。结果：离散代数方程必须联立求解（直接解法或迭代解法），而不可能把区域中某一部分的值求得后再去确定其区域上的值。抛物型方程：因变量与时间有关，或问题中存在类拟于时间的自变量，与单向坐标相联系，能量方程有非稳态项。对应于：物理学上一类非平衡问题，或非稳态问题，步进问题。特点：过区域中任一点P有一条实特征线，其方向与单向坐标相垂直，如图，P点的依赖区域和影响区域以特征线为分界线。对非稳态问题：某一瞬间物体中的温度（分布）取决于该瞬时以前的情况及边界条件，而与该瞬时以后将要发生的情况无关，反之，某一时刻的温度只影响此后的温度分布边界层类型的流动与换热：忽略主流方向的扩散作用，下游的物理量（u,v,t）取决于上游，上游 只会影响下游的物理量。结果：不必将单向坐标上所有位置处的离散方程联立求解！只需从某一初始值出发，结合边界条件，沿单向坐标一步步向前推进步进法步进问题。大量节省计算机的内存及计算时间，二维存储一维，三维存储二维。双曲型方程：依赖和影响区域与抛物型方程相同，即：依赖区域位于运动方向的上游，影响区域位于下游。不同之处：某点的依赖区域或影响区域，不是其上游或下游区域的全部，而只是过该点的两条特征线之间的区域，如图。例：无粘流体的非稳态流动，无粘流体的稳态超音速流动。大多数工程导热、对流换热问题都属于椭圆型或抛物型方程，本课程只限于讨论这两类方程。（2）从物理现象上划分在传热相应于椭圆型与抛物型方程的物理过程，有专门的称谓。抛物型有一个空间坐标是单向坐标边界层型（流动或换热）问题；椭圆型所有空间坐标都是双向的回流型（ ）问题。三、一维模型方程控制方程由四部分组成：非稳态项、对流项、扩散项、源项。在研究建立离散方程和方法时，为了避免复杂化，不必着眼于完全的方程，只需把同一类型的项取出一项作为代表来分析，得到一个维非稳态对流扩散方程，即在对流换热问题数值计算研究中广泛采用的所谓一维模型方程。非守恒型：守恒型：广义变量：u、v、w、h、c等广义扩散系数，对于传热S广义源项，代表一切不能归并到非稳、对流，扩散项中去的量，不定是物理上的真正源项。对于导热问题：C=常数，由对流换热C=常数，h=cT第三章 离散化方法3-1传热问题数值求解的基本步骤一、数值解法的基本思想数值解法是一种离散近似的计算方法。这种方法得到的是求解区域中某些代表性位置上未知（待求）物理量（速度、温度、浓度等）的近似值，而不是象分析解或近似分析解那样的连续函数，“数值解法”一词即由此而得名。例：大平板（-LxL）在第三类B.C.下的一维非稳态导热B.C.：i.C.：用分离变量法求解可得：令式中，分别是下列（特征）方程的根：，毕渥数，物理意义，内部导热热阻/外部换热热阻对不同的Bi值，有不同的值系列。结论：采用计算机进行数值计算不仅是求解偏微分方程的有力工具，而且对一些经验公式和用无穷级数表示的分析解，也常常需要用计算机来获得数值结果。如何进行数值计算，基本思想？把原来在时间，空间坐标上连续的物理量场（速度场，温度场，浓度场等），用求解区域中有限个离散点上的值的集合来代替，并按一定方式建立起关于这些值的代数方程，求解代数方程以获得物理量场的离散近拟解。二、基本步骤例：长方柱体中的导热，已知：稳态，四个侧表面各自维持均匀温度，用数值方法求柱体中的温度分布。1建立物理模型，对问题作必要的简化。(1) 假定涉及的温度变化范围不大，常物性；(2) 柱体长度方向端部效应忽略不计，二维问题稳态：二维、常物性、无内热源的导热问题，1st B.C.。2. 建立相应的数学模型：控制方程+单值性条件有的工程实际问题的模型化工作较难，需要对实际问题进行仔细分析，还需要经验。例：离散电阻片热源特点：温度分布在一定的深度的不均匀，畸变温度分布的不均匀性呈现周期性。所以只研究一个周期区间即可！3. 区域离散化：在计算区域中配置需要计算温度的地点（位置）（称为节点）这一步骤称为区域离散化。4. 控制方程的离散化：按照一定的原则（通常是守恒原理），建立每个节点上未知量（，温度）与其邻点上未知量之间的代数关系式（称为离散方程）控制方程的离散化。例如，由能量方程可以得出每一节点温度与相邻节点温度之间的代数关系。5. 求解所得到的代数方程，获得节点上的未知量之值。6. 对所获得的数值结果进行分析、比较和讨论。基本步骤框图（流程图）：三、传热问题的主要数值方法1、有限差分法（Finte Difference Method, FDM）2、有限元法（Finite Element method, FEM）3、边界元法（Boundary Element Method, BEM）4、有限分析法（Finite Analysis Method, FAM）应用最多的是有限差分法，与有限元法相比，差分法占有一定优势：方法发展的成熟程度；实施的难易程度；应用的广泛性。3-2 区域离散化方法定义：用一系列与坐标轴平行的曲线将计算区域划分成很多互不重叠的子区域，并选定每个子区域中的节点的过程。每一个节点可视为相应微小容积（称为控制容积）的代表，控制容积的边界称为界面（常用虚线表示）。沿坐标轴方向联结相邻两节点而形成的曲线簇称为网格线（常用实线表示）。一、分类 根据节点在子区域中的位置进行划分1、外节点法节点位于子区域的顶点用直线簇划分出子区域。节点：子区域的顶点即直线簇的突点网线线：直线簇界面：网格间距（相邻两节点间的网格线段）的中垂线；控制容积：不是子区域特点：子区域与控制容积不重合，先确定节点，后确定界面，practiceA，方法A2、内节点法节点不在控制容积中心，节点位于子区域的中心界面：划分子区域的直线簇控制容积：=子区域节点：控制容积的中心网格线：特点：先定界面，后d定节点,practiceB方法B。二、网格系统的标记方法以二维为例1 i-j-n系统2 P、E、W、N、S系统节点间距（网格间距），界面间距对于均分网格：，内、外节点法的节点分布在区域内部趋于一致。例：一维网格系统三、两种区域离散化方法的比较1、边界节点所代表的控制容积不相同A：，能更好地考虑边界节点之间的传热作用，更符合实际；B：，不能考虑边界节点之间的传热作用，（适合于今后处理2、3类B.C.的附加源项法）2、网格不均匀时界面位置A：界面位于相邻两节点正中间：B、否对导热量计算的准确度不相同 节点位置A：不在控制容积中心B：在控制容积中心、对求解区域内材料物性突变的适应能力不同A:B: 易方便地将物性发生阶跃变化的交界面作为控制容积的界面，从而使同一控制容积内的物性保持均匀一致。4、对B、C突变的适应能力不同A：不便于适应B：便于适应结论：两种方法都得到应用，但当有材料物性突变时，建议采用内节点法。3-3 建立离散方程的方法建立离散方程的常用方法有四种，即泰勒（Taylor）级数展开法、多项式拟合法、控制容积积分法及平衡法，这里主要介绍泰勒级数展开法和控制容积积分法。一、泰勒（Taylor）级数展示法定义：把控制方程中的各阶导数用相应的差分表达式来代替而形成离散方程的方法。由于各阶导数的差分表达式可以由泰勒级数展开而导出，故名曰泰勒级数展开法。以一维模型方程为例（1）为简便计，假定为常数（2）1、导数的差分表达式建立均匀时间空间网格：在x-坐标系的求解区域中，用划分网格，把在时一空网格中节点（i，n）处之值记为t(i,n), 相应地有、将作泰勒级数展开：.（3）（4）由式（3）移项可得同理，由式（4）移项可得由式（3）式（4）可得说明：O（rder）表示数量级之意，O（x）,O(x)2表示未写出的高阶项之和的数量级，常称为截断误差。O(x)表示截断误是x的数量级，O(x)2则表示截断误差是（x2）的数量级。特别注意，O（x）并未给出截断误差的准确值，而只表明截断误差是如何随x0而变小的，例如O（x）即表示截断误差随x0而按x减小，O(x)2则表示截断误差随x0而按（x2）减小。所以O(x)2比O（x）高一个x的阶，当x足够小时，截断误差为O(x)2的表达式比截断误差为O（x）的表达式更准确。这里，t(i,n)、t(i+1,n)和t（i-1,n）是函数在节点（i,n）、(i+1,n)和（i-1,n）处的精确值，它是未知、待求的。在进行数值计算时，只能获得其近似值，分别记为且它们应满足以下关系式(或按以下关系式来确定近似解：向前差分一阶精度向后差分一阶精度中心差分二阶精度向前、向后差分双涉及所论节点一侧的节点函数值，称为单侧差分或偏差分中心差分涉及所论节点两侧的节点函数值几何表示：一般规律：用到的节点函数值越多，涉及的节点区域范围越大，提供的特定变量的变化信息越充分、完整，差分表达式的准确度越高！请自证：同理，对非稳态项的一阶偏导数也有类似的差分表达式：对二阶导数项亦有：；自证：差分表达式正确性的检验；各阶导数的差分表达式应当满足两个基本要求（1）函数为常数时，差分表达式亦应成立（差分表达式中分子各项系数的代数和为零）；（2）差分表达式的量纲必须与导数的量纲一致。2、控制方程的差分表达式离散方程的建立有各阶导数的差分表达式后，如何得到控制方程的差分表达式；把控制方程中的所有导数用同一节点处的相应差分表达式来代替，并将所有各项的截断误差相加，从而得到要求的差分离散方程和相应的截断误差。问题：在一个时间步上，按步进方式采用向前差分很容易获得相应的差分表达式如何处理？在间隔内，也是随时间变化的。怎样选择展开点？有三种常用的选择三种常用格式（1）显格式：，选在时层的开始时刻，用中心差分；（2），全隐格式，展开点选在时层终了时刻仍用中心差分；（3）C-N格式，展开点选在时层的正中间，仍用中心差分；关键是，将在（）处展开。两式相加，所以注意：的差分表达式，对节点（）展开是中心差分，对再用中心差分，所以思考：以上三种典型格式的物理意义。二、多项式拟合法导数的差分表达式也可以通过多项式的拟合来获得，这相当于对未知函数的局部变化曲线（型线）采用多项式来逼近。1、线性拟合：假定函数在节点（i、n）附近对x的变化是近似的线性关系（线性化），则有：取（i，n）的x坐标，则由此二式得的向前差分为2、二次拟合，且仍取于是有由此解得：可得一、二阶导数的差分表达式分别为：在流动与传热问题的数值计算中，多项式拟合法主要用来处理边界条件，例如当由已知的节点温度确定边界上的热流或由给定的热流确定节点温度。例：在图示情形中，已知边界（节点）和区域内部节点的温度，材料导热系数，试用多项式拟合法确定边界上的热流。解：（1）线性拟合边界附近温度成线性变化，则（2）边界附近温度分布的二次曲线拟合，解得于是得与有二阶精度的单侧差分格式一致对于采用离散化方法B的网格：此二式相减，反之，如果给定了边界上的热流密度qB，则可用多项式拟合法由qB及内节点温度来计算边界温度。例如，从上面两个二次曲线拟合所得的qB表达式，可以立即写出qB和所表达的：方法A方法B对方法A的网格，采用三次多项式拟合时，此式的截断误差可通过将对节点（i,1）作Taylor展开来获得。三、控制容积积分法控制容积积分法又称有限容积法，是传热问题数值计算中广为采用的方法。控制容积积分法的主要实施步骤可以归纳如下：1、步骤（1）将节点视为控制容积的代表，把守恒型控制方程在任一控制容积和时间间隔内对空间变量和时间变量作积分。（2）选定函数（温度）随空间坐标和时间坐标的变化规律，即分布曲线（常称形状函数、型线）。（3）根据选定的型线，作出方程中的各项积分，并整理成关于未知节点函数值（温度值）的代数方程。在举例说明控制容积积分法的具体实施过程以前，先介绍两各常用的型线，即阶梯式分布和分段线性分析。2、两种常用型线（a）t-x (b)t-在t-的阶梯式分布中，又有显式与隐式两种情形。显式：在的整个时间间隔内，函数（温度）均取为开始时刻之值，仅在该间隔的结束时刻才取得终了时刻（+）之值；隐式：在整个时间间隔内，函数（温度）均取为终了时刻之值。CN：在间隔内，从开始时刻值直线变化到终了时刻之值。3、实施过程控制容积P，y、z方向均为单位长度，（1）将模型方程对控制容积P在时间间隔内积分体现了间隔内控容P上的能量平衡关系。（2）选择的型线，完成积分计算非稳态项，t及对x作阶梯式变化，扩散项，对作显式、阶梯式变化，间隔内按时刻计算，得式中，必须转换为用节点温度值来表示，涉及t随x的变化方式，取为分段线性分布，则有式中，分别是界面e,上的导热系数。源项 S对及x均作阶梯式变化，且对为显式变化。（3）将各项积分结果代入能量平衡式，整理得到离散方程（两端同除以）式中主节点系数；邻点系数，b-常数采用均匀网格时，常物性改用（i、n）网格标记，则上式将与泰勒级数展开法得到的显格式结果相同。关于型线选择的说明：型线是一种辅助关系式，一旦离散方程建立起来，型线就失去了意义，因此不必追求一致性，即对不同的物理量可以选择不同的型线；对同一物理量在不同的项中，可以选择不同的型线。选择的原则：便于积分运算；使离散方程有较好的数值特性（显，隐？）四、控制容积平衡法把守恒定律直接应用于所研究的控制容积，并把节点视为控制容积的代表，即可导出节点上未知物理量值之间的代数关系式。对于模型方程所表示的一维扩散问题，守恒定律（热力学第一定律）表现为：时间间隔内控制容积P中能量（热能）的增量，等于同一时间间隔内扩散进入该控制容积的净热量及内热生成量之和，故有进一步假定：1、上式右端各项取时刻之值（显式）；2、界面导数亦按分段线性计算；3、对x均作阶梯式变化与控容积分法的型线相同则经整理得到的节点温度代数方程与控制容积积分法相同。五、建立离散方程方法的比较1、泰勒级数展开法和多项式拟合法偏重于从数学角度进行推导，把控制方程中的各阶导数用相应的差分表达式代替，易于对离散方程的数学特性进行分析（如分析截断误差等）。缺点：非均匀空间网格时的离散方程形式比较复杂；导出过程的物理概念不够清晰； 不能完全保证所得差分方程具有守恒特性。2、控制容积积分法和平衡法：着重于从物理观点分析，推导过程的物理概念清晰，离散方程系数有一定的物理意义，所得到的离散方程都是有限大小容积上某种物理量与守恒的表达式，具有守恒特性。缺点：不便对方程进行数学特性分析。3-4 差分方程（离散方程）的特性分析无论采用四种建立离散方程方法中的哪一种方法，在推导过程都作了近似处理，必然会引入误差，直接影响数值解的准确度，从数学的角度看，这些误差包括：差分方程的截断误差相容性问题差分方程解的截断误差收敛性问题计算过程中引入的舍入误差稳定性问题另一方面，希望所获得的数值解能保持热物理过程和现象原有的一些基本属性，例如在有限容积范围内的守恒特性，在纯对流问题中扰动仅沿流动方向传递的特性。一、数学特性1、截断误差和相容性讨论误差，必须引入差分方程精确解的概念。差分方程精确解离散方程求解过程中不引入任何误差的解，其所采用的计算机是无限位字长的。考虑最简单的一维扩散问题：（1）符号表示对函数在点（i，n）作规定微分运算的算子，例如，对上述一维扩散问题，有微分算子：而，就是节点（i、n）处的方程（1）符号表示对节点函数作某些差分运算的算子，例如（2）于是，就代表了方程（1）的差分格式。一个差分方程的截差是指其差分算子与相应微分算子之差，记为，差分方程的截差可以通过对差分方程的精解解作Taylor展开来导出（特别对控容积分法和平衡法导出的差分方程，并假设差分方程精确解满足作Taylor展开的条件）。例如，对方程（1）的差分格式，（2）把在点（i，n）作Taylor展开，并代入该差分方程，整理得到：所以 （当时）定义：当时间和空间网格的步长（间距）均趋于零时，如果差分方程的截断误差0，则称此差分方程与微分方程相容。相容性意味着当时空网格间距均0时，差分方程逼近微分方程，差分方程所数值模拟的确系原控制微分方程。相容性条件：截差呈的形式，且m、n均大于零时，差分方程具有相容性。注意：当截差表达式中含有项时，相容性仅在一定条件下才能满足。例：一维非稳态导热的Dufort-Frankel格式控制方程Richardson格式：三层格式思路：中心差分二阶精度对时间t的一阶偏导数用中心差分。，意指时刻满足， =其中差分方程理查森格式是一个显格式，但与前面提到的前三种格式不同：前三种格式：由n层（n+1）时层，涉及两个时层，称为二层格式。Richardson格式：由n-1、nn+1，涉及三个时层，称为三层格式。时，相容，该差分方程数值模拟的确是，但该格式是完全不稳定格式，无实用价值。如何使不稳定稳定？Dufort,Frankel格式将用代替，则差分方程变为无条件稳定的三层显式格式。出现的问题：格式的精度；要求的速度比的速度快，为什么所以讨论：当（常数）当时，且0，相容，数值模拟的是导热方程，但截差，则格式逼近的是如下双曲型方程不满足相容性要求的差分方程模拟的是别的微分方程，数值解没有意义。2、离散误差和收敛性离散误差：网格点上（i，n）差分方程的精确解与该点处微分方程的精确解的偏差称为差分方程解在该点处的离散误差，记为；当时间和空间网格间距和均0时，如果所有网格节点（i、n）上的离散误差均0，则称该差分方程的解是收敛的，或者说解具有收敛性。相容性是收敛性的必要条件，但不是充分条件。舍入误差：计算机的字长总是有限的，在运算过程中必然有四舍五入，产生舍入误差，节点（i，n）上实际求得的数值解与其对应的差分方程的精确解之间的偏差，记为；数值解的总偏差；计算实践表明，是的主要成分。的影响因素，m、n均0（相容性要求！）一般而言，m、n之值越大（即截差阶数越高），；在一定截差下，网格加密（），。但决不是截差越高，网格越密，数值解的准确性就越高！网格加密，运算次数，机时，对解的准确度造成影响；解的准确性受所有节点上截差的影响，而在近边界的内节点上通常不易获得高阶截差的离散方程，故无论网格如何加密，解的准确性也受此约束！所以通常采用具有二阶截差的格式究竟收敛性的条件或者判据是什么？这里不拟直接就收敛性进行证明，而引述一个定理，Lax等价性定理：对于线性初边值问题，如果问题是适定的，并且差分方程（格式）满足相容性条件，那么差分格式的稳定性就是该格式收敛性的充分而必要的条件。所以集中注意力格式的稳定性讨论！3、误差传播（累积）和稳定性用步进法求解非稳态问题（广而言之，抛物型问题）时，数值解是在单向坐标上按步长（对时间单向坐标是按时层），逐渐推进的，以时间坐标为例，如果在某一时层（n）上引入误差，则会影响到其后的各时层的计算结果，计算过程误差的引入是不可避免的，例如：舍入误差，初始条件、边界条件给定的误差等，定义：对于一定的差分格式，如果在某一时间层上引入的误差在其后逐层的计算中能得到控制，即逐步消失或保持有界，则此差分格式是稳定的，否则是不稳定的。不稳定的差分格式没有实用价值，因为在这种情形下，累积误差将随时层的而逐步，以至差分方程的真解终将被误差所淹没。对常系数差分方程而言，任意时层上引入误差产生的影响最终可以归结为在初始时层上引入的误差对数值解所造成的影响，即初值不稳定性。稳定性与否是格式本身所固有的？举例说明显格式，均匀网格，引入，则格式变为（网格傅里叶数）为简化计算，取，则问题的差分格式为假定：（1）B.C.精确，无误差（2）在初始层（n=0，初始条件），某节点（i,0）产生误差，而在初始层其它节点上无误差；（3）在以后各层的计算中，完全准确无误。这是一种最轻微的情况，看这唯一的误差如何影响节点温度的发展与传播？差方程的精解解设为，则偏差，它将满足下列误差方程：看由引起的误差传播50000040000300020010n=0ni误差逐层减小，并将最终消失，所以当F=1/2时，显格式是稳定的。改取，则有误差传播？图（自己完成），结论：误差逐层增大，所以F=1时，显格式是不稳定的（I，B，C）。4、差分格式的稳定性理论分析（1）矩阵分析法显格式：A全隐格式：C-N格式：上述三种典型格式的矩阵形式可统一表示为：式中，H是三对角矩阵，其元素依赖于FA，称为增长矩阵，b,是相应的常数向量。（N-1）（N-1）假定B.C.准确给定，只是在初始层上引入了误差，而在其后各时层的计算中未再引入其它误差，则由此而得到的解应满足下列方程：近似解可以导出误差：应当满足的方程为描述误差的传播规律，称为误差方程，显然，它是差分格式的齐次形式。显格式：右端的不等式显然对一切，都成立，左端：所以，时，显格式稳定！全隐格式，H=B-1，B=2I-A设A的特征向量为，则： =所以B的特征值 所以H=B-1的特征值，为：恒1，无条件稳定C-N格式：H=（I+A）（I+B）1恒1，无条件稳定。加权格式：，无条件稳定，条件稳定，稳定条件为（2）傅里叶级数法von.Neumann方法一维问题：0，L，利用误差的网格函数（定义在网格节点上），在0，L上定义一个函数，B、C准确给定，定义为阶梯函数。取为将周期性地延拓到整个x轴，并展开成复数形式的傅里叶级数，得系数是时层数n的函数，具有谐波波幅的意义，用来表示：因为所讨论的问题是线性的，满足迭加原理，故考察误差随时层的变化时，只需取傅里叶级中