计算流体力学一维稳态导热编程作业

课后习题第一章1.计算流体动力学的基本任务是什么计算流体动力学是通过计算机数值计算和图像显示，对包含有流体流动和热传导等相关物理现象的系统所做的分析。

2.什么叫控制方程？常用的控制方程有哪几个？各用在什么场合？流体流动要受物理守恒定律的支配，基本的守恒定律包括：质量守恒定律、动量守恒定律、能量守恒定律。

如果流动包含有不同组分的混合或相互作用，系统还要遵守组分守恒定律。

如果流动处于湍流状态，系统还要遵守附加的湍流输运方程。

控制方程是这些守恒定律的数学描述。

常用的控制方程有质量守恒方程、动量守恒方程、能量守恒方程、组分质量守恒方程。

质量守恒方程和动量守恒方程任何流动问题都必须满足，能量守恒定律是包含有热交换的流动系统必须满足的基本定律组分质量守恒方程，在一个特定的系统中，可能存在质的交换，或者存在多种化学组分，每种组分都需要遵守组分质量守恒定律。

4.研究控制方程通用形式的意义何在？请分析控制方程通用形式中各项的意义。

建立控制方程通用形式是为了便于对各控制方程进行分析，并用同一程序对各控制方程进行求解。

各项依次为瞬态项、对流项、扩散项、源项。

6.CFD商用软件与用户自行设计的CFD程序相比，各有何优势？常用的商用CFD软件有哪些？特点如何？由于CFD的复杂性及计算机软硬件条件的多样性，用户各自的应用程序往往缺乏通用性。

CFD商用软件的特点是功能比较全面、适用性强。

具有比较易用的前后处理系统和其他CAD及CFD软件的接口能力，便于用户快速完成造型、网格划分等工作。

具有比较完备的容错机制和操作界面，稳定性高。

可在多种计算机、多种操作系统，包括并行环境下运行。

常用的商用CFD软件有PHOENICS、CFX、SRARCD、FIDAP、FLUENT。

PHOENICS除了通用CFD软件应该拥有的功能外，PHOENICS软件有自己独特的功能：开放性、CAD接口、运动物体功能、多种模型选择、双重算法选择、多模块选择。

1. 已知有限体积法求解的通用控制方程为()()()u div div grad S tρφρφφ∂+=Γ+∂ 其中φ为通用变量，可代表u 、v 、w 、T 等求解变量。

（1）试说明通用控制方程中各项的物理意义；（2）对于特定的方程，φ、Γ、S 具有特定的形式，对应于质量守恒方程、动量守恒方程、能量守恒方程、（多种化学组分的）组分质量守恒方程，试写出φ、Γ、S 的具体表达式。

解答：（1） 方程中的各项（从左到右）分别是瞬态项、对流项、扩散项及源项。

（2）2. 简述有限体积法建立离散方程时应遵守的四条基本原则。

解答：1) 控制体积界面上的连续性原则当一个面为相邻的两个控制体积所共有时，在这两个控制体积的离散方程中，通过该界面的通量（包括热通量、质量通量、动量通量）的表达式必须相同。

即：通过某特定界面从一个控制体积所流出的热通量，必须等于通过该界面进入相邻控制体积的热通量，否则，总体平衡就得不到满足。

2) 正系数原则中心节点系数aP 和相邻节点系数anb 必须恒为正值。

该原则是求得合理解的重要保证。

当违背这一原则，结果往往是得到物理上不真实的解。

例如，如果相邻节点的系数为负值，就可能出现边界温度的增加引起的相邻节点温度降低。

3) 源项的负斜率线性化原则源项斜率为负可以保证正系数原则。

从式(C2)中看到，当相邻节点的系数皆为正值，但有源项Sp 的存在，中心节点系数aP 仍有可能为负。

当我们规定Sp ≤0，便可以保证aP 为正值。

4) 系数aP 等于相邻节点系数之和原则当源项为0时，我们发现中心节点系数等于相邻节点系数之和，而当有源项存在时也应该保证这一原则，如果不能满足这个条件，可以取SP 为0。

3．什么是对流质量流量F 、扩散传导量D 以及Pelclet 数Pe ，试用定义式表达之。

解答：F 表示通过界面上单位面积的对流质量通量(convective mass flux)，简称对流质量流量；D 表示界面的扩散传导性(量)(diffusion conductance)。

第二章部分答案-稳态导热2-46. 一厚度为7cm 的大平壁，一侧绝热，另一侧暴露于温度为30℃的流体中，其内热源热量为5103⨯W/m 3。

已知该平壁材料的导热系数为18K)W/(m ∙，平壁与流体间的对流表面传热系数为450)K W/(m 2⋅，试确定该平壁中的最高温度位置及其温度值？解：(1) 该题为具有内热源的一维平壁稳态导热问题，导热微分方程式为：022=Φ+λ dxt d 边界条件为：0=x ，0=dxdt ； δ=x ，∞+Φ==t t t w λ2 （根据热平衡求得：δΦ=-∞ )(t t h w ) 解方程，并代入边界条件得温度场为：∞+Φ+-Φ=t h x t δδλ )(222(2) 该平壁中最高温度在0=x 处(即0=dxdt )：117.5 30450)107()103()107(182103225252=+⨯⨯⨯+⨯⨯⨯⨯=+Φ+Φ=--∞t h t δδλ ℃ 2-47 核反应堆的辐射防护壁因受γ射线的照射而发热，这相当于防护壁内有ax e -Φ=Φ0 的内热源，其中0Φ 是X=0的表面上的发热率，a 为已知常数。

已知x=0处t=t1,x=δ处t=2t ,试导出该防护壁中温度分布的表达式及最高温度的所在位置。

导热系数λ为常数。

解：由题意导热微分方程0022=Φ+-ax e dx t d λ又x=0处t=t1,x=δ处t=2t积分并结合边界条件可得λδλλλδ2012020210a t x a e a t t a e t a ax Φ++Φ-Φ+--Φ=-- 令0=dx dt 可得：当()⎥⎦⎤⎢⎣⎡-+Φ--=-δδλδa e t t a a x a 1ln 1021时，t 最大。

2-48 核反应堆中一个压力容器的器壁可以按厚为δ的大平壁处理。

内表面（x=0处）绝热，外表面维持在恒定温度2t 。

γ射线对该容器的加热条件作用可以用一个当量热源Φ 来表示，且ax e -Φ=Φ0 ，a 为常数，x 是从加热表面起算的距离。

1. 已知有限体积法求解的通用控制方程为()()()u div div grad S tρφρφφ∂+=Γ+∂ 其中φ为通用变量，可代表u 、v 、w 、T 等求解变量。

（1）试说明通用控制方程中各项的物理意义；（2）对于特定的方程，φ、Γ、S 具有特定的形式，对应于质量守恒方程、动量守恒方程、能量守恒方程、（多种化学组分的）组分质量守恒方程，试写出φ、Γ、S 的具体表达式。

解答：（1） 方程中的各项（从左到右）分别是瞬态项、对流项、扩散项及源项。

（2）2. 简述有限体积法建立离散方程时应遵守的四条基本原则。

解答：1) 控制体积界面上的连续性原则当一个面为相邻的两个控制体积所共有时，在这两个控制体积的离散方程中，通过该界面的通量（包括热通量、质量通量、动量通量）的表达式必须相同。

即：通过某特定界面从一个控制体积所流出的热通量，必须等于通过该界面进入相邻控制体积的热通量，否则，总体平衡就得不到满足。

2) 正系数原则中心节点系数aP 和相邻节点系数anb 必须恒为正值。

该原则是求得合理解的重要保证。

当违背这一原则，结果往往是得到物理上不真实的解。

例如，如果相邻节点的系数为负值，就可能出现边界温度的增加引起的相邻节点温度降低。

3) 源项的负斜率线性化原则源项斜率为负可以保证正系数原则。

从式(C2)中看到，当相邻节点的系数皆为正值，但有源项Sp 的存在，中心节点系数aP 仍有可能为负。

当我们规定Sp ≤0，便可以保证aP 为正值。

4) 系数aP 等于相邻节点系数之和原则当源项为0时，我们发现中心节点系数等于相邻节点系数之和，而当有源项存在时也应该保证这一原则，如果不能满足这个条件，可以取SP 为0。

3．什么是对流质量流量F 、扩散传导量D 以及Pelclet 数Pe ，试用定义式表达之。

解答：F 表示通过界面上单位面积的对流质量通量(convective mass flux)，简称对流质量流量；D 表示界面的扩散传导性(量)(diffusion conductance)。

1 提出问题[问题描述]Sod 激波管问题是典型的一类Riemann 问题。

如图所示，一管道左侧为高温高压气体，右侧为低温低压气体，中间用薄膜隔开。

t=0 时刻，突然撤去薄膜，试分析其他的运动。

Sod 模型问题：在一维激波管的左侧初始分布为：0 ,1 ,1111===u p ρ，右侧分布为：0 ,1.0 ,125.0222===u p ρ，两种状态之间有一隔膜位于5.0=x 处。

隔膜突然去掉，试给出在14.0=t 时刻Euler 方程的准确解，并给出在区间10≤≤x 这一时刻u p , ,ρ的分布图。

2 一维Euler 方程组分析可知，一维激波管流体流动符合一维Euler 方程，具体方程如下： 矢量方程：0U ft x∂∂+=∂∂ (0.1)分量方程：连续性方程、动量方程和能量方程分别是：222,,p u ρ()()()()2000u tx u u pt x x u E p E tx ρρρρ∂⎧∂+=⎪∂∂⎪⎪∂∂∂⎪++=⎨∂∂∂⎪⎪∂+⎡⎤∂⎣⎦+=⎪∂∂⎪⎩ (0.2)其中 22v u E c T ρ⎛⎫=+ ⎪⎝⎭对于完全气体，在量纲为一的形式下，状态方程为：()2p T Ma ργ∞=(0.3)在量纲为一的定义下，定容热容v c 为：()211v c Ma γγ∞=- (0.4)联立(1.2)，(1.3)，(1.4)消去温度T 和定容比热v c ，得到气体压力公式为：()2112p E u γρ⎛⎫=-- ⎪⎝⎭(0.5)上式中γ为气体常数，对于理想气体4.1=γ。

3 Euler 方程组的离散3.1 Jacibian 矩阵特征值的分裂Jacibian 矩阵A 的三个特征值分别是123;;u u c u c λλλ==+=-，依据如下算法将其分裂成正负特征值：()12222k k k λλελ±±+=(0.6)3.2 流通矢量的分裂这里对流通矢量的分裂选用Steger-Warming 分裂法，分裂后的流通矢量为()()()()()()()12312322232121212122f u u c u c u u c u c w γλλλργλλλγλλγλ⎛⎫⎪-++ ⎪=-+-++ ⎪ ⎪ ⎪-+-+++ ⎪⎝⎭＋＋＋＋＋＋＋＋＋＋＋(0.7)()()()()()()()12312322232121212122f u u c u c u u c u c w γλλλργλλλγλλγλ⎛⎫⎪-++ ⎪=-+-++ ⎪ ⎪ ⎪-+-+++ ⎪⎝⎭－－－－－－－－－－－(0.8)其中：()()()223321c w γλλγ±±±-+=- c 为量纲为一的声速：22Tc Ma ∞=(0.9)联立(1.3)，(1.9)式，消去来流马赫数得：ργp c =3.3 一阶迎风显示格式离散Euler 方程组 10n n i i x i x i U U f f t xδδ+-++--++=∆∆ (0.10)得到()()n+1nj j 11U =U j j j j t f f f f x++---+∆⎡⎤--+-⎣⎦∆ 算法如下：① 已知初始时刻t=0的速度、压力及密度分布000,,j j j u P ρ，则可得到特征值分裂值0k λ±，从而求出流通矢量0j f ±；② 应用一阶迎风显示格式可以计算出1t t =∆时刻的组合变量1j U ，从而得到1t t =∆时刻的速度、压力及密度分布111,,j j j u P ρ；③ 利用1t t =∆时刻的速度、压力及密度分布111,,j j j u P ρ可得特征值分裂值1k λ±，从而求出流通矢量1j f ±；④ 按照步骤2的方法即可得到2t t =∆时刻的速度、压力及密度分布222,,j j j u P ρ；⑤ 循环以上过程即可得到()1t n t =+∆时刻的速度、压力及密度分布n+1n+1n+1,,j j j u P ρ。

一维稳态热传导方程
一维稳态热传导方程
一维稳态热传导方程是热传导问题的基本方程，其中，热传导系数K是用温度未变量的函数。

热传导方程的求解是热储存问题的基本算法，常用Finite element法和Finite Difference法求解。

一维稳态热传导方程
一维稳态热传导方程用来描述传统的单调方程热导率, 它表示
时间和空间的变化：
/d2T/dx2 + q_e = 0
其中：/d2T/dx2 为温度对x的二阶导数，q_e 为热源或热损失。

一般应用条件
1. 物质性质不变：
流体恒定的密度和热容，固体恒定的力学强度，导热系数和热容。

2. 境界条件：边界条件需要定义，自由表面是定义化学反应率的温度，或热损失。

热传导方程的求解
常用Finite element法和Finite Difference法求解。

Finite element法：采用有限元素方法，把区域分成小的等边形或三角形元素，然后求解每一个元素对应的温度场。

Finite Difference法：将方程化为一组一次方程组，然后由矩阵的求解得到方程组的解。

结论
一维稳态热传导方程是热传导问题的基本方程，其中，热传导系数K是用温度未变量的函数，然后常用Finite element法和Finite Difference法求解。

一维稳态导热方程推导导热方程简介在理解一维稳态导热方程之前，我们先来了解导热方程的概念。

导热方程描述了热量在介质中的传播规律，通过该方程可以推导出温度在空间中的分布情况。

一维稳态导热方程的定义一维稳态导热方程描述了在一个维度上，介质中的热量分布不随时间变化的情况下的热传导行为。

该方程可以用数学形式表示为：$$\\frac{{d}}{{dx}}(k(x)\\frac{{dT}}{{dx}}) = 0$$其中，k(x)是介质的导热系数，T(x)是温度分布函数，x表示一维空间坐标。

推导过程为了推导一维稳态导热方程，我们需要考虑热传导过程中的热流量和热量守恒。

假设我们有一个长度为L的介质，其中的热传导是沿着x方向进行的。

首先，我们考虑介质内的一个微小段dx，该段温度为T(x)，导热系数为k(x)。

根据热量守恒定律，该微小段内的热量变化率应与进入和离开该段的热流量之和相等。

我们可以将该微小段的热量变化率表示为：$$\\frac{{d}}{{dt}}(Q(x)) = -\\frac{{d}}{{dx}}(F(x))$$其中，Q(x)表示该微小段内的热能，F(x)表示通过该微小段的热流量。

考虑微小段dx与其邻近的微小段之间的热传导，我们可以使用傅里叶定律来表示热流量F(x)：$$F(x) = -k(x) \\frac{{dT}}{{dx}} dx$$代入热量守恒定律的表达式中，我们得到：$$\\frac{{d}}{{dx}}(k(x)\\frac{{dT}}{{dx}}) = 0$$这就是一维稳态导热方程。

方程的物理意义一维稳态导热方程描述了热量在一个维度上的传播规律。

方程表明，在稳态下，介质中的热传导是由温度梯度驱动的。

温度梯度越大，热传导越强。

导热系数k(x)描述了介质内导热的特性，它可以反映介质的性质和结构。

通过求解一维稳态导热方程，我们可以确定介质内不同位置的温度分布。

这对于许多工程问题和科学研究都是非常重要的，例如热传导问题、传热设备设计等。

作业(11月11日课堂上递交)。

题1: 一维非稳态无源导热问题 考虑一维非稳态无源导热问题：)(dxdT k dx d dt dT C P =ρ。

其中导热系数为s)K J/(m 10⋅⋅=k ，P C ρ为常数K)J/(m 1037⋅=P C ρ。

初始条件：t =0 s, 0≤x ≤0.02 m, T (x )=200 ºC边界条件：t >0 s, x =0 m, ∂T/∂x =0; x =0.02 m, T =100 ºC请用显式格式、C-N 格式和全隐式格式分别求解该问题(内节点法，均匀网格，∆x =0.004 m ，即五个网格)。

只需要提交以下四个图(请将四个图打印在一张纸上)：图1：画出t=40s 时的温度空间分布(线为精确解，符号1为数值解∆t =2 s ，符号2为数值解∆t =4 s ，符号3为数值解∆t =8 s)。

数值计算格式统一为显式格式。

图2：画出t=40s, 80s, 120s 三个时刻的温度空间分布(线为精确解，符号为数值解)。

数值计算的时间步长统一为∆t =2 s 。

数值计算格式统一为C-N 格式。

图3：画出t=40s, 80s, 120s 三个时刻的温度空间分布(线为精确解，符号为数值解)。

数值计算的时间步长统一为∆t =2 s 。

数值计算格式统一为全隐式格式。

图4：画出t=40s 时的温度空间分布(线为精确解，符号1为显示格式数值解 ，符号2为C-N 格式数值解，符号3为全隐式格式数值解)。

数值计算的时间步长统一为∆t =8 s 。

提示：精确解为∑∞=+---+=121)cos()ex p(12)1(400100),(n n n n x t n t x T γαγπ，)/(,02.0,2)12(P n C k m L L n ραπγ==-=题2: 准一维稳态有源导热问题提示：(见提示)(只需要写出内节点P 以及左右边界节点分别对应的三个离散方程 )(只需提交一张图，不用给出计算方法和计算程序) w e。

计算传热学大作业一维稳态矩形直肋问题一维非稳态无限大平壁导热问题一维稳态矩形直肋问题问题描述：等截面直肋稳态导热问题，图中t0 =100℃，t f =20℃，表面传热系数h= 50W /( m2·K )，导热系数λ=50W /( m·K )，肋高l=45mm，肋厚δ=10mm 。

1.加密网格，肋端绝热边界条件下计算程序编写矩形直肋的一维稳态、无内热源、常物性导热问题计算程序；计算等截面直肋的肋片效率。

2.肋端第三类边界条件下计算程序编写矩形直肋的一维稳态、无内热源、常物性导热问题计算程序；计算等截面直肋的肋片效率。

一．肋端绝热边界条件下1. 数学模型该问题属于一维稳态导热问题，常物性，无内热源。

其导热微分方程为0单值性条件为x=0,t0=100℃,肋端绝热。

2. 计算区域离散化时间离散（一维稳态，不存在时间离散）空间离散，划分多少N=45个区域.有N+1=46个点.3. 离散方程组对于内部节点（2≤i≤N+1）对绝热边界节点（i=N+1）4. 方程求解对内部节点（2≤i≤N+1）对绝热边界节点（i=N+1）求解：雅可比迭代5.肋片精确解及肋片效率C程序如下：#include <stdio.h>#include <math.h>void main (){int N=45,K=100000,i,N1=N+1,IT=0,TP;floatEPS=0.00001,T0=100.0,TF=20.0,h=50.0,LAMD=50.0,DT=0.01,T1[3000],T2[3000],L=0.045,TI=100,D TX=L/N,T3[3000];//给参数赋初值double m=sqrt(2*h/LAMD/DT),YT;//精确解求解公式printf("N=%d K=%d EPS=%6.5f T0=%6.2f TF=%6.2f h=%6.2f LAMD=%6.2f L=%6.2fDT=%6.2f DTX=%6.2f\n",N,K,EPS,T0,TF,h,LAMD,L,DT,DTX);//打印参数，方便查看for(i=1;i<=N+1;i++){T1[i]=T0;//内节点迭代计算初值}do{for(i=2;i<=N;i++){T2[i]=T1[i];//保留旧值T1[i]=((T1[i-1]+T1[i+1])*LAMD*DT+2*h*TF*DTX*DTX)/(2*LAMD*DT+2*h*DTX*DTX);//计算出内部各节点的温度}T1[N+1]=(DT*LAMD*T1[N]+h*DTX*DTX*TF)/(LAMD*DT+h*DTX*DTX);//计算出绝热边界点的温度TP=0;for(i=2;i<=N;i++){if(fabs(T2[i]-T1[i])>EPS) TP=1;//误差校核}if(TP==0) break;IT++;//进入下一次迭代}//完成do循环while(IT<=100000);if(IT==100001) printf("NO CONVERGENCE\n");else{printf("NO.ITERATIONS=%d\n",IT);//输出迭代次数总数for(i=1;i<=N1;i++){printf(" %6.2f",T1[i]);}printf("\n");}//输出每个节点温度值数值解YT= tanh(m*L)/m/L;//求肋片效率printf(" %6.2f",YT);//输出肋片效率printf(" \n");T3[1]=T0;for(i=2;i<=N;i++){T3[i]=0;}for(i=2;i<=N1;i++){T3[i]=TF+(T0-TF)*(cosh(m*(L-(i-1)*DTX)))/cosh(m*L);//求内部各节点的理论解}for(i=2;i<=N1;i++){printf(" %6.2f",T3[i]);//输出每个节点的理论解}} //结束运行结果如下迭代次数为5264次，肋片效率η=0.946. 解的分析将上述结果以折线图表示由分析可知，数值解与理论精确解的误差随深入肋片的距离而增加最大误差为6.88%，存在误差的主要原因是因为该理论精确解的计算公式主要针对长而薄的肋片，而题目中给出肋片为短而粗的肋片。

#include<math.h>#include<stdio.h>#define NT 201 //可计算的最大节点数#define Nt 200 //可输出的温度的最大时刻点数//子程序说明//void tdma(intN,double A[],double B[],double C[],double f[],double x[]); double lmda(double Tem);//主程序void main(void){//变量定义intN,i,j,Iter,It,Iout,Nout;double A[NT],B[NT],C[NT],f[NT],Tw[NT];double Tw0[NT],T0[NT],Tout[Nt][NT],Time[Nt];double dlt,TL,Tf,Ti,h,dens,Cp;double Tao,dtao,Taomax,dx,ap0;double lW,lE,lP,le,lw;double Err,Er,Eps;FILE *fp1,*fp2;//已知参数dlt=0.25; //壁厚TL=350.0; //左侧温度Tf=25.0; //流体温度Ti=25.0; //初始温度h=1500.0; //表面换热系数dens=150.0; //材料密度Cp=750.0; //材料的比热N=51; //网格节点数dtao=100.0; //时间步长-秒Taomax=20.0; //最大计算时间Eps=1.e-5; //计算精度Nout=10; //温度输出控制//计算参数dx=dlt/(N-1); //空间步长ap0=dens*Cp*dx/dtao; //系数ap0for(i=0;i<N;i++){T0[i]=Ti; //初始温度Tw[i]=Ti;Tw0[i]=Ti; //假设温度场}//初始温度场Tao=0.0; //初始时刻It=0; //时间步长的计数器Iout=0; //输出的个数Time[Iout]=Tao;for(i=0;i<N;i++)Tout[0][i]=Tw[i]; //初始温度//后续时刻的温度场计算do{Tao=Tao+dtao; //计算时刻It=It+1;Iter=0; //导热系数随温度变化的非线性迭代计算//左边界节点loop:A[0]=0.0;B[0]=1.0;C[0]=0.0;f[0]=TL; //左侧温度已知//内部节点for(i=1;i<N-1;i++){//节点处导热系数lW=lmda(Tw0[i-1]);lP=lmda(Tw0[i]);lE=lmda(Tw0[i+1]);//调和平均法计算界面的导热系数lw=2.*lW*lP/(lW+lP);le=2.*lE*lP/(lE+lP);//各系数A[i]=lw/dx;C[i]=le/dx;B[i]=A[i]+C[i]+ap0;f[i]=ap0*T0[i];}//右边界节点lW=lmda(Tw0[N-2]); //近右边界节点的导热系数lP=lmda(Tw0[N-1]); //右边界节点的导热系数lw=2*lW*lP/(lW+lP); //调和平均A[N-1]=lw/dx;B[N-1]=ap0/2.+h+lw/dx;C[N-1]=0.0;f[N-1]=h*Tf+ap0*T0[N-1]/2.;//调用TDMA计算新的温度tdma(N,A,B,C,f,Tw);//计算迭代间的最大误差Er=0.0;for(i=0;i<N;i++){Err=fabs(Tw[i]-Tw0[i]); //各节点的误差if(Err>Er) Er=Err; //各节点的最大误差}Iter=Iter+1; //迭代次数if(Er>Eps) //不满足精度，重新迭代计算{for(i=0;i<N;i++)Tw0[i]=Tw[i]; //取最新的温度值计算导热系数goto loop;}//满足精度输出本时刻参考点温度printf("Tao=%6.3fIter=%3dT2=%7.3fT8=%7.3f\n",Tao/3600.,Iter,Tw[1],Tw[7]);//保存给定时刻的温度if(fmod(It,Nout)==0.0) //每隔一定时间点输出该时刻的温度{Iout=Iout+1;Time[Iout]=Tao/3600.;for(i=0;i<N;i++)Tout[Iout][i]=Tw[i];}//设置下一时刻计算的温度初始值for(i=0;i<N;i++)T0[i]=Tw[i];}while((Tao/3600.)<Taomax);//将计算结果输出到作图文件中fp1=fopen("result1.dat","w");fprintf(fp1," ");for(i=0;i<Iout;i++)fprintf(fp1,"%8.3f",Time[i]); //时刻fprintf(fp1,"\n");for(i=0;i<N;i++){fprintf(fp1,"%6.4f",i*dx); //位置for(j=0;j<Iout;j++)fprintf(fp1,"%8.3f",Tout[j][i]); //温度fprintf(fp1,"\n");}fclose(fp1);//将计算结果输出到文本中fp2=fopen("result2.dat","w");fprintf(fp2," ");for(i=1;i<N/5+1;i++)fprintf(fp2,"%6.3f",i*dx*5); //节点位置fprintf(fp2,"\n");for(j=0;j<Iout;j++){fprintf(fp2,"%7.2f",Time[j]); //时刻for(i=1;i<N/5+1;i++)fprintf(fp2,"%7.2f",Tout[j][i*5]); //温度fprintf(fp2,"\n");}fclose(fp2);}//追赶法void tdma(intN,double a[],double b[],double c[],double f[],double x[]) {int i;double P[NT],Q[NT];P[0]=c[0]/b[0];Q[0]=f[0]/b[0];for(i=1;i<N;i++){P[i]=c[i]/(b[i]-P[i-1]*a[i]+1.e-20);Q[i]=(f[i]+a[i]*Q[i-1])/(b[i]-P[i-1]*a[i]+1.e-20);}x[N-1]=Q[N-1];for(i=N-2;i>=0;i--)x[i]=P[i]*x[i+1]+Q[i];return;}//导热系数与温度的关系double lmda(double T){double k;k=0.038+0.00023*T;return k;}。

《计算流体力学》练习题 Copyright by 李新亮一维算例问题1：正弦波的线性传播问题[问题描述]： 一匀速传播的正弦波，试用数值解给出t=T 时刻的速度场。

[控制方程]：)2sin()0,(0x x u xu t u π==∂∂+∂∂]1,0[∈x ；周期性边界条件[目的] 通过算例考核差分格式的耗散及色散误差 [计算网格] 均匀网格，20个网格点 [要求]A. 计算出t=2, 5, 10时刻的速度场B. 为了避免时间离散误差的影响，时间推进采用较高精度的方法（推荐采用如下三阶TVD 型 Runge-Kutta 方法）： 设)(u Q tu =∂∂，则)()1(nnu tQ uu∆+=)(4/14/14/3)1()1()2(u tQ u uu n∆++=)(3/23/23/1)2()2(1utQ uuunn ∆++=+其中 n u 和1+n u 分别为n 时间步及n+1时间步的u 值。

为了尽量减小时间离散误差的影响，建议采用很小的时间步长（例如 001.0=∆t ，或005.0=∆t ）C. 采用多种差分方法计算，对这些方法的结果进行比较 （同时画出精确解）。

问题2：方波线性传播问题（Jiang GS & Shu CW, J. Comput. Phys. 126,202-228,1996）.图1. 方波、尖波的传播[控制方程])()0,(00x u x u xu t u ==∂∂+∂∂]1,1[-∈x ，周期性边界条件其中初值为：该函数的曲线如图1所示，包含了一个方波，一个尖波及另外两个相对光滑的波。

[目的] 考察差分格式对于（导数）间断问题的分辨率。

[要求] 采用200个网格点，给出t=8时刻u 的分布，并与精确解进行对比。

为了减小时间误差的影响，推荐采用三阶Runge-Kutta 方法（公式见上文）或更高阶方法进行时间推进。

计算的CFL 数选为0.1或更小。

采用多种差分方法计算，对这些方法的结果进行比较 （同时画出精确解）。