（制冷及低温工程专业论文）热声谐振管中声振荡的格子boltzmann模拟.pdf

华 中 科 技 大 学 硕 士 学 位 论 文 i 摘要 热声制冷机是一种新型环保的制冷机，所具有的诸多优点，使其应用装置得到 广泛的研究。然而由于热声制冷机工作于非线性的热声耦合机制、多尺度效应及多 孔介质结构环境下，使得热声理论和传统的数值模拟具有一定的挑战性，急需开发 新的模拟分析方法。作者所在的课题组采用格子气对热声问题进行了探索性研究， 并取得了一定的进展，本文在此基础上采用格子 boltzmann 方法对热声谐振管中的 介质振荡及传热进行研究。 首先，针对热声谐振管的具体物理模型选择合适的格子 boltzmann 模型（d2q9 双分布热格子 boltzmann 模型），并根据相应的物理边界条件选择合适的格子 boltzmann 边界处理格式，同时将非平衡态外推边界处理方法应用到等温边界中。 其次， 采用d2q9总能双分布热格子boltzmann模型对声波衰减进行了数值模拟， 通过与分析解的比较，验证了所选择的格子 boltzmann 模型的合理性。并对活塞驱 动的谐振管稳定后的介质振荡特征进行了研究，分析了驱动频率对介质振荡密度幅 值的影响，重点分析了谐振管中介质振荡的非线性现象，并用实验结果进行验证， 模拟结果与实验结果的一致性证明了格子bolzmann方可用于研究谐振管中介质振荡 问题。 最后，对内置板叠的谐振管进行了研究，通过板叠前后介质振荡相位差及幅值 变化来分析板叠阻抗对流动介质振荡所产生的影响。研究表明板叠的抗对介质振荡 影响不大，阻的影响比较明显。 关键词：关键词：格子 boltzmann，谐振管，板叠，热声制冷机，振荡 华 中 科 技 大 学 硕 士 学 位 论 文 ii abstract thermoacoustic refrigerator is a type of environmentally friendly new refrigerator, its advantages make it have been researched all over the word. however, the working environmental of thermoacoustic refrigerator is non-linear coupling of acoustic and thermal, the space and time multi-scale effects and the flows with porous medium, which make it very difficult in theory research and numerical simulation. therefore, it is necessary to investigate new simulation analysis method. our team has been using lattice gas method to study the thermoacoustic engine and made some progress. the aim of present thesis work is to use the lattice boltzmann(lb) method to study the characteristics of medium oscillation in resonator. suitable model (d2q9 total energy distribution function lattice boltzmann model) was selected for the resonant. the characteristics of the media oscillation and heat transfer are studied by numerical simulation. suitable boundary condition for lattice boltzmann simulations also are selected for the physical boundary. at the same time, non-equilibrium extrapolation method is applied to dealing with isothermal boundary. the acoustic attenuation is simulated by d2q9 total energy distribution function lattice boltzmann model, which show the reliability of the lattice boltzmann model. the media oscillation characteristics in stable state of the resonant which was drived by a piston are analyzed. the effect of the driven frequency on the density amplitude is analyzed. the nonlinear oscillation phenomena of the resonator is analyzed preliminarily, the consistency of the results predicted with the experimental results shows the reliability of using lb to simulate media oscillation problem in resonator. finally, the resonance which has stacks in it is analyzed. the effects of the impedance of stacks to media oscillation are analyzed by comparing the discrepancy of phase angle and amplitude. at last, we concluded that the effects of the reactance can be neglected and the effects of the resistance is very obvious. keywords: lattice boltzmann, resonance, stack, thermoacoustic refrigerator, oscillation 独创性声明独创性声明 本人声明所呈交的学位论文是我个人在导师指导下进行的研究工作及取得的研 究成果。尽我所知，除文中已经标明引用的内容外，本论文不包含任何其他个人或 集体已经发表或撰写过的研究成果。对本文的研究做出贡献的个人和集体，均已在 文中以明确方式标明。本人完全意识到，本声明的法律结果由本人承担。 学位论文作者签名： 日期： 年 月 日 学位论文版权使用授权书学位论文版权使用授权书 本学位论文作者完全了解学校有关保留、使用学位论文的规定，即：学校有权 保留并向国家有关部门或机构送交论文的复印件和电子版，允许论文被查阅和借阅。 本人授权华中科技大学可以将本学位论文的全部或部分内容编入有关数据库进行检 索，可以采用影印、缩印或扫描等复制手段保存和汇编本学位论文。 保 密 在_年解密后适用本授权书 本论文属于 不保密 （请在以上方框内打“”） 学位论文作者签名： 指导教师签名： 日期： 年 月 日 日期： 年 月 日 华 中 科 技 大 学 硕 士 学 位 论 文 1 1 绪论 1.1 选题的背景及意义 热声制冷机是以热声效应作为工作原理，实现热和声的相互转换的一种新型环 保制冷机。热声制冷机具有诸多优点：结构简单、无污染、无或较少运动部件、无 润滑可靠性高以及可以利用太阳能之类的可再生能源及废热之类的低品位能源。由 于以上这些优点，在目前提倡节能和环保的大环境下，热声发动机和制冷机受到了 学术界和工业界的极大关注，得到广泛的应用研究。 但是由于热声转换的内部热力过程非常复杂，热声学科及热声问题的研究涉及 到热流体力学、声学、波动、振荡、非线性等许多分支学科的前沿领域和交叉学科 的基础研究。导致了实际热声制冷机的理论和模拟研究较为复杂。同时又由于热声 制冷机内部声场的非线性、热声的交叉耦合作用、热声元件的多尺度流动效应及 多孔介质结构。以上这些因素使热声理论的研究及热声装置系统的模拟具有一定的 挑战性。还没有完善的非线性热声理论和非线性热声模型求解的高效方法。传统的 数值微积分计算方法模拟热声的能力有限且具有局限性。因此，寻找高效的数值方 法和模拟技术，发展热声模型理论，对促进热声技术的应用和发展具有更重要的意 义。 随着电子信息时代的到来，高存储、高速度计算机的出现，人们开始利用高效 率的计算机来解决越来越复杂的物理问题。相比较，实验方法与高效计算机的出现 促进了计算流体力学的快速发展，并且在各个领域都得到了广泛的应用。数值模拟 方法一般分为两大类： （1）从宏观角度出发，用数值计算法求解基于连续介质假设的 navier-stokes 方程和描述热流体问题的偏微分方程62-64。 （2）从微观角度出发，采用分子动力学方法，基于动力学理论，通过跟踪介观 粒子的演化来实现对流动宏观现象的模拟。目前，该方法在复杂流体流动问题中得 到了极其广泛的应用65。 这两种方法在目前的科学研究中应用都比较多，两者在不同领域中发挥着各自 华 中 科 技 大 学 硕 士 学 位 论 文 2 独特的作用。但是，基于 navier-stokes 方程的传统方法（如有限体积、有限元、有 限差分法等）在模拟热声制冷机内复杂热流体问题时存在局限性：（1）边界问题复 杂，很难满足数值稳定性；（2）基本难以实现并行；（3）热声制冷机中非线性的 热声耦合机制及多孔介质元件内的多尺度流动效应理论不够成熟，较难用 navier-stokes 宏观方程实现对物理问题的合适描述。鉴于以上传热方法的局限性， 发展高效数值计算方法，利用其来研究热声制冷机及热声装置内部复杂流体流动和 传热问题是十分必要和有意义的。 格子 boltzmann 方法（lattice boltzmann method, lbm）是一种高效的流体模拟 的介观方法，作为一种微观方法，它起源于格子气自动机，同时也可以从连续 boltzmann 方程经过差分离散得到1,2。它与以连续微分方程为基础的宏观数值计算 方法有本质区别，lbm 是在空间、时间和速度上均离散的介观模型。与传统的计算 流体力学方法比较，lbm 具有很多不可取代的优势，其边界容易实现、具有完全并 行性、计算效率高及完全离散48。除此之外，格子方法打破了传统的建模观念，为 复杂系统的模型建立提供了新的途径。 由于 lbm 在研究复杂流动时有着许多其他计 算方法无法取代的优点，将 lbm 应用到热声谐振管介质振荡的研究中，并为推动该 方法在热声领域的应用做一些有益的探索与尝试，这是工作的主要研究目的。 1.2 研究现状 1.2.1 热声的模拟和研究 热声制冷机的主要结构包括谐振管、冷热端换热器及回热器。回热器是热声制 冷机的关键部件，常见热声制冷机的回热器一般都是各种不同形式的多孔介质。在 多孔介质内部实现热声效应，当声波在回热器中传播时，回热器的固体结构与气体 介质之间会产生能量形式的转换和转移，在回热器内部沿着或者逆着声波传播的方 向会产生一个热流，进而实现热流与声功的转换。此过程即热声效应。 在 19 世纪中叶人们相继发现 sondhauss 管3和 rijke 管4后，人们开始着手于热 声理论的研究。1868 年，feldman 和 kirchhoff 对管道中由于固体等温管壁和气体之 间的热传导引起声波衰减进行了研究3-5。 华 中 科 技 大 学 硕 士 学 位 论 文 3 直到 1887 年，lord rayleigh 提出 rayleigh 准则6，第一次对热声效应进行了合 理定性的解释。紧接着 1949 年，taconis 发现了著名的 taconis 振荡7，该发现对于 正确认识低温系统中振荡现象具有非常重要的意义。 日本的富永昭教授8从热动力学 观点出发较深入地研究了热声理论。 以上的所有研究都只是针对一维线性的热声现象和理论进行的研究，由于一维 线性理论有很大的局限性，只适合于压力振幅小，输入功率低时的情况。否则由于 高次谐波分量的增加、非线性效应增强会导致线性理论与实验结果偏离的情况发生。 随着发展和应用的需求不断提高，要求进一步提高热声热机的效率和功率，迫切需 要进一步研究非线性热声理论和大声幅热声现象。 swift 等人9首次对多孔介质中的湍流、非线性热声现象及声流现象的影响进行 了研究。紧接着，国外的大量学者采用 fluent、cfx、flow23d 等商用软件对 rijke 管内的热声振荡非线性现象进行了大量的数值研究。 美国 losalamos 国家实验室 cao n 等人通过求解可压缩流体的 navier stokes 方程对回热器板叠间隙间的热力过程 进行了数值模拟研究10。加拿大的 bauwens 采用多时间尺度法对一端封闭、一端输 入速度波的圆管内流动与换热进行了二维求解11,12。并且指出了采用一维方法得到 的一系列结果仅仅适合于管道大于热边界层的情况。美国 john hopkins 大学的 karpov 和 prosperetti 等在对热声制冷机中热声板叠周围低马赫数非稳态流动过程进 行了数值研究13-15之后，在此基础上他们建立了一套在截面平均假设条件下的准一 维近似热声声学方程16,17，该模型方程是弱非线性、时域表达的，同时适用于沿波 传播方向截面变化的复杂情况。接下来，他们又在原来的基础上提出一个更通用的、 时域的、非线性模型18-20。美国 texas 大学的 hamilton 等人建立了一套研究非线性 热声效应的二维模型和高效数值算法21,22。该模型的数值积分稳定条件与谐振管长 度相关，与板叠空隙直径无关。因此谐振管横截面可以随着位置的变化而变化，对 板叠长度与位置无限制。利用该模型，他们对起振过程和波形随谐振管形状变化的 非线性现象进行了数值研究。 国内的非线性研究，理论分析只局限在定性分析层次上，给出的大多是数值积 分模型。国内马大猷提出了一维非线性行波和驻波的严格解并利用该严格解对热声 华 中 科 技 大 学 硕 士 学 位 论 文 4 系统进行了多节推演。南京大学的韩非对 rijke 管内的非线性效应进行了深入研究 23,24，他们可以计算出管内声波从小扰动增长到有限幅值的振荡过程。西安交通大 学的刘继平等人在研究管内受热气体层流流动时发现了 rijke 管振荡的原因。 以上所有的非线性分析模型和数值研究，都可以看成是利用各种数值原理推导 出能描述非线性效应的偏微分控制方程，再利用各种数值方法对方程求解。这些模 型在不同程度上体现了热声系统的非线性本质特性，可以估计热声装置中板叠位置、 流固间温度梯度和横截面积变化等设计变量的变化趋势和影响。 但是这些模型和数值研究方法存在一些缺陷和局限性如算法复杂、模型控制方 程推导过程和数值积分求解过程复杂繁琐、边界条件难以处理、模型基于近似条件， 其解经常不适用，范围有限等。针对这些问题，近年来国内外有些学者开始尝试着 应用介观方法如格子方法来对热声问题进行研究。 1.2.2 格子方法应用于热声的研究 由于用传统的数值方法对热声进行研究存在诸多问题，而高效格子方法的发展 为研究非线性复杂热流动问题的研究提供了一种新的、有前景的途径25，其方法的 优势得到了广泛的关注，因此，近些年，一些学者提出使用格子方法来对热声发动 机和制冷机内的流动与换热过程进行研究。 格子气方法首次应用于热声装置的研究是清华大学的陈宇和刘旭以及华中科技 大学的张晓青等人26,27。他们用带温度的二维九速（d2q9）格子气模型对热声谐振 腔中的自激振荡及温度场的分布和随时间的变化过程进行了模拟。同时也研究了热 声回热器在声场中的位置及长度对声幅的影响。华中科技大学的张晓青及课题组等 人在前人的基础上，深入研究了热声振荡过程中粘性的影响、完善了边界处理方法 及格子气碰撞规则。研究了在不同热源温度下热声振荡的起振特性，对热声回热器 结构参数和空间位置及回热器两端温差随时间的响应进行了探索28,29,66。 同时他们开 发了新的热声格子气模型，将二维十三速格子气模型应用到热声装置的研究中，并 提出了相应的温度边界处理方法，开发了程序算法的实现技术30。 随着研究的深入，为了改善格子气带来的缺陷及保留原有优点，人们在格子气 华 中 科 技 大 学 硕 士 学 位 论 文 5 的基础上将格子 boltzmann 发展起来。较之格子气，格子 boltzmann 具有很多优点， 用实数分布函数取代整数布尔分布，消除了统计噪声，简化了碰撞算子的指数复杂 性同时稳定性增强。它在保留了格子气自动机的诸多优点的同时克服了前者的诸多 缺陷31,48。 格子 boltzmann 方法的应用涉及到很多领域，在微尺度流、量子力学、纳米流 体、声学光学等很多领域都有应用。在复杂边界和复杂流体的应用中取得了很大的 成功和突破，如多孔介质、悬浮颗粒流动、磁流体等48。 密西西比大学的 carl jensen 首次将热格子 boltzmann 模型应用到热声模拟中 32,33，针对多孔介质回热器的不同结构，研究多孔介质中热流动随时间的演化过程， 借助高性能计算机对不规则热声边界条件内多孔介质微尺度流动现象进行了模拟和 研究。 国内清华大学采用热格子 boltzmann 模型对共轭换热进行了研究34,35， 将热边界 滑移修正温度边界处理法应用到平行平板间的共轭换热物理模型中，分析了振荡流 共轭换热的速度场、温度场及热流分布与流固体导热系数比的关系。同时他们对多 孔介质中声波传播特性采用格子 boltzmann 进行了数值研究36，采用四参数方法生 成孔隙率已知的多孔介质，分析了多孔介质两端压差及振荡流的阻力系数与 renoylds 数之间的关系，同时分析了振荡压差与速度之间的相位差多孔介质空隙率 之间的关系。 西安交通大学的何雅玲，王勇等人采用热格子 boltzmann 方法对二位管道内不 可压的交变流动与传热进行了数值模拟37-39，分析了声波在管中的传播特性，得到 一些定性的有效结论。 所有以上可查阅的文献报道的研究，目前都只是在热声格子方法及模拟研究方 面的探索，其研究结果是鼓舞人的，但研究工作由于格子方法在热声领域的研究具 有较大的开创性，也具有一定的挑战性，还有待在广度和深度上进一步探索。 1.2.3 谐振腔的格子方法研究 谐振腔是热声制冷机的一个基本，重要组成部分，提供声振荡的条件，因而大 华 中 科 技 大 学 硕 士 学 位 论 文 6 量学者对声波在管中的传播进行了大量实验和数值模拟计算研究来掌握谐振腔的声 特征。 英国的 buick40采用 d2q7 格子 boltzmann 模型对声波在无限长管的衰减进行 了数值模拟。a. alexeev41对封闭腔中气体的周期性振荡进行实验和数值模拟研究， 实现了谐振腔中压力和温度的可视化，发现管内的气体温度随时间发生明显的变化。 结果显示在谐振腔体中气体的流动趋势与无限长管中的流动趋势相反。文中提出一 个研究振荡流的重要的参数（l l=，是一个很小的参数），振荡流中的热动 力学变化和速度的变化都是该参数的函数，一维气体的振荡受该参数很大的影响。 ilgamov et al 42对封闭腔的气体振荡进行大量的数值计算和实验研究， 他们发现随着 参数的增大， 实验数值与理论值的偏差增大。 通过大量的实验得出结论当0.1时， 现有的声学理论就不太适合应用在谐振腔里面了。 近年来国内不少学者也对声波在谐振腔中的振荡进行了数值计算研究，西安交 通大学的何雅林，王勇等人用格子 boltzmann 对谐振腔进行了数值模拟。模拟结果 显示用格子 boltzmann 模拟的结果与以往学者的实验数值以及传统模拟结果都吻合 良好37-39。 1.3 主要内容 由于热声问题的复杂性，存在非线性理论不完善及解析解难以求得的弊端，采 用传统的数值计算方法和实验研究都存在较大的局限性。因而作者将格子方法应用 到热声的研究中。格子方法存在分子动力学的优点，边界易于处理、具有完全并行 性及编程容易实现。从理论上分析可以得出，格子方法避免了传统数值方法的诸多 缺陷，将其应用到热声的研究中很有发展前景。在本论文中，将采用 d2q9 格子 boltzmann 模型对热声谐振腔中介质振荡特性进行研究。 第一章绪论部分介绍了热声效应的基本理论，热声效应的发展历史和研究现状， 重点概述了格子方法应用于热声的研究现状，以及热声热机中谐振腔的研究状况。 同时通过了解历史和分析现状，阐述论文工作的意义所在。 第二章从驻波场热声制冷角度出发，通过分析热声制冷基本原理简单介绍热声 效应。 同时介绍格子boltzmann的发展历史和应用及几个典型的bgk格子boltzmann 华 中 科 技 大 学 硕 士 学 位 论 文 7 模型。 第三章根据谐振腔物理模型， 选择适合该物理问题的 d2q9 格子 boltzmann 模型 （总能双分布热格子 boltzmann 模型）及边界处理方式。同时将非平衡态外推法应 用到等温边界处理中，并应用格子 boltzmann 方法对谐振腔的模拟流程进行描述。 第四章为谐振腔的格子 boltzmann 模拟和分析。首先对声波在无限长管道内的 传播进行研究。分析声波的速度和压力与位置，粘性及频率之间的关系，并解模拟 结果与理论分析解进行比较，对格子 boltzmann 方法及程序进行验证。其次对声波 在封闭的谐振腔体中的传播特性进行研究，分析谐振腔中介质的振动特征，以及压 力、速度随位置的变化关系。并对谐振频率下声波的非线性积累现象进行重点研究。 最后研究谐振腔中加入板叠对声波传播的影响，分析板叠孔隙率对声波速度和压力 的影响及其相互关系。 第五章对全文进行总结，并对未来的研究工作进行展望。 华 中 科 技 大 学 硕 士 学 位 论 文 8 2 热声基本原理及格子 boltzmann 理论基础 由于热声学是热学、声学和流体力学等多门学科的交叉学科，属于前沿科学。 包含了声学、热声和流体力学的高深理论知识；格子 boltzmann 方法是新兴的高效 数值算法，包含了统计物理和计算数学的各方面的知识，理论基础比较深厚。在开 展后续工作之前对热声的基本理论及格子 boltzmann 方法的理论基础给予简单介绍。 2.1 热声基本原理 分析热声系统中的热力循环可以看出，热力学循环与声场有关。不同声场中的 热力学循环过程不同。结合文章的主要工作，下面从驻波场的声制冷角度来分析热 声效应基本原理。 2.1.1 声波特性 首先分析声波特性，由于任何线性声波都可以表示成简谐波的叠加，我们可以 通过分析简谐声场来得到声波的特性。 频率相同传播方向相反的两列简谐声波沿着x 轴传播，前行波与反行波的声波压强波动量可以用复数分别表示成43： ()it kx iia pp e =， ()it kx rra pp e + = （2.1） 两列声波叠加以后得到： () 2cos() i tit kx irraiara ppppkxeppe = +=+ （2.2） 其中， p 表示声场压强的波动量。式（2.2）中右端第一项为驻波场，驻波场每个位 置的媒质均做同相位振动，振幅大小随位置发生变化。第二项为平面行波场。质点 振动速度为： 1p udt x = （2.3） 当式（2.2）中两列声波振幅相同时，第二项为零，叠加声场为纯驻波场。令 0 p 为静 压， iaraa ppp=得到驻波场声压压强为 0 2cos i t a pppkxe =+ （2.4） 华 中 科 技 大 学 硕 士 学 位 论 文 9 令流体静止密度为 0 ，得到声场压强与振动速度之间的关系为： 0 pu xt = （2.5） 结合（2.4）及（2.5）可得到驻波场流体振动速度为： 00 21 sin i t a ip kp udtkxe x = = （2.6） 位移为： ti a kxe kp udts sin 2 0 2 = （2.7） 对式（2.4），（2.6），（2.7）取实部，可以得到流体压强p，速度u及位移s随时 间t的变化关系。其中压强p与位移s之间相差180相角；速度u比压强p落后90 相角。 2.1.2 热力循环过程 热声热机中的热力循环过程一般有三种情况：（1）绝热循环过程；（2）流体 微元与固体板叠等温循环过程；（3）一般的热力循环过程，即固体壁面与外热源间 的有限热交换过程67。其中（1）和（2）两个循环过程分别表示流体微元与板叠及 壁面间不进行热量交换，以及流体微元与板叠瞬时进行热力交换的过程，也就是导 热系数为分别0和的两种理想情况。实际中的热力循环过程都是处于两者之间的 一个过程，导热系数为有限正实数。 在加了板叠的soundhauss管的板叠两端加上两个温度不同的热源 （如图2.1） 时， 发现可以获得显著的热声效应，令高温端热源温度为 h t，低温端热源温度为 c t，管 中板叠长度为 1 l。考虑板叠处的温度梯度为() 1hc ttl线性的情况。 华 中 科 技 大 学 硕 士 学 位 论 文 10 图 2.1 加入板叠的 soundhauss 管 先对绝热循环过程进行分析，当流体振动过程是一个绝热的过程时，则根据热 力学可以得到流体的温度表达式： 1 tcp = （2.8） 在小振幅振动下，式（2.8）可以化简为： 11 0 00 00 2cos11 (1)(1) i t a pppkx tcpcpe pp =+=+ （2.9） 根据式（2.7）和式（2.9）可以得到绝热循环过程下流体温度与位移的关系表达式： 2 00 0 0 1t ttctgkx s p k = （2.10） 其中x表示水平坐标。将式（2.10）画在图上，可以得到如图2.2所示的直线。当 0 0ctgkx时，直线为负，反之则为正。 由热力学知识可知，绝热过程中 1 tvc =，根据式（2.10），当 0 0ctgkx时， 气体微团体积随位移的增大而增大，反之则变小。由于绝热过程中 1 tvc =，该曲 线在pv图中包围的面积为0 ，得出绝热过程不做功的结论。 接下来对等温过程进行分析。当流体微元与板叠等温时，温度与位移的关系表 达式为： () 0211 tts ttl= （2.11） 将式（2.11）在温度位移关系图上表示时，可以得到直线。为实际循环过程中温 度与位移的关系。 华 中 科 技 大 学 硕 士 学 位 论 文 11 2.1.3 声制冷效应 如图2.2所示，图中 0 x表示气体微团平衡位置； 0 t表示静止状态下气体微团温 度。 从图中可以看出来，1状态点处为温度最大值和负向位移最大值， 体积最小值处； 2，4状态点的板叠温度与气体温度一致；3状态点为气体微团位移最大值，温度最 小值和体积最大值处。1点附近处的板叠温度低于气体微团温度，气体微团放热。当 气体微团从1点往平衡位置点运动时，温度降低，对外做功；当气体微团从平衡位 置往1点运动时，温度升高，外界对其做功。由能量守恒定理可以得到，气体微团 从平衡位置向状态点1运动时，温度随位移增大而升高的值要低于气体微团从1点 往平衡位置运动温度降低的值。从而整个循环方向为逆向的。 0 t 0 xs 0 x 0 xs+ 图 2.2 声制冷中温度与位移的关系 如图2.3所示，当在412 过程中时，气体微团向板叠放热；在234过 程中，气体微团从板叠吸收热量；在123过程中，气体微团体积膨胀，对外做 功；在341过程中，气体微团体积变小，外界对其做功。 如图2.4所示，为声制冷循环过程中的pv图。和为等温线，分别表示实 际循环中，温度的最大值和最小值；为实际循环过程中压强与体积的关系曲线。1 状态点处，流体微元的温度最大，体积最小，气体微元向板叠放热；3状态点处，气 体微元体积最大，温度最小，气体微元从板叠处吸收热量。 华 中 科 技 大 学 硕 士 学 位 论 文 12 图 2.3 声制冷循环效应 从声制冷热力效应的整个过程看来，气体微团从低温热源吸收热量，往高温热 源放出热量，并且从图2.4可以看出外界对气体微团做功。在1状态点气体微团被压 缩到最小并向板叠的高温端放热；在3状态点气体微团膨胀最大前从板叠低温端吸 收热量，根据瑞利准则，系统声振动会衰减，需要外界不断补充声功维持声振荡。 该循环反映了声制冷的基本原理。 1 3 p v 图 2.4 声制冷过程中压强与体积的关系 华 中 科 技 大 学 硕 士 学 位 论 文 13 2.2 格子 boltzmann 理论基础 流体是自然界最常见的一种物质。自然界不同流体之间形态各异，各有特色， 但追根溯源却发现流体的各种宏观现象都是由于分子之间的热运动频繁碰撞所造成 的。在物理上，流体可以看成是由大量分子所构成的一个离散系统，分子之间的热 量和动量交换都是通过分子之间的热运动达到的。从微观角度来看流体，它是一个 时间和空间非常复杂的系统，特性为：随机性、离散性和不均匀性。而从宏观角度 来看，则是一个可以用宏观特征量如密度、压力、速度、温度及内能来描述的系统， 具有确定性、连续性和均匀性。流体的宏观运动及其宏观特征量都是在微观运动统 计平均下的结果48。 在流体的连续介质假设条件下，流体的宏观特征和运动可以用非线性偏微分方 程描述出来，即包含流体质量、动量和能量守恒的navier-stokes方程62-64： ().0u t += （2.12 a） () ()(). 3 u uupuug t += + （2.12 b） () ()(). e uetpu t += + （2.12 c） 其中，、t、u、e分别表示密度、温度 、速度和内能，压力表示prt=， 动 力 学 粘 性 系 数=，为 运 动 粘 性 系 数 ，为 导 热 系 数 ， () 2 2:2.3suu = 为粘性热耗散函数。 上述方程成立的条件是流体为连续介质，适用于所有的连续流体。一般的情况 下，大多数的流体都属于连续介质，因此，流体的宏观运动都能够用navier-stokes 方程正确的表达出来。对于不同的物理问题，可以对（2.12）系列方程进行变化，进 一步得到euler方程、boussinesq方程或者其他粘性不可压方程。但是在一些复杂流 体流动中，如多相流、多孔介质内的复杂流动及多组分系统中要求解该偏微分方程 很困难或者说在目前的基础上无法精确求解，只有给与一定的假设条件，求其近似 解。同时，有些更为复杂的流动系统，我们无法得到其相应的数学模型。对于这些 华 中 科 技 大 学 硕 士 学 位 论 文 14 复杂流动系统，使用分子动力学方法或者介观方法进行描述更为恰当和可行。 分子动力学方法是一种微观的模型，从流体的微观结构出发研究分子运动，通 过跟踪分子的碰撞演化，统计各个分子信息来得到流体的宏观特征量。但是该方法 需要跟踪多个粒子，计算量和存储量都很大。因此到目前为止都只能用该方法对简 单的物理问题进行研究。复杂流动问题由于受到庞大计算量的限制而无法涉及。分 子动力学理论中的boltzmann方程是一个相对于navier-stokes方程更为复杂的方程， 无法用常规的方法来实现它的完全求解。在这种情况下，格子boltzmann方法成了 一种求解复杂流动问题有效的方法和途径48。 2.2.1 格子boltzmann方法及其发展历史 格子boltzmann方法可以从格子气自动机(lattice gas automata,lga)发展而来也 可以从统计力学的boltzmann方程离散得到。 但是最初的格子boltzmann方法是由格 子气自动机发展而来的。格子气自动机方法是借鉴离散运动论的分子运动和元胞自 动机演化过程提出来的44,45,48，是计算流体力学的一种新发展。 lga是一个完全离散的动力学系统：在空间上，整个流场离散成有限个规则格 点； 时间上， 根据时间步长离散为一时间序列； 流体离散为大量的介观微粒。 在lga 中微粒只能沿着网格线运动，并且只能在一个时间步内从一个格点线性迁移到相邻 的格点48。 20世纪70年代，hardy,de pazzis和pomeau提出第一个lga模型 46,47 。该模 型为一个二维正方形格子lga模型（如图2.5（1）所示），称之为hpp模型。在 模型中，每个格点状态由四个布尔变量表示：() 1234 ,nn n n n=。hpp模型能够满足 质量和动量守恒但是无法满足navier-stokes方程，因为其没有足够的对称性。 为了得到正确的宏观动力学方法，需要提出更高对称性的格子模型。1986年 frish,hasslacher和pomeau（fhp）提出了一个格子气模型（fhp模型）45，该模型 为一个二维的正六边形模型（如图2.5中的（2）所示）。fhp模型能成功恢复到 navier-stokes方程，为开创lga的研究打开了新局面。 华 中 科 技 大 学 硕 士 学 位 论 文 15 1 2 3 4 (1) 图 2.5 hpp 模型格点和 fhp 模型格点 但是格子气自动机存在不足：（1）存在统计噪声；（2）碰撞算子存在指数复 杂性；（3）不满足伽利略不变性。格子气自动机的统计噪声可以通过采用时间或者 空间平均法加以克服，但是这对电脑的计算和存储有很高的要求。为了消除噪声， 1988年，mcnamara和zanetti提出一种想法，即用粒子分布函数来替代布尔变量进 行演化，从而开创了流体计算的一个新的方向，格子boltzmann方法48。 mcnamara和zanetti提出的这个模型（简称mz模型）即最早的格子boltzmann 模型。在mz模型中，粒子分布函数按照与格子气自动机类似的方式来演化： ()()()(), iiiii fxct ttfx tfx t+=+ （2.13） 其中 ii cce=， i cxt= 为离散速度， i e为离散速度单位矢量， i 为碰撞算子。该模 型中 i 与格子气自动机的碰撞算子形式一致，但用粒子分布函数替代了布尔变量。 由于布尔型的粒子数直接用实数型的粒子分布函数替代，mz模型从本质上克服了格 子气自动机方法的噪声。但是该模型采用的碰撞算子任然具有指数复杂性，格子气 自动机的其他缺点它不能克服48。 higuera和jimenez在1989年对mz模型进行了改进，得到了一个新的模型（简 称hj模型）。他们将碰撞算子进行简化，用一个线性算子近似。引入一个平衡态分 布函数 eq i f，同时假设： eqneq iii fff=+ （2.14） 其中 neq i f表示非平衡态部分，量级为( )（是chapmann-enskog多尺度展开中的 小参数） ，该模型的碰撞方式与格子气自动机的碰撞方式一致，则 eq i f与一般的格子 华 中 科 技 大 学 硕 士 学 位 论 文 16 气自动机模型的平衡态分布函数一致。对碰撞算子 i 进行泰勒展开，并略去和马 赫数（ma）的高阶项，可以得到： () ( ) () 0neqeq iiijjijjj ffaff= = （2.15） 其中 ( )0 ijij a = 为碰撞矩阵。 hj模型的碰撞矩阵和平衡态分布函数虽然任然依赖原格子气自动机模型，但其 计算复杂度与mz模型的指数形式相比被大大降低，也提高了计算效率。同时，hj 模型首次在演化过程中提出了平衡态分布函数，为格子boltzmann方法后来的发展 打好了基础。但是hj模型却有一个不足，即数值稳定性不好。同时碰撞算子与格子 气自动机的碰撞规则有关，且平衡态分布函数本质上是femi-dirac分布，这些特点 使得适用范围大大受到限制48。 为了使得格子boltzmann模型更具有通用性，higuera、succi和benzi构造了一 种新的格子boltzmann模型（简称hsb模型）。该模型利用maxwell-boltzmann平 衡态分布函数来替代格子气自动机模型中的femi-dirac平衡态分布函数，同时设计 适合的碰撞矩阵，使得相应的宏观方程得到恢复。hsb模型的碰撞算子为many矩 阵，为了更简化，在1991年前后学者们提出了一种单松弛模型49,50。 ()() 1 eq iiii fff = （2.16） 其中 eq i f表示平衡态分布函数，为无量纲松弛时间，即粒子分布函数从非平衡态达 到平衡态所需要的时间。 该模型也称为lbgk（lattice bhatnagar-gross-krook） 模型。 在lbgk模型中，最关键的地方是平衡态分布函数 eq i f的确定，选择合适的 eq i f，可 以恢复到相应的navier-stokes方程。一般的，平衡态分布函数 eq i f，需要满足如下的 质量和动量守恒条件： eq ii ii ff= ， eq iiii ii uc fc f= （2.17） lbgk模型大大地提高了计算效率， 并且成为目前格子boltzmann方法中使用最广泛 的lbm模型。 华 中 科 技 大 学 硕 士 学 位 论 文 17 2.2.2 几个典型的格子boltzmann模型介绍 bgk模型是目前应用最为广泛的一种模型，也是格子boltzmann方法的主要模 型。下面要介绍的是几种典型的bgk格子boltzmann模型，针对几个典型的物理问 题来简单介绍一下二维九速（d2q9）的不可压等温格子boltzmann模型和非耦合的 双分布热格子boltzmann模型。 论文中使用的格子boltzmann模型将在第三章重点介 绍。 （1）不可压等温格子boltzmann模型48 郭照立在前人的基础上构造了一个能够模拟一般不可压navier-stokes方程的二 维九速格子boltzmann模型51 。该模型是一个不可压模型，平衡态分布函数满足： eq i i fconst= 。模型的演化方程为： ()()()() 1 , eq iiiii fxct ttfx tfx tfx t += （2.18） 其中(), i fx t表示粒子密度分布函数，为松弛时间， i c表示粒子的离散速度 （cxt= ，x为空间步长，t为时间步长）。 平衡态分布函数定义为： () ( ) ( ) ( ) 00 2 2 2 4,0 ,1 4 ,58 eq ii i p sui c p fx ts ui c p s ui c += =+= += （2.19） 其中，(), eq i fx t表示平衡态分布函数； 0 为静止时流体的密度；p为流体压力；、 和为模型参数，需要满足条件：+=及21 2+=；( ) i s u的表达式为： ( ) () 2 2 242 . 22 i i ii sss c uc uu s u ccc =+ （2.20） 其中， i 为权系数，对于d2q9格子模型，权系数分别为 0 4 9=， 1 4 1 9 =， 5 8 1 36 =； s c为声速（ 22 3 s cc=，其中c为粒子速度）；u为宏观流速