下载文档原格式
Fluent软件的燃烧模型介绍Fluent软件中包含多种燃烧模型、辐射模型及与燃烧相关的湍流模型,适用于各种复杂情况下的燃烧问题,包括固体火箭发动机和液体火箭发动机中的燃烧过程、燃气轮机中的燃烧室、民用锅炉、工业熔炉及加热器等。
燃烧模型是FLUENT软件优于其它CFD软件的最主要的特征之一。
下面对Fluent软件的燃烧模型作一简单介绍:一、气相燃烧模型·有限速率模型这种模型求解反应物和生成物输运组分方程,并由用户来定义化学反应机理。
反应率作为源项在组分输运方程中通过阿累纽斯方程或涡耗散模型。
有限速率模型适用于预混燃烧、局部预混燃烧和非预混燃烧。
应用领域:该模型可以模拟大多数气相燃烧问题,在航空航天领域的燃烧计算中有广泛的应用。
PDF模型该模型不求解单个组分输运方程,但求解混合组分分布的输运方程。
各组分浓度由混合组分分布求得。
PDF模型尤其适合于湍流扩散火焰的模拟和类似的反应过程。
在该模型中,用概率密度函数PDF来考虑湍流效应。
该模型不要求用户显式地定义反应机理,而是通过火焰面方法(即混即燃模型或化学平衡计算来处理,因此比有限速率模型有更多的优势。
应用领域:该模型应用于非预混燃烧(湍流扩散火焰,可以用来计算航空发动机的环形燃烧室中的燃烧问题及液体/固体火箭发动机中的复杂燃烧问题。
非平衡反应模型层流火焰模型是混合组分/PDF模型的进一步发展,从而用来模拟非平衡火焰燃烧。
在模拟富油一侧的火焰时,典型的平衡火焰假设失效。
该模型可以模拟形成Nox的中间产物。
应用领域:该模型可以模拟火箭发动机的燃烧问题和RAMJET及SCRAMJET 的燃烧问题。
预混燃烧模型该模型专用于燃烧系统或纯预混的反应系统。
在此类问题中,充分混合的反应物和反应产物被火焰面隔开。
通过求解反应过程变量来预测火焰面的位置。
湍流效应可以通过层流和湍流火焰速度的关系来考虑。
应用领域:该模型可以用来模拟飞机加力燃烧室中的复杂流场模拟、气轮机、天然气燃炉等。
译者:wyxpuma时间有限,错漏之处在所难免.如有发现可站内 message 我: )使用非预混燃烧模型 Using the Non-Premixed Combustion Model 概要 先决条件 问题描述 使用 prePDF 的准备工作 Step 1: Define the Preliminary Adiabatic System in prePDF Step 2: Compute and Review the Adiabatic System prePDF Look-Up Tables Step 3: Create and Compute the Non-Adiabatic prePDF System Preparation for FLUENT Calculation Step 4: Grid Step 5: Models: Continuous (Gas) Phase Step 6: Models: Discrete Phase Step 7: Materials: Continuous Phase Step 8: Materials: Discrete Phase Step 9: Boundary Conditions Step 10: Solution Step 11: Postprocessing Step 12: Energy Balances and Particle Reporting SummaryIntroduction 概要对煤粉燃烧的模拟包括对连续气相的模拟及其对煤颗粒分散相的相互作用 的模拟.煤颗粒穿过气相时,会进行脱挥发分和焦炭的燃烧,从而成为进行气相 反应的燃料的源.对反应的模拟可通过组分传输模型(the species transport model)或非预混燃烧模型(the non-premixed combustion model)进行模拟. 在此实例指南中, 你将使用非预混燃烧模型模拟一个简化的燃煤炉体内的化学反 应. 在此实例指南中,你将学习: 使用 prePDF 预处理器为煤粉燃料(的燃烧)准备一个 PDF 表 为非预混燃烧化学反应模型定义 FLUENT 的输入参数 为煤颗粒定义一个离散相 求解一个包含反应离散相的煤颗粒的数值模型非预混燃烧模型是采用这样一种建模方法, 即求解一个或两个守恒标量以及混合 分数的的传输方程.多化学组分时,包括根和中间产物,可以采用这种方法被考译者:wyxpuma时间有限,错漏之处在所难免.如有发现可站内 message 我: )虑进来,并且他们的浓度能够通过预测出的混合分数分布计算出来.通过化学数 据库,可以得出组分的正确数据.湍流和化学反应的相互作用可通过使用β或双 δ的概率密度函数 (PDF) 来模拟. 关于非预混燃烧的模型详见 the User's GuidePrerequisitesThis tutorial assumes that you are familiar with the menu structure in FLUENT, and that you have solved Tutorial 1 or its equivalent. Some steps in the setup and solution procedure will not be shown explicitly.问题描述(Problem Description )本例程考虑的煤粉燃烧系统时一个简单的二维管道,截面为 10m×1m,如图 13.1 所示.因为管道对称,图中只显示了宽度的一半.该二维管道的入口分成 两股来流.靠近管道中心的高速来流速度为 50m/s,跨度为 0.125m.另一股来流 速度为 15m/s,跨度为 0.375m.两股来流都是 1500K 的空气.靠近中心高速流股 进入炉膛的煤颗粒,其质量流率为 0.1kg/s(煤颗粒的总流率为 0.2kg/s).管 道的壁温为 1200K 的常温. 基于入口尺寸和平均入口速度的雷诺数约为 100,000. 因此,流动为湍流. 关于煤组成和粒度分布的详细信息将在 step 5(Models: Continuous (Gas) Phase)和 Step 8( Materials: Discrete Phase)中给出.Figure 13.1: 2D Furnace with Pulverized Coal Combustion译者:wyxpuma时间有限,错漏之处在所难免.如有发现可站内 message 我: )对 prePDF 的准备(Preparation for prePDF )1. 启动 prePDF.当你使用非预混燃烧模型,你要先使用前处理器—prePDF-生成一个 PDF 文件. 该 PDF 文件包含了关于各组分浓度,温度与混合分数值之间的联系的信息,并被 FLUENT 在求解过程中用来获得这些标量(组分浓度,温度等) .Step 1: 在 prePDF 中定义初步的绝热系统 (Define the Preliminary Adiabatic System in prePDF )1. 定义 prePDF 模型类别.你可以定义单一燃料流或者燃料流加上另一个流股.激活第二个流股可以让你明了 两个混合分数.对于煤颗粒的燃烧,采用两个流股的方法可以让你明确地将挥发分 (第二个流股)流股同焦炭流股(燃料流)分开.在该例程中,我们将不采用两个 流股的方法,而采用单一混合分数法.SetupCase...(a)译者:wyxpuma时间有限,错漏之处在所难免.如有发现可站内 message 我: )在 Heat transfer options, 保持默认设置 Adiabatic.该例程中讨论的煤粉燃烧器是一个非绝热系统, 在气相与煤颗粒相间, 以及燃 烧器壁面上,都有热量传递.因此,必须在 prePDF 中考虑非 绝热系统.由于非绝热系统较之绝热系统在计算时要耗费更多的时间,因此你 可在启动 prePDF 时先考虑系统为绝热系统. 根据对绝热系统的 PDF/化学 平衡的计算结果,你将确定大致的系统参数,这将使得在非绝热系统下的 计算更具效率.而且,绝热系统的计算将提供绝热火焰温度,混合分数的 化学当量以及各单独组分对化学反应系统的重要性的信息. 这种首先在绝 热系统中计算的方法将被用于所有的考虑在非绝热条件下的 PDF 计算.(b) 在 Chemistry models, 保持默认设置 Equilibrium Chemistry. 对大多数基于 PDF 的数值模拟,推荐 使用 Equilibrium Chemistry 选项. Stoichiometric Reaction ( 混合即完全反应 ) 选项计算量小,但不准确; Laminar Flamelets 选项能够考虑空气拉伸的非平衡现象,例如超平衡基浓度和亚平衡温度. 这对于 NOx 的预测非常重要,但在此例中不必考虑. (c) 保持默认设置 PDF models. 一般的,推荐使用β函数的 PDF 积分,因为他比δ函数的 PDF 积分方法更准确.(d) 在 Empirically Defined Streams, 选择 Fuel stream 选项. 这将允许你采用 empirical 输入选项来定义燃料流股(即根据煤的 ultimate analysis 定义) Empirical 输入选项允许你指定燃料中 H,C,N 和 O 的元素组成 (DAF) 低位发热量和比热. , 当元素分析和热值已知时, 这种指定方法是很有效的. (e) 点击 Apply ,关闭面板 2. 定义系统中的化学组分指定系统中包含哪些组分依赖于燃料类型和燃烧系统. 关于指定组分的指南在FLUENT User's Guide 中有介绍.在这里,你将假定该平衡系统中包含 13 种组分: C, C(s), CH , CO, CO , H, H ,H O, N, N , O, O , and OH.之所以包含 C, H, O, 和 N 是因为我们采用``empirical''输入方法 对燃料进行定义的,因此要用到元素组分.! 当使用 empirical input 选项时,你应当将 C and C(S) 都选入系统SetupSpeciesDefine...译者:wyxpuma时间有限,错漏之处在所难免.如有发现可站内 message 我: )(a) 将 Maximum # of Species 设为 13. 使用上,下箭头来设置最大组分数目,或在数 字框中直接键入数字然后摁<ENTER>. (b) 在 Defined Species 列表中选择最上面的组分(初始标志为 UNDEFINED). (c) 在 Database Species 下拉表单中,使用滚动条拖动表单,并选择 C. 此时 Defined Species 表单中的最顶部组分将显示为 C . (d) 在 Defined Species 列表中选择第二个组分 (或在 Species # 选项中将数目增大为 2). (e) 在 Database Species 下拉 表 单 中 , 使 用 滚 动 条 拖 动 表 单 , 并 选 择 下 一 个 组 分 ( C(S)). (f) 重复 步骤 (d) 和 (e),直到 13 个组分都被定义 (g) 点击 Apply,关闭面板 Note:在其他燃烧系统中,你也许想添加其他化学组分,但你不要添加慢反应的化学组分, 例如 NOx..3. 确定燃料组分输入.本例中考虑的燃料是已知的,根据工业分析,知道含 28%的挥发分,64%的焦炭和 8%的灰分.你要根据这些信息,以及以下给出的元素分析结果,在 prePDF 中定义 煤的组成.燃料流的组成(焦炭和挥发分)可计算如下:译者:wyxpuma时间有限,错漏之处在所难免.如有发现可站内 message 我: )首先将工业分析数据(干燥基)转化为可燃基(干燥无灰基).工业分析 Wt % Wt % Proximate Analysis (dry) Volatiles Char (C(s)) Ash 28 64 8 (DAF) 30.4 69.6 -煤的元素分析如下表:Element Wt % (DAF) C H O N S 89.3 5.0 3.4 1.5 0.8为简单起见,煤中的硫分可并入氮的质量分数,因此最终数据如下:Element Wt % (DAF) C H O N S 89.3 5.0 3.4 2.3 -我们可根据工业分析和元素分析的数据得出挥发分的元素组成如下: can combine Wethe proximate and ultimate analysis data to yield the following elemental composition of the volatile stream: (译者注:这里不是挥发分流股的元素组成,而是整个燃料流股的元素组成,因为该模型采用单混合分数模型,而且表中数据也显示是单流股的摩 尔组成)译者:wyxpuma时间有限,错漏之处在所难免.如有发现可站内 message 我: )Element Wt % Moles Mole Fraction C H O N Total 89.3 5.0 3.4 2.3 7.44 5 0.21 0.16 12.81 0.581 0.390 0.016 0.013你将用到上表中最后一列中的数据.prePDF 将使用这些信息以及煤的热值来定义 燃料中的组分. 煤的低位热值为(DAF) : LCV= 35.3 MJ/kg煤的热容和密度分别为 1000 J/kg-K 和 1 kg/m4. 输入燃料和氧化剂的组分SetupSpeciesComposition...(a) 激活氧化剂流股组分的输入框氧化剂(空气)的组成为 21%的 O2 和 79% 的 N2译者:wyxpuma时间有限,错漏之处在所难免.如有发现可站内 message 我: )i. 在 Stream 下, 选择 Oxidiser. ii. 在 Specify Composition In, 保持默认选项 Mole Fractions. iii. 在 Defined Species 选项中选择 O2 并在 Species Fraction 中输入 0.21. iv. 在 Defined Species 选项中选择 N2 并在 Species Fraction 中输入 0.79 (b) 激活燃料流股组分的输入框 Note: 因为燃料流股已经选择了 empirical 输入选项,您将被提示输入 C,H,O 和 N 的原子摩尔分数,以及燃料的热值和热容.译者:wyxpuma时间有限,错漏之处在所难免.如有发现可站内 message 我: )i. 在 Stream 下, 选择 Fuel. ii. 在 Specify Composition In 中,保持默认选项 Mole Fractions. iii. 在 Defined Species 列表中选择 C 并在 Atom Fraction 中键入 0.581. iv. 在 Defined Species 列表中选择 H 并在 Atom Fraction 中键入 0.390. v. 在 Defined Species 列表中选择 N 并在 Atom Fraction 中键入 0.016. vi. 在 Defined Species 列表中选择 O 并在 Atom Fraction 中键入 0.013. vii. 在 Lower Caloric Value 中键入 3.53e+07 J/kg 在 Specific Heat 中键入 1000 J/kg-K . viii. 点击 Apply ,关闭面板. 5. 定义固体碳的密度.译者:wyxpuma时间有限,错漏之处在所难免.如有发现可站内 message 我: )这里,假定固定碳密度为 1300 kg/m3SetupSpeciesDensity...(a) 在 Defined Species 中选择 C(S) . (b) 将 Density 设为 1300. (c) 点击 Apply ,关闭面板. Note:在计算燃料的混合密度时,prePDF 会用到这个信息.你应当定义固定焦炭的密度, 这个值可能和 FLUENT 中定义的煤的密度不一致, FLUENT 中定义的是包含灰 的煤粒的显密度.6. 定义系统操作条件.在化学平衡计算中需要用到系统压力和入口流温度. 煤燃烧情况下, 燃料流股 的入口温度应当是开始脱挥发分的温度.氧化剂流股的入口温度应当对应于空气入 口温度.在本例程中,煤脱挥发分的温度设置为 400 K,空气入口温度为 1500K,系 统压力为 1 大气压.SetupOperating Conditions...(a)在Fuel和Oxidiser中分别将温度设为400 K 和1500 K。
2005 Fluent 中国用户大会论文集采用FLUENT软件研究旋流煤粉燃烧器燃烧特性由长福（清华大学热能工程系，北京 100084）摘要：本文FLUENT软件研究了实际电站锅炉单个双调风旋流燃烧器附近区域的煤粉燃烧过程。
并分别研究了内二次风旋流强度，外二次风风率，一次风风率和三次风风率等因素对燃烧性能的影响。
各工况计算结果表明，总体上在燃烧器出口处形成了高温区和高煤粉浓度区，燃烧器出口一定距离后的炉内温度呈逐渐上升趋势，炉膛温度分布均匀。
中心高温区出现迟的工况，后期分级燃烧充分。
表明该燃烧器具有高效稳燃和变工况运行稳定的性能。
关键词：旋流燃烧器；数值计算；燃烧性能引 言当前国内使用的电站锅炉，80%是四角切圆煤粉燃烧锅炉，不到10%采用旋流燃烧锅炉[1]。
和四角切圆煤粉锅炉相比，旋流燃烧器锅炉是一种新型的锅炉，结构复杂得多。
已有较多学者采用数值模拟方法研究旋流燃烧器燃烧性能的例子[1-4]，这些例子的计算结果都详细预报了由于测量困难而不能充分获得的炉膛内部的温度场，速度场，燃烧产物各组分的浓度分布和污染物的分布，其中文献[2]和[3]还与实验数据比较，比较结果表明，模拟结果与锅炉热态试验数据吻合情况较好，为数值模拟的更广应用提供了依据。
简图如图1燃烧器中心通一股直流的三次风，风量较小。
针对该燃烧器的结构，本文研究了内二次风的旋流强度，二次风的配比，一次风和中心风的风率对燃烧性能的影响。
作者：由长福（1969），男（汉族），黑龙江，副教授，博士，清华大学热能工程系1 计算方法1.1 计算对象和网格生成计算域为单个旋流燃烧器附近的区域，大致为两个燃烧器之间的水冷壁和炉膛。
根据旋流燃烧器出口附近的流场特性，采用二维轴对称结构模拟该区域。
在计算区域的出口采用了倾斜一定角度的斜面以避免由于回流产生的压力计算不准确。
由于要计算旋转流动，为了得到较好的收敛结果，对燃烧器喉部壁面附近、水冷壁附近进行了网格细分。
利⽤fluent组分输运模型模拟锅炉混煤燃烧过程fluent组分输运模拟混煤燃烧之前⽤组分输运做过⼀些混煤燃烧的⼯作，因为⾃⼰⼀开始接触组分输运的时候也遇到很多困难，⽤组分输运做混煤模拟更是⼀⽆所知，后来在之前课题组基础上，加上⾃⼰的摸索，对⽤组分输运做混煤模拟的套路⼤概了解了，所以就把这个“套路”总结了⼀下写了出来，希望可以帮到有需要的朋友。
当然，下⾯的内容更多的是做混煤模拟的⼀个过程的描述，具体⾥⾯的有些参数的设置我也不是太懂，尤其是⼀些涉及到化学反应的参数，所以这篇⽂章只是告诉⼤家设置的“套路”，具体的参数还是要⼤家查阅相关⽂献或书籍。
另外⼤家也没必要死搬硬套我这个套路，我这篇⽂章只是希望能给想做混煤模拟的朋友⼀些启发，⼤家应当在我这个⽂章的基础上多去琢磨，搞清楚每⼀步的设置都是在做什么，这样⾃⼰遇到⼀些我⽂中没有提到的问题时也能⾃⼰解决。
⽂中若有什么错误或未描述清楚的地⽅，欢迎互相交流。
1.打开species⾯板，选择species transport（组分输运），Reactions勾选上Volumetric，表⽰组分输运在某体积内有化学反应，Turbulence-Chemistry Interaction点选Finite-Rate/Eddy-Dissipation，表⽰化学反应是有限反应速率的，反应速率受化学反应本⾝与湍流混合⼆者共同控制。
2.点Finite-Rate/Eddy-Dissipation后，下⾯会出现coal calculator，⽤于对煤的反应进⾏计算，点coal calculator，弹出如下界⾯，根据煤质分析结果，填⼊相应数据，这⾥假设有两种煤，⼀种中等挥发份，取名为coal-mv，⼀种为⾼挥发份，取名为coal-hv，相应结果如下图。
3.点完Apply后点OK，会弹出如下界⾯。
4.这样coal-hv就设置好了，然后继续点coal-calculator，以同样的⽅法设置coal-mv。
第六章，FLUENT中的燃烧模拟6.1 燃烧模拟的重要性●面向实际装置（如锅炉、内燃机、火箭发动机、火灾等）●面向实际现象（如点火、熄火、燃烧污染物生成等）6.2 FLUENT燃烧模拟方法概要●FLUENT可以模拟宽广范围内的燃烧（反应流）问题。
然而，需要注意的是：你必须保证你所使用的物理模型要适合你所研究的问题。
FLUENT在燃烧模拟中的应用可如下图所示：●气相燃烧模型一般的有限速率形式（Magnussen 模型）守恒标量的PDF模型（单或二组分混合物分数）层流火焰面模型（Laminar flamelet model）Zimont 模型●离散相模型煤燃烧与喷雾燃烧●热辐射模型DTRM, P-1, Rosseland 和Discrete Ordinates模型●污染物模型NO x 模型，烟（Soot）模型6.3 气相燃烧模型6.3.1 燃烧的化学动力学模拟实际中的燃烧过程是湍流和化学反应相互作用的结果，燃烧的化学反应速率是强非线性和强刚性的。
通常的化学反应机理包含了几十种组分和几百个基元反应，而且这些组分之间的反应时间尺度相差很大（10－9～102秒），因此在实际问题的求解过程中计算量和存储量极大，目前应用尚不现实。
在FLUENT 中，针对不同的燃烧现象，采用了不同的化学动力学处理手段，以减少计算成本，如下：● 有限速率燃烧模型——>预混、部分预混和扩散燃烧● 混合物分数方法（平衡化学的PDF 模型和非平衡化学的层流火焰面模型）——>扩散燃烧● 反应进度方法（Zimont 模型）——>预混燃烧● 混合物分数和反应进度方法的结合——>部分预混燃烧6.3.2一般的有限速率模型● 化学反应过程一般采用总包机理（即简化化学反应，如单步反应）进行描述 ● 求解组分的输运方程，得到每种组分的时均质量分数值，如下：6－1其中组分j 的反应源项为所有K 个反应中，组分j 的净生成速率：6－2 式中，反应k 中的组分j 的反应速率可按照Arrhenius 公式、混合（mixing ）速率或 “eddy breakup” 速率的方法求解。
计算流体力学作业FLUENT 模拟燃烧问题描述：长为2m、直径为0.45m的圆筒形燃烧器结构如图1所示，燃烧筒壁上嵌有三块厚为0.0005 m，高0.05 m的薄板，以利于甲烷与空气的混合。
燃烧火焰为湍流扩散火焰。
在燃烧器中心有一个直径为0.01 m、长为0.01 m、壁厚为0.002 m的小喷嘴，甲烷以60 m/s的速度从小喷嘴注入燃烧器。
空气从喷嘴周围以0.5 m/s的速度进入燃烧器。
总当量比大约是0.76(甲烷含量超过空气约28%)，甲烷气体在燃烧器中高速流动，并与低速流动的空气混合，基于甲烷喷嘴直径的雷诺数约为5.7×103。
假定燃料完全燃烧并转换为:CH4+2O2→CO2+2H2O反应过程是通过化学计量系数、形成焓和控制化学反应率的相应参数来定义的。
利用FLUENT的finite-rate化学反应模型对一个圆筒形燃烧器内的甲烷和空气的混合物的流动和燃烧过程进行研究。
1、建立物理模型，选择材料属性，定义带化学组分混合与反应的湍流流动边界条件2、使用非耦合求解器求解燃烧问题3、对燃烧组分的比热分别为常量和变量的情况进行计算，并比较其结果4、利用分布云图检查反应流的计算结果5、预测热力型和快速型的NO X含量6、使用场函数计算器进行NO含量计算一、利用GAMBIT建立计算模型第1步启动GAMBIT，建立基本结构分析：圆筒燃烧器是一个轴对称的结构，可简化为二维流动，故只要建立轴对称面上的二维结构就可以了，几何结构如图2所示。
(1)建立新文件夹在F盘根目录下建立一个名为combustion的文件夹。
(2)启动GAMBIT(3)创建对称轴①创建两端点。
A(0,0,0)，B(2,0,0)②将两端点连成线(4)创建小喷嘴及空气进口边界①创建C、D、E、F、G点C D E F Gx 0 0.01 0.01 0 0y 0.005 0.005 0.007 0.007 0.225②连接AC、CD、DE、DF、FG。
计算流体力学作业FLUENT 模拟燃烧问题描述：长为2m、直径为0.45m的圆筒形燃烧器结构如图1所示，燃烧筒壁上嵌有三块厚为0.0005 m，高0.05 m的薄板，以利于甲烷与空气的混合。
燃烧火焰为湍流扩散火焰。
在燃烧器中心有一个直径为0.01 m、长为0.01 m、壁厚为0.002 m的小喷嘴，甲烷以60 m/s的速度从小喷嘴注入燃烧器。
空气从喷嘴周围以0.5 m/s的速度进入燃烧器。
总当量比大约是0.76(甲烷含量超过空气约28%)，甲烷气体在燃烧器中高速流动，并与低速流动的空气混合，基于甲烷喷嘴直径的雷诺数约为5.7×103。
假定燃料完全燃烧并转换为:CH4+2O2→CO2+2H2O反应过程是通过化学计量系数、形成焓和控制化学反应率的相应参数来定义的。
利用FLUENT的finite-rate化学反应模型对一个圆筒形燃烧器内的甲烷和空气的混合物的流动和燃烧过程进行研究。
1、建立物理模型，选择材料属性，定义带化学组分混合与反应的湍流流动边界条件2、使用非耦合求解器求解燃烧问题3、对燃烧组分的比热分别为常量和变量的情况进行计算，并比较其结果4、利用分布云图检查反应流的计算结果5、预测热力型和快速型的NO X含量6、使用场函数计算器进行NO含量计算一、利用GAMBIT建立计算模型第1步启动GAMBIT，建立基本结构分析：圆筒燃烧器是一个轴对称的结构，可简化为二维流动，故只要建立轴对称面上的二维结构就可以了，几何结构如图2所示。
(1)建立新文件夹在F盘根目录下建立一个名为combustion的文件夹。
(2)启动GAMBIT(3)创建对称轴①创建两端点。
A(0,0,0)，B(2,0,0)②将两端点连成线(4)创建小喷嘴及空气进口边界①创建C、D、E、F、G点C D E F Gx 0 0.01 0.01 0 0y 0.005 0.005 0.007 0.007 0.225②连接AC、CD、DE、DF、FG。
Fluent-UDF_Coal-Combustion-Multiphase-Flow-Processes ---Codes for various kinds of heterogeneous reactions【第二部分-Part2】DEFINE_HET_RXN_RATE(SteamGasif_Rev_MGAS,c,t,hr,mw,yi,rr,rr_t){Thread **pt = THREAD_SUB_THREADS(t);Thread *tp = pt[0]; /* gas phase */int index_phase = Get_Phase_Index(hr);Thread *ts = pt[index_phase]; /* solid phase */*rr = 0;double direction = 0.0, mol_weight, y_carbon, rr_turb = 1e+20;/*CO + H2 ---> H2O + 1/25 SootThe reverse steam gasification reaction, which is CO + H2 ---> 1/25 Soot + H2O is written as the reactionshown above. So a negative rate means CO and H2 is consumed and H2O and Soot are generated. Note that noC(s) is generated, the stoich coeff for C(s) in the above reation is zero.Set the phase and species indices. Ash species index is initialized to zero, with all other indices.Ash species index is used as a flag to execute SetSpeciesIndex only once. This is done by the firstreaction, defined in the heterogeneous reaction panel in FLUENT GUI.*/if(IS_ASH == 0)SetSpeciesIndex();if(MGAS_Gasif){double RoRT = C_R(c,tp) * UNIVERSAL_GAS_CONSTANT * C_T(c,tp);double p_h2o = RoRT * yi[IP_H2O][IS_H2O]/mw[IP_H2O][IS_H2O]/ 101325.;double p_co = RoRT * yi[IP_CO][IS_CO]/mw[IP_CO][IS_CO] / 101325.;double p_h2 = RoRT * yi[IP_H2][IS_H2]/mw[IP_H2][IS_H2] / 101325.;y_carbon = yi[IP_SOOT][IS_SOOT];mol_weight = mw[IP_SOOT][IS_SOOT];if (rp_ke)rr_turb = Turbulent_rr(c, t, hr, yi);if(C_VOF(c, ts) >= eps_s_small){*rr = rr_steam_gasif(c, t, ts, tp, p_h2o, p_co, p_h2, y_carbon, mol_weight, &direction); /* mol/(cm^3 .s) */if( direction > 0.0) /* positive value implies C(s) + H2O ---> CO + H2 */*rr = 0.0;else /* negative value implies CO + H2 ---> H2O + 1/25 Soot */{*rr = abs(*rr);*rr = MIN(*rr, rr_turb);}}}}double rr_steam_gasif(cell_t c, Thread *t, Thread *ts, Thread *tp, double p_h2o, double p_co, double p_h2, double y_carbon, double mol_weight, double* direction){double rate, prod, T_g = MIN((MAX(TMIN,C_T(c,tp))),TMAX);double p_h2o_star = p_h2 * p_co / ( exp(17.29 - 16326/T_g) );if(MGAS_Gasif) *direction = p_h2o - p_h2o_star;if(PCCL_Gasif) *direction = pow(p_h2o, N_steam_gasification)/(1.+K_steam_gasification*p_h2);prod = y_carbon*(C_R(c,ts)*1e-03)/mol_weight*C_VOF(c,ts); /*1e-3 is to convert density from kg/m^3 to g/cm^3 */if(MGAS_Gasif && *direction < 0.0) /* this implies reverse H2O gasification */prod = y_carbon*(C_R(c,tp)*1e-03)/mol_weight*C_VOF(c,tp); /*1e-3 is to convert density from kg/m^3 to g/cm^3 */rate = A_steam_gasification*exp(-E_steam_gasification/Rgas/T_g)* Annealing_steam_gasification * prod * *direction; /* mol/cm^3.s */rate *= 1000.; /* kmol/(m^3 .s) */return rate;}DEFINE_HET_RXN_RATE(Co2Gasif,c,t,hr,mw,yi,rr,rr_t){Thread **pt = THREAD_SUB_THREADS(t);Thread *tp = pt[0]; /* gas phase */int index_phase = Get_Phase_Index(hr);Thread *ts = pt[index_phase]; /* solid phase */*rr = 0;double direction = 0.0, mol_weight, y_carbon;/*C(s) + CO2 ---> 2COSet the phase and species indices. Ash species index is initialized to zero, with all other indices.Ash species index is used as a flag to execute SetSpeciesIndex only once. This is done by the firstreaction, defined in the heterogeneous reaction panel in FLUENT GUI.*/if(IS_ASH == 0)SetSpeciesIndex();double RoRT = C_R(c,tp) * UNIVERSAL_GAS_CONSTANT * C_T(c,tp);double p_co = RoRT * yi[IP_CO][IS_CO]/mw[IP_CO][IS_CO] / 101325.;double p_co2 = RoRT * yi[IP_CO2][IS_CO2]/mw[IP_CO2][IS_CO2] / 101325.;SolidFuel_Reactant(c, t, hr, &y_carbon, &mol_weight);if(C_VOF(c, ts) >= eps_s_small){*rr = rr_co2_gasif(c, t, ts, tp, p_co, p_co2, y_carbon, mol_weight, &direction); /* mol/(cm^3 .s) */if( direction < 0.0) /* negative implies reverse steam gasification */*rr = 0.0;}}DEFINE_HET_RXN_RATE(Soot_CO2_Gasif,c,t,r,mw,yi,rr,rr_t){Thread **pt = THREAD_SUB_THREADS(t);Thread *tp = pt[0]; /* gas phase */*rr = 0;double rr_turb = 1e+20;double T_g = MIN((MAX(TMIN,C_T(c,tp))),TMAX);/*1/25 Soot + CO2 ---> 2COSet the phase and species indices. Ash species index is initialized to zero, with all other indices.Ash species index is used as a flag to execute SetSpeciesIndex only once. This is done by the firstreaction, defined in the heterogeneous reaction panel in FLUENT GUI.*/if(IS_ASH == 0)SetSpeciesIndex();double RoRT = C_R(c,tp) * UNIVERSAL_GAS_CONSTANT * C_T(c,tp);double p_co = RoRT * yi[IP_CO][IS_CO]/mw[IP_CO][IS_CO]/ 101325.;double p_co2 = RoRT * yi[IP_CO2][IS_CO2]/mw[IP_CO2][IS_CO2] / 101325.;if (rp_ke)rr_turb = Turbulent_rr(c, t, r, yi);double prod = yi[IP_SOOT][IS_SOOT]*(C_R(c,tp)*1e-03)/mw[IP_SOOT][IS_SOOT]*C_VOF(c,tp); /*1e-3 is to convert density from kg/m^3 to g/cm^3 */*rr = A_soot_co2_gasification*exp(-E_soot_co2_gasification/Rgas/T_g)* Annealing_soot_co2_gasification * prod *pow(p_co2, N_soot_co2_gasification)/(1.+K_soot_co2_gasification*p_co); /* mol/cm^3.s */*rr *= 1000.; /* kmol/(m^3 .s) */*rr = MIN(*rr, rr_turb);}DEFINE_HET_RXN_RATE(Co2Gasif_Rev_MGAS,c,t,hr,mw,yi,rr,rr_t){Thread **pt = THREAD_SUB_THREADS(t);Thread *tp = pt[0]; /* gas phase */int index_phase = Get_Phase_Index(hr);Thread *ts = pt[index_phase]; /* solid phase */*rr = 0;double direction = 0.0, mol_weight, y_carbon, rr_turb = 1e+20;/*2CO ---> CO2 + 1/25 SootThe reverse CO2 gasification reaction, which is 2CO ---> 1/25 Soot + CO2 is written as the reactionshown above. So a negative rate means CO is consumed and CO2 and Soot are generated. Note that noC(s) is generated, the stoich coeff for C(s) in the above reation is zero.Set the phase and species indices. Ash species index is initialized to zero, with all other indices.Ash species index is used as a flag to execute SetSpeciesIndex only once. This is done by the firstreaction, defined in the heterogeneous reaction panel in FLUENT GUI.*/if(IS_ASH == 0)SetSpeciesIndex();if(MGAS_Gasif){double RoRT = C_R(c,tp) * UNIVERSAL_GAS_CONSTANT * C_T(c,tp);double p_co = RoRT * yi[IP_CO][IS_CO]/mw[IP_CO][IS_CO] / 101325.;double p_co2 = RoRT * yi[IP_CO2][IS_CO2]/mw[IP_CO2][IS_CO2] / 101325.;y_carbon = yi[IP_SOOT][IS_SOOT];mol_weight = mw[IP_SOOT][IS_SOOT];if (rp_ke)rr_turb = Turbulent_rr(c, t, hr, yi);if(C_VOF(c, ts) >= eps_s_small){*rr = rr_co2_gasif(c, t, ts, tp, p_co, p_co2, y_carbon, mol_weight, &direction); /* mol/(cm^3 .s) */if( direction > 0.0) /* positive value implies C(s) + CO2 ---> 2CO */*rr = 0.0;else /* negative value implies 2CO ---> CO2 + 1/25 Soot */{*rr = abs(*rr);*rr = MIN(*rr, rr_turb);}}}}double rr_co2_gasif(cell_t c, Thread *t, Thread *ts, Thread *tp, double p_co, double p_co2, double y_carbon, double mol_weight, double* direction){double T_g = MIN(MAX(TMIN,C_T(c,tp)), TMAX), prod;double p_co2_star = p_co * p_co/(exp(20.92 - 20282/T_g));if(MGAS_Gasif) *direction = p_co2-p_co2_star;if(PCCL_Gasif) *direction = pow(p_co2, N_co2_gasification)/(1. + K_co2_gasification * p_co);prod = y_carbon*C_R(c,ts)*1.e-3/mol_weight* C_VOF(c,ts); /*1e-3 is to convert density from kg/m^3 to g/cm^3 */if(MGAS_Gasif && *direction < 0.0) /* this implies reverse CO2 gasification */prod = y_carbon*(C_R(c,tp)*1e-03)/mol_weight*C_VOF(c,tp); /*1e-3 is to convert density from kg/m^3 to g/cm^3 */double rate = A_co2_gasification*exp(-E_co2_gasification/Rgas/T_g)*Annealing_co2_gasification * prod * (*direction); /* mol/cm^3.s */rate *= 1000.; /* kmol/(m^3 .s) */return rate;}DEFINE_HET_RXN_RATE(H2Gasif,c,t,hr,mw,yi,rr,rr_t){Thread **pt = THREAD_SUB_THREADS(t);Thread *tp = pt[0]; /* gas phase */int index_phase = Get_Phase_Index(hr);Thread *ts = pt[index_phase]; /* solid phase */*rr = 0;double direction = 0.0, mol_weight, y_carbon;/*1/2 C(s) + H2 ---> 1/2 CH4Set the phase and species indices. Ash species index is initialized to zero, with all other indices.Ash species index is used as a flag to execute SetSpeciesIndex only once. This is done bythe firstreaction, defined in the heterogeneous reaction panel in FLUENT GUI.*/if(IS_ASH == 0)SetSpeciesIndex();double RoRT = C_R(c,tp) * UNIVERSAL_GAS_CONSTANT * C_T(c,tp);double p_h2 = RoRT * yi[IP_H2][IS_H2]/mw[IP_H2][IS_H2] / 101325.;double p_ch4 = RoRT * yi[IP_CH4][IS_CH4]/mw[IP_CH4][IS_CH4] / 101325.;SolidFuel_Reactant(c, t, hr, &y_carbon, &mol_weight);if(C_VOF(c, ts) >= eps_s_small){*rr = rr_h2_gasif(c, t, ts, tp, p_h2, p_ch4, y_carbon, mol_weight, &direction); /* mol/(cm^3 .s) */if( direction < 0.0) /* negative implies reverse steam gasification */*rr = 0.0;}}DEFINE_HET_RXN_RATE(Soot_H2_Gasif,c,t,r,mw,yi,rr,rr_t){Thread **pt = THREAD_SUB_THREADS(t);Thread *tp = pt[0]; /* gas phase */*rr = 0;double rr_turb = 1e+20;double T_g = MIN((MAX(TMIN,C_T(c,tp))),TMAX);/*1/25 Soot + 2H2 ---> CH4Set the phase and species indices. Ash species index is initialized to zero, with all other indices.Ash species index is used as a flag to execute SetSpeciesIndex only once. This is done by the firstreaction, defined in the heterogeneous reaction panel in FLUENT GUI.*/if(IS_ASH == 0)SetSpeciesIndex();double RoRT = C_R(c,tp) * UNIVERSAL_GAS_CONSTANT * C_T(c,tp);double p_h2 = RoRT * yi[IP_H2][IS_H2]/mw[IP_H2][IS_H2] / 101325.;if (rp_ke)rr_turb = Turbulent_rr(c, t, r, yi);double prod = yi[IP_SOOT][IS_SOOT]*(C_R(c,tp)*1e-03)/mw[IP_SOOT][IS_SOOT]*C_VOF(c,tp); /*1e-3 is to convert density from kg/m^3 to g/cm^3 */*rr = A_soot_h2_gasification*exp(-E_soot_h2_gasification/Rgas/T_g)* Annealing_soot_h2_gasification * prod *pow(p_h2, N_soot_h2_gasification); /* mol/cm^3.s */*rr *= 1000.; /* kmol/(m^3 .s) */*rr = MIN(*rr, rr_turb);}DEFINE_HET_RXN_RATE(H2Gasif_Rev_MGAS,c,t,hr,mw,yi,rr,rr_t){Thread **pt = THREAD_SUB_THREADS(t);Thread *tp = pt[0]; /* gas phase */int index_phase = Get_Phase_Index(hr);Thread *ts = pt[index_phase]; /* solid phase */*rr = 0;double direction = 0.0, mol_weight, y_carbon, rr_turb = 1e+20;/*1/2 CH4 ---> H2 + (0.5)*1/25 SootThe reverse H2 gasification reaction, which is 1/2 CH4 ---> 1/25 Soot + H2 is written as the reactionshown above. So a negative rate means CH4 is consumed and H2 and Soot are generated. Note that noC(s) is generated, the stoich coeff for C(s) in the above reation is zero.Set the phase and species indices. Ash species index is initialized to zero, with all other indices.Ash species index is used as a flag to execute SetSpeciesIndex only once. This is done by the firstreaction, defined in the heterogeneous reaction panel in FLUENT GUI.*/if(IS_ASH == 0)SetSpeciesIndex();if(MGAS_Gasif){double RoRT = C_R(c,tp) * UNIVERSAL_GAS_CONSTANT * C_T(c,tp);double p_h2 = RoRT * yi[IP_H2][IS_H2]/mw[IP_H2][IS_H2] / 101325.;double p_ch4 = RoRT * yi[IP_CH4][IS_CH4]/mw[IP_CH4][IS_CH4] / 101325.;y_carbon = yi[IP_SOOT][IS_SOOT];mol_weight = mw[IP_SOOT][IS_SOOT];if(C_VOF(c, ts) >= eps_s_small){if (rp_ke)rr_turb = Turbulent_rr(c, t, hr, yi);*rr = rr_h2_gasif(c, t, ts, tp, p_h2, p_ch4, y_carbon, mol_weight, &direction); /* mol/(cm^3 .s) */if( direction > 0.0) /* positive value implies 1/2 C(s) + H2 ---> 1/2 CH4 */*rr = 0.0;else /* negative value implies 1/2 CH4 ---> H2 + (0.5)*1/25 Soot */{*rr = abs(*rr);*rr = MIN(*rr, rr_turb);}}}}double rr_h2_gasif(cell_t c, Thread *t, Thread *ts, Thread *tp, double p_h2, double p_ch4, double y_carbon, double mol_weight, double* direction){double rate = 0.0, prod;double T_g = MIN((MAX(TMIN,C_T(c,tp))), TMAX);double p_h2_star = pow ((p_ch4/(exp(-13.43 + 10999/T_g))), 0.5);prod = y_carbon*C_R(c,ts)*1.e-3/mol_weight * C_VOF(c,ts); /*1e-3 is to convert density from kg/m^3 to g/cm^3 */if(MGAS_Gasif){*direction = p_h2-p_h2_star;if(*direction < 0.0) /* this implies reverse H2 gasification */prod = y_carbon*(C_R(c,tp)*1e-03)/mol_weight*C_VOF(c,tp); /*1e-3 is to convert density from kg/m^3 to g/cm^3 */rate = exp( -7.087 - 8078/T_g )* prod * *direction ; /* mol/cm^3.s */}if(PCCL_Gasif){*direction = p_h2;rate = A_h2_gasification*exp(-E_h2_gasification/Rgas/T_g)*Annealing_h2_gasification * prod * *direction; /* mol/cm^3.s */}rate *= 1000.; /* kmol/(m^3 .s) */return rate;}DEFINE_HET_RXN_RATE(coal_combustion,c,t,hr,mw,yi,rr,rr_t){Thread **pt = THREAD_SUB_THREADS(t);Thread *tp = pt[0]; /* gas phase */int index_phase = Get_Phase_Index(hr);Thread *ts = pt[index_phase]; /* solid phase */double mol_weight, y_carbon, y_ash;*rr = 0.0;/* Set the phase and species indices. Ash species index is initialized to zero, with all other indices.Ash species index is used as a flag to execute SetSpeciesIndex only once. This is done by the firstreaction, defined in the heterogeneous reaction panel in FLUENT GUI.*/if(IS_ASH == 0)SetSpeciesIndex();if( C_YI(c,tp,IS_O2) >= spe_small){SolidFuel_Reactant(c, t, hr, &y_carbon, &mol_weight);y_ash = yi[index_phase][IS_ASH];*rr = rr_combustion(c, t, ts, tp, yi[IP_O2][IS_O2], y_ash, y_carbon); /* mol/(cm^3 .s) */*rr *= 1000.; /* kmol/(m^3 .s) */}}double rr_combustion(cell_t c, Thread *t, Thread *ts, Thread *tp, double yi_O2, double y_ash,double y_carbon){double rd, k_f, k_r, factor, k_a, rate = 0.0, vrel;double Pt = MAX(0.1, (op_pres+C_P(c,t))/101325);double gas_constant = 82.06; /* atm.cm^3/mol.K */double T = C_T(c,tp), T_s = C_T(c,ts), D_p = C_PHASE_DIAMETER(c,ts)*100.;double p_o2 = C_R(c,tp)*UNIVERSAL_GAS_CONSTANT* T *yi_O2/mw[IP_O2][IS_O2] / 101325.; /* atm */if(fc_ar > 0.){if (y_ash > 0.){rd = pow( (y_carbon * ash_ar/100.)/(y_ash * fc_ar/100.), (1./3.) );rd = MIN(1., rd);}elserd = 1.;}elserd = 0.;double diff = MAX((4.26 * pow((T/1800.),1.75)/Pt), 1.e-10); /* cm^2/s */double Sc1o3 = pow(C_MU_L(c,tp)/(C_R(c,tp) * diff * 1.e-4), 1./3.);#if RP_2Dvrel = pow(( (C_U(c,tp)-C_U(c,ts))*(C_U(c,tp)-C_U(c,ts)) +(C_V(c,tp)-C_V(c,ts))*(C_V(c,tp)-C_V(c,ts))), 0.5);#endif#if RP_3Dvrel = pow(( (C_U(c,tp)-C_U(c,ts))*(C_U(c,tp)-C_U(c,ts)) +(C_V(c,tp)-C_V(c,ts))*(C_V(c,tp)-C_V(c,ts)) +(C_W(c,tp)-C_W(c,ts))*(C_W(c,tp)-C_W(c,ts)) ), 0.5);#endifdouble Re = C_VOF(c,tp) * D_p/100. * vrel * C_R(c,tp)/(C_MU_L(c,tp)+SMALL_S);double N_sherwood = (7. - 10. * C_VOF(c,tp) + 5. * C_VOF(c,tp) * C_VOF(c,tp) )*(1. + 0.7 * pow(Re, 0.2) * Sc1o3) +(1.33 - 2.4 * C_VOF(c,tp) + 1.2 * C_VOF(c,tp) * C_VOF(c,tp)) *pow(Re, 0.7) * Sc1o3;if ( rd <= 0. || C_VOF(c, ts) <= 0. ){rate = 0.;}else{k_f = diff * N_sherwood / (D_p * gas_constant/mw[IP_O2][IS_O2] * T ); /* g/(atm.cm^2.s) */k_r = A_c_combustion * exp( -E_c_combustion/Rgas/T_s ) * rd * rd;if ( rd >= 1.){rate = 1. / (1./k_f + 1./k_r);}else{k_a = 2. * rd * diff * f_ep_a / (D_p * (1.-rd) * gas_constant/mw[IP_O2][IS_O2] * T_s );rate = 1. / (1./k_f + 1./k_r + 1./k_a);}factor = y_carbon / (y_carbon + 1.e-6);rate *= p_o2 * 6. * C_VOF(c,ts) * factor / (D_p * 32.); /* mol/(cm^3 .s) */ }return rate;}#if !RP_NODE || !PARALLELvoid volatile_mass_fractions(){read_c3m_data();/* pan2 : Oct 2012 ... added CX_Messages for debugging */CX_Message("PCCL_Devol = %d\n",PCCL_Devol);CX_Message("MGAS_Devol = %d\n",MGAS_Devol);CX_Message("CPD_Devol = %d\n",CPD_Devol);CX_Message("FGDVC_Devol = %d\n",FGDVC_Devol);CX_Message("HPTR_Devol = %d\n",HPTR_Devol);CX_Message("MGAS_Moisture = %d\n",MGAS_Moisture);CX_Message("PCCL_Moisture = %d\n",PCCL_Moisture);CX_Message("MGAS_TarCracking = %d\n",MGAS_TarCracking);CX_Message("PCCL_2nd_Pyro = %d\n",PCCL_2nd_Pyro);CX_Message("MGAS_gasif = %d\n",MGAS_Gasif);CX_Message("PCCL_gasif = %d\n",PCCL_Gasif);CX_Message("PCCL_TarCracking = %d\n",PCCL_TarCracking);CX_Message("MGAS_WGS = %d\n",MGAS_WGS);CX_Message("PCCL_soot_gasif = %d\n",PCCL_soot_gasif);CX_Message("MGAS_char_combustion = %d\n",MGAS_char_combustion);CX_Message("PCCL_char_combustion = %d\n",PCCL_char_combustion);CX_Message("PCCL_soot_oxidation = %d\n",PCCL_soot_oxidation);CX_Message("TAR_oxidation = %d\n",TAR_oxidation);CX_Message("MGAS_gas_phase_oxidation = %d\n",MGAS_gas_phase_oxidation);CX_Message("fc_ar = %f\n",fc_ar);CX_Message("vm_ar = %f\n",vm_ar);CX_Message("ash_ar = %f\n",ash_ar);CX_Message("moist_ar = %f\n",moist_ar);CX_Message("a1_devolatilization = %f\n",A1_devolatilization);CX_Message("e1_devolatilization = %f\n",E1_devolatilization);CX_Message("a2_devolatilization = %f\n",A2_devolatilization);CX_Message("e2_devolatilization = %f\n",E2_devolatilization);CX_Message("a_tar_cracking = %f\n",A_tar_cracking);CX_Message("e_tar_cracking = %f\n",E_tar_cracking);CX_Message("A_steam_gasification = %f\n",A_steam_gasification);CX_Message("E_steam_gasification = %f\n",E_steam_gasification);CX_Message("K_steam_gasification = %f\n",K_steam_gasification);CX_Message("N_steam_gasification = %f\n",N_steam_gasification);CX_Message("Annealing_steam_gasification = %f\n",Annealing_steam_gasification);CX_Message("A_co2_gasification = %f\n",A_co2_gasification);CX_Message("E_co2_gasification = %f\n",E_co2_gasification);CX_Message("K_co2_gasification = %f\n",K_co2_gasification);CX_Message("N_co2_gasification = %f\n",N_co2_gasification);CX_Message("Annealing_co2_gasification = %f\n",Annealing_co2_gasification);CX_Message("A_h2_gasification = %f\n",A_h2_gasification);CX_Message("E_h2_gasification = %f\n",E_h2_gasification);CX_Message("N_h2_gasification = %f\n",N_h2_gasification);CX_Message("Annealing_h2_gasification = %f\n",Annealing_h2_gasification);CX_Message("A_soot_steam_gasification = %f\n",A_soot_steam_gasification);CX_Message("E_soot_steam_gasification = %f\n",E_soot_steam_gasification);CX_Message("K_soot_steam_gasification = %f\n",K_soot_steam_gasification);CX_Message("N_soot_steam_gasification = %f\n",N_soot_steam_gasification);CX_Message("Annealing_soot_steam_gasification= %f\n",Annealing_soot_steam_gasification);CX_Message("A_soot_co2_gasification = %f\n",A_soot_co2_gasification);CX_Message("E_soot_co2_gasification = %f\n",E_soot_co2_gasification);CX_Message("K_soot_co2_gasification = %f\n",K_soot_co2_gasification);CX_Message("N_soot_co2_gasification = %f\n",N_soot_co2_gasification);CX_Message("Annealing_soot_co2_gasification= %f\n",Annealing_soot_co2_gasification);CX_Message("A_soot_h2_gasification = %f\n",A_soot_h2_gasification);CX_Message("E_soot_h2_gasification = %f\n",E_soot_h2_gasification);CX_Message("N_soot_h2_gasification = %f\n",N_soot_h2_gasification);CX_Message("Annealing_soot_h2_gasification = %f\n",Annealing_soot_h2_gasification);CX_Message("A_soot_combustion = %f\n",A_Soot_Combustion);CX_Message("E_soot_combustion = %f\n",E_Soot_Combustion);CX_Message("A_c_combustion = %f\n",A_c_combustion);CX_Message("E_c_combustion = %f\n",E_c_combustion);CX_Message("N_c_combustion = %f\n",N_c_combustion);CX_Message("Annealing_c_combustion = %f\n",Annealing_c_combustion);CX_Message("A_moisture_release = %f\n",A_moisture_release);CX_Message("E_moisture_release = %f\n",E_moisture_release);CX_Message("wg3 = %f\n",wg3);CX_Message("moisture_flux = %f\n",Moisture_Flux);/* f_ep_a is used in shrinking core model, in coal combustion model */double ep_a = 0.25 + 0.75*(1-ash_ar/100.);f_ep_a = pow(ep_a,2.5);}/* pan c3m start */static void SetBooleanValue(char * var , char * svalue) /* pan : Oct 2012 ... new function */ {cxboolean value;if ( strcmp(svalue,"true") == 0) /* pan2 : Oct 2012 : correction */value = TRUE;elsevalue = FALSE;if (strcmp(var,"pccl_devol") == 0) PCCL_Devol = value;if (strcmp(var,"mgas_devol") == 0) MGAS_Devol = value;if (strcmp(var,"cpd_devol") == 0) CPD_Devol = value;if (strcmp(var,"fgdvc_devol") == 0) FGDVC_Devol = value;if (strcmp(var,"hptr_devol") == 0) HPTR_Devol = value;if (strcmp(var,"mgas_moisture") == 0) MGAS_Moisture = value;if (strcmp(var,"pccl_moisture") == 0) PCCL_Moisture = value;if (strcmp(var,"mgas_tarcracking") == 0) MGAS_TarCracking = value;if (strcmp(var,"pccl_2nd_pyro") == 0) PCCL_2nd_Pyro = value;if (strcmp(var,"mgas_gasif") == 0) MGAS_Gasif = value;if (strcmp(var,"pccl_gasif") == 0) PCCL_Gasif = value;if (strcmp(var,"pccl_tarcracking") == 0) PCCL_TarCracking = value;if (strcmp(var,"mgas_wgs") == 0) MGAS_WGS = value;if (strcmp(var,"pccl_soot_gasif") == 0) PCCL_soot_gasif = value;if (strcmp(var,"mgas_char_combustion") == 0) MGAS_char_combustion = value;if (strcmp(var,"pccl_char_combustion") == 0) PCCL_char_combustion = value;if (strcmp(var,"pccl_soot_oxidation") == 0) PCCL_soot_oxidation = value;if (strcmp(var,"tar_oxidation") == 0) TAR_oxidation = value;if (strcmp(var,"mgas_gas_phase_oxidation") == 0) MGAS_gas_phase_oxidation = value;}static void SetValue(char * var , char * svalue) /* pan : oct 2012 ... replace entire function */{char * pEnd;double value = strtod(svalue,&pEnd);if (strcmp(var,"fc_ar") == 0) fc_ar = value;if (strcmp(var,"vm_ar") == 0) vm_ar = value;if (strcmp(var,"ash_ar") == 0) ash_ar = value;if (strcmp(var,"moist_ar") == 0) moist_ar = value;/* NOTE:If the two stage devolatilization model is not used, A2 and E2 have to be set to zero*/if (strcmp(var,"a1_devolatilization") == 0) A1_devolatilization = value;if (strcmp(var,"e1_devolatilization") == 0) E1_devolatilization = value;if (strcmp(var,"a2_devolatilization") == 0) A2_devolatilization = value;if (strcmp(var,"e2_devolatilization") == 0) E2_devolatilization = value;/* if (strcmp(var,"c3m_aem") == 0) c3m_aem = value; should not be used pan : Oct 2012 */if (strcmp(var,"a_tar_cracking") == 0) A_tar_cracking = value;if (strcmp(var,"e_tar_cracking") == 0) E_tar_cracking = value;if (strcmp(var,"a_steam_gasification") == 0) A_steam_gasification = value;if (strcmp(var,"e_steam_gasification") == 0) E_steam_gasification = value;if (strcmp(var,"annealing_steam_gasification") == 0) Annealing_steam_gasification = value;if (strcmp(var,"k_steam_gasification") == 0) K_steam_gasification = value;if (strcmp(var,"n_steam_gasification") == 0) N_steam_gasification = value;if (strcmp(var,"a_co2_gasification") == 0) A_co2_gasification = value;if (strcmp(var,"e_co2_gasification") == 0) E_co2_gasification = value;if (strcmp(var,"annealing_co2_gasification") == 0) Annealing_co2_gasification = value;if (strcmp(var,"k_co2_gasification") == 0) K_co2_gasification = value;if (strcmp(var,"n_co2_gasification") == 0) N_co2_gasification = value;if (strcmp(var,"a_h2_gasification") == 0) A_h2_gasification = value;if (strcmp(var,"e_h2_gasification") == 0) E_h2_gasification = value;if (strcmp(var,"annealing_h2_gasification") == 0) Annealing_h2_gasification = value;if (strcmp(var,"n_h2_gasification") == 0) N_h2_gasification = value;if (strcmp(var,"a_soot_steam_gasification") == 0) A_soot_steam_gasification = value;if (strcmp(var,"e_soot_steam_gasification") == 0) E_soot_steam_gasification = value;if (strcmp(var,"annealing_soot_steam_gasification") == 0) Annealing_soot_steam_gasification = value;if (strcmp(var,"k_soot_steam_gasification") == 0) K_soot_steam_gasification = value;if (strcmp(var,"n_soot_steam_gasification") == 0) N_soot_steam_gasification = value;if (strcmp(var,"a_soot_co2_gasification") == 0) A_soot_co2_gasification = value;if (strcmp(var,"e_soot_co2_gasification") == 0) E_soot_co2_gasification = value;if (strcmp(var,"annealing_soot_co2_gasification") == 0) Annealing_soot_co2_gasification = value;if (strcmp(var,"k_soot_co2_gasification") == 0) K_soot_co2_gasification = value;if (strcmp(var,"n_soot_co2_gasification") == 0) N_soot_co2_gasification= value;if (strcmp(var,"a_soot_h2_gasification") == 0) A_soot_h2_gasification = value;if (strcmp(var,"e_soot_h2_gasification") == 0) E_soot_h2_gasification = value;if (strcmp(var,"annealing_soot_h2_gasification") == 0) Annealing_soot_h2_gasification = value;if (strcmp(var,"n_soot_h2_gasification") == 0) N_soot_h2_gasification = value;if (strcmp(var,"a_soot_combustion") == 0) A_Soot_Combustion = value;if (strcmp(var,"e_soot_combustion") == 0) E_Soot_Combustion = value;if (strcmp(var,"a_c_combustion") == 0) A_c_combustion = value;if (strcmp(var,"e_c_combustion") == 0) E_c_combustion = value;if (strcmp(var,"n_c_combustion") == 0) N_c_combustion = value;if (strcmp(var,"annealing_c_combustion") == 0) Annealing_c_combustion = value;if (strcmp(var,"a_moisture_release") == 0) A_moisture_release = value;if (strcmp(var,"e_moisture_release") == 0) E_moisture_release = value;if (strcmp(var,"c3m_wg3") == 0) wg3 = value;if (strcmp(var,"moisture_flux") == 0) Moisture_Flux = value;/* if (strcmp(var,"c3m_ae5") == 0) c3m_ae5 = value; pan : Oct 2012 ... should not be used */}void read_c3m_data(){CX_Message("start of read_c3m_data \n");FILE * pFile;char line[80];char field1[80];char field2[80];char *pch;int i , field_index;pFile = fopen("fluent_c3m_udf.inp","r");if (pFile != NULL){CX_Message("fopen OK \n");char * pEnd = line;while (pEnd != NULL){/* CX_Message("read a line - start \n"); */pEnd = fgets(line,80,pFile);if (pEnd != NULL){for (i=0; i<80; ++i) line[i] = tolower(line[i]);pch = strtok(line," ,\t\n=");field_index = 0;while (pch != NULL){if (field_index == 0){strcpy(field1,pch);field_index = 1;/* CX_Message("token %s \n",field1); */}else{strcpy(field2,pch);field_index = 0;SetValue(field1,field2);SetBooleanValue(field1,field2); /* pan : Oct 2012 *//*CX_Message("\n\n");CX_Message("token %s \n",field1);CX_Message("token %s \n",field2);*/}pch = strtok(NULL," ,\t\n=");}}}fclose(pFile);}。
Fluent-UDF_Coal-Combustion-Multiphase-Flow-Processes ---Codes for various kinds of heterogeneous reactions【第二部分-Part2】DEFINE_HET_RXN_RATE(SteamGasif_Rev_MGAS,c,t,hr,mw,yi,rr,rr_t){Thread **pt = THREAD_SUB_THREADS(t);Thread *tp = pt[0]; /* gas phase */int index_phase = Get_Phase_Index(hr);Thread *ts = pt[index_phase]; /* solid phase */*rr = 0;double direction = 0.0, mol_weight, y_carbon, rr_turb = 1e+20;/*CO + H2 ---> H2O + 1/25 SootThe reverse steam gasification reaction, which is CO + H2 ---> 1/25 Soot + H2O is written as the reactionshown above. So a negative rate means CO and H2 is consumed and H2O and Soot are generated. Note that noC(s) is generated, the stoich coeff for C(s) in the above reation is zero.Set the phase and species indices. Ash species index is initialized to zero, with all other indices.Ash species index is used as a flag to execute SetSpeciesIndex only once. This is done by the firstreaction, defined in the heterogeneous reaction panel in FLUENT GUI.*/if(IS_ASH == 0)SetSpeciesIndex();if(MGAS_Gasif){double RoRT = C_R(c,tp) * UNIVERSAL_GAS_CONSTANT * C_T(c,tp);double p_h2o = RoRT * yi[IP_H2O][IS_H2O]/mw[IP_H2O][IS_H2O]/ 101325.;double p_co = RoRT * yi[IP_CO][IS_CO]/mw[IP_CO][IS_CO] / 101325.;double p_h2 = RoRT * yi[IP_H2][IS_H2]/mw[IP_H2][IS_H2] / 101325.;y_carbon = yi[IP_SOOT][IS_SOOT];mol_weight = mw[IP_SOOT][IS_SOOT];if (rp_ke)rr_turb = Turbulent_rr(c, t, hr, yi);if(C_VOF(c, ts) >= eps_s_small){*rr = rr_steam_gasif(c, t, ts, tp, p_h2o, p_co, p_h2, y_carbon, mol_weight, &direction); /* mol/(cm^3 .s) */if( direction > 0.0) /* positive value implies C(s) + H2O ---> CO + H2 */*rr = 0.0;else /* negative value implies CO + H2 ---> H2O + 1/25 Soot */{*rr = abs(*rr);*rr = MIN(*rr, rr_turb);}}}}double rr_steam_gasif(cell_t c, Thread *t, Thread *ts, Thread *tp, double p_h2o, double p_co, double p_h2, double y_carbon, double mol_weight, double* direction){double rate, prod, T_g = MIN((MAX(TMIN,C_T(c,tp))),TMAX);double p_h2o_star = p_h2 * p_co / ( exp(17.29 - 16326/T_g) );if(MGAS_Gasif) *direction = p_h2o - p_h2o_star;if(PCCL_Gasif) *direction = pow(p_h2o, N_steam_gasification)/(1.+K_steam_gasification*p_h2);prod = y_carbon*(C_R(c,ts)*1e-03)/mol_weight*C_VOF(c,ts); /*1e-3 is to convert density from kg/m^3 to g/cm^3 */if(MGAS_Gasif && *direction < 0.0) /* this implies reverse H2O gasification */prod = y_carbon*(C_R(c,tp)*1e-03)/mol_weight*C_VOF(c,tp); /*1e-3 is to convert density from kg/m^3 to g/cm^3 */rate = A_steam_gasification*exp(-E_steam_gasification/Rgas/T_g)* Annealing_steam_gasification * prod * *direction; /* mol/cm^3.s */rate *= 1000.; /* kmol/(m^3 .s) */return rate;}DEFINE_HET_RXN_RATE(Co2Gasif,c,t,hr,mw,yi,rr,rr_t){Thread **pt = THREAD_SUB_THREADS(t);Thread *tp = pt[0]; /* gas phase */int index_phase = Get_Phase_Index(hr);Thread *ts = pt[index_phase]; /* solid phase */*rr = 0;double direction = 0.0, mol_weight, y_carbon;/*C(s) + CO2 ---> 2COSet the phase and species indices. Ash species index is initialized to zero, with all other indices.Ash species index is used as a flag to execute SetSpeciesIndex only once. This is done by the firstreaction, defined in the heterogeneous reaction panel in FLUENT GUI.*/if(IS_ASH == 0)SetSpeciesIndex();double RoRT = C_R(c,tp) * UNIVERSAL_GAS_CONSTANT * C_T(c,tp);double p_co = RoRT * yi[IP_CO][IS_CO]/mw[IP_CO][IS_CO] / 101325.;double p_co2 = RoRT * yi[IP_CO2][IS_CO2]/mw[IP_CO2][IS_CO2] / 101325.;SolidFuel_Reactant(c, t, hr, &y_carbon, &mol_weight);if(C_VOF(c, ts) >= eps_s_small){*rr = rr_co2_gasif(c, t, ts, tp, p_co, p_co2, y_carbon, mol_weight, &direction); /* mol/(cm^3 .s) */if( direction < 0.0) /* negative implies reverse steam gasification */*rr = 0.0;}}DEFINE_HET_RXN_RATE(Soot_CO2_Gasif,c,t,r,mw,yi,rr,rr_t){Thread **pt = THREAD_SUB_THREADS(t);Thread *tp = pt[0]; /* gas phase */*rr = 0;double rr_turb = 1e+20;double T_g = MIN((MAX(TMIN,C_T(c,tp))),TMAX);/*1/25 Soot + CO2 ---> 2COSet the phase and species indices. Ash species index is initialized to zero, with all other indices.Ash species index is used as a flag to execute SetSpeciesIndex only once. This is done by the firstreaction, defined in the heterogeneous reaction panel in FLUENT GUI.*/if(IS_ASH == 0)SetSpeciesIndex();double RoRT = C_R(c,tp) * UNIVERSAL_GAS_CONSTANT * C_T(c,tp);double p_co = RoRT * yi[IP_CO][IS_CO]/mw[IP_CO][IS_CO]/ 101325.;double p_co2 = RoRT * yi[IP_CO2][IS_CO2]/mw[IP_CO2][IS_CO2] / 101325.;if (rp_ke)rr_turb = Turbulent_rr(c, t, r, yi);double prod = yi[IP_SOOT][IS_SOOT]*(C_R(c,tp)*1e-03)/mw[IP_SOOT][IS_SOOT]*C_VOF(c,tp); /*1e-3 is to convert density from kg/m^3 to g/cm^3 */*rr = A_soot_co2_gasification*exp(-E_soot_co2_gasification/Rgas/T_g)* Annealing_soot_co2_gasification * prod *pow(p_co2, N_soot_co2_gasification)/(1.+K_soot_co2_gasification*p_co); /* mol/cm^3.s */*rr *= 1000.; /* kmol/(m^3 .s) */*rr = MIN(*rr, rr_turb);}DEFINE_HET_RXN_RATE(Co2Gasif_Rev_MGAS,c,t,hr,mw,yi,rr,rr_t){Thread **pt = THREAD_SUB_THREADS(t);Thread *tp = pt[0]; /* gas phase */int index_phase = Get_Phase_Index(hr);Thread *ts = pt[index_phase]; /* solid phase */*rr = 0;double direction = 0.0, mol_weight, y_carbon, rr_turb = 1e+20;/*2CO ---> CO2 + 1/25 SootThe reverse CO2 gasification reaction, which is 2CO ---> 1/25 Soot + CO2 is written as the reactionshown above. So a negative rate means CO is consumed and CO2 and Soot are generated. Note that noC(s) is generated, the stoich coeff for C(s) in the above reation is zero.Set the phase and species indices. Ash species index is initialized to zero, with all other indices.Ash species index is used as a flag to execute SetSpeciesIndex only once. This is done by the firstreaction, defined in the heterogeneous reaction panel in FLUENT GUI.*/if(IS_ASH == 0)SetSpeciesIndex();if(MGAS_Gasif){double RoRT = C_R(c,tp) * UNIVERSAL_GAS_CONSTANT * C_T(c,tp);double p_co = RoRT * yi[IP_CO][IS_CO]/mw[IP_CO][IS_CO] / 101325.;double p_co2 = RoRT * yi[IP_CO2][IS_CO2]/mw[IP_CO2][IS_CO2] / 101325.;y_carbon = yi[IP_SOOT][IS_SOOT];mol_weight = mw[IP_SOOT][IS_SOOT];if (rp_ke)rr_turb = Turbulent_rr(c, t, hr, yi);if(C_VOF(c, ts) >= eps_s_small){*rr = rr_co2_gasif(c, t, ts, tp, p_co, p_co2, y_carbon, mol_weight, &direction); /* mol/(cm^3 .s) */if( direction > 0.0) /* positive value implies C(s) + CO2 ---> 2CO */*rr = 0.0;else /* negative value implies 2CO ---> CO2 + 1/25 Soot */{*rr = abs(*rr);*rr = MIN(*rr, rr_turb);}}}}double rr_co2_gasif(cell_t c, Thread *t, Thread *ts, Thread *tp, double p_co, double p_co2, double y_carbon, double mol_weight, double* direction){double T_g = MIN(MAX(TMIN,C_T(c,tp)), TMAX), prod;double p_co2_star = p_co * p_co/(exp(20.92 - 20282/T_g));if(MGAS_Gasif) *direction = p_co2-p_co2_star;if(PCCL_Gasif) *direction = pow(p_co2, N_co2_gasification)/(1. + K_co2_gasification * p_co);prod = y_carbon*C_R(c,ts)*1.e-3/mol_weight* C_VOF(c,ts); /*1e-3 is to convert density from kg/m^3 to g/cm^3 */if(MGAS_Gasif && *direction < 0.0) /* this implies reverse CO2 gasification */prod = y_carbon*(C_R(c,tp)*1e-03)/mol_weight*C_VOF(c,tp); /*1e-3 is to convert density from kg/m^3 to g/cm^3 */double rate = A_co2_gasification*exp(-E_co2_gasification/Rgas/T_g)*Annealing_co2_gasification * prod * (*direction); /* mol/cm^3.s */rate *= 1000.; /* kmol/(m^3 .s) */return rate;}DEFINE_HET_RXN_RATE(H2Gasif,c,t,hr,mw,yi,rr,rr_t){Thread **pt = THREAD_SUB_THREADS(t);Thread *tp = pt[0]; /* gas phase */int index_phase = Get_Phase_Index(hr);Thread *ts = pt[index_phase]; /* solid phase */*rr = 0;double direction = 0.0, mol_weight, y_carbon;/*1/2 C(s) + H2 ---> 1/2 CH4Set the phase and species indices. Ash species index is initialized to zero, with all other indices.Ash species index is used as a flag to execute SetSpeciesIndex only once. This is done bythe firstreaction, defined in the heterogeneous reaction panel in FLUENT GUI.*/if(IS_ASH == 0)SetSpeciesIndex();double RoRT = C_R(c,tp) * UNIVERSAL_GAS_CONSTANT * C_T(c,tp);double p_h2 = RoRT * yi[IP_H2][IS_H2]/mw[IP_H2][IS_H2] / 101325.;double p_ch4 = RoRT * yi[IP_CH4][IS_CH4]/mw[IP_CH4][IS_CH4] / 101325.;SolidFuel_Reactant(c, t, hr, &y_carbon, &mol_weight);if(C_VOF(c, ts) >= eps_s_small){*rr = rr_h2_gasif(c, t, ts, tp, p_h2, p_ch4, y_carbon, mol_weight, &direction); /* mol/(cm^3 .s) */if( direction < 0.0) /* negative implies reverse steam gasification */*rr = 0.0;}}DEFINE_HET_RXN_RATE(Soot_H2_Gasif,c,t,r,mw,yi,rr,rr_t){Thread **pt = THREAD_SUB_THREADS(t);Thread *tp = pt[0]; /* gas phase */*rr = 0;double rr_turb = 1e+20;double T_g = MIN((MAX(TMIN,C_T(c,tp))),TMAX);/*1/25 Soot + 2H2 ---> CH4Set the phase and species indices. Ash species index is initialized to zero, with all other indices.Ash species index is used as a flag to execute SetSpeciesIndex only once. This is done by the firstreaction, defined in the heterogeneous reaction panel in FLUENT GUI.*/if(IS_ASH == 0)SetSpeciesIndex();double RoRT = C_R(c,tp) * UNIVERSAL_GAS_CONSTANT * C_T(c,tp);double p_h2 = RoRT * yi[IP_H2][IS_H2]/mw[IP_H2][IS_H2] / 101325.;if (rp_ke)rr_turb = Turbulent_rr(c, t, r, yi);double prod = yi[IP_SOOT][IS_SOOT]*(C_R(c,tp)*1e-03)/mw[IP_SOOT][IS_SOOT]*C_VOF(c,tp); /*1e-3 is to convert density from kg/m^3 to g/cm^3 */*rr = A_soot_h2_gasification*exp(-E_soot_h2_gasification/Rgas/T_g)* Annealing_soot_h2_gasification * prod *pow(p_h2, N_soot_h2_gasification); /* mol/cm^3.s */*rr *= 1000.; /* kmol/(m^3 .s) */*rr = MIN(*rr, rr_turb);}DEFINE_HET_RXN_RATE(H2Gasif_Rev_MGAS,c,t,hr,mw,yi,rr,rr_t){Thread **pt = THREAD_SUB_THREADS(t);Thread *tp = pt[0]; /* gas phase */int index_phase = Get_Phase_Index(hr);Thread *ts = pt[index_phase]; /* solid phase */*rr = 0;double direction = 0.0, mol_weight, y_carbon, rr_turb = 1e+20;/*1/2 CH4 ---> H2 + (0.5)*1/25 SootThe reverse H2 gasification reaction, which is 1/2 CH4 ---> 1/25 Soot + H2 is written as the reactionshown above. So a negative rate means CH4 is consumed and H2 and Soot are generated. Note that noC(s) is generated, the stoich coeff for C(s) in the above reation is zero.Set the phase and species indices. Ash species index is initialized to zero, with all other indices.Ash species index is used as a flag to execute SetSpeciesIndex only once. This is done by the firstreaction, defined in the heterogeneous reaction panel in FLUENT GUI.*/if(IS_ASH == 0)SetSpeciesIndex();if(MGAS_Gasif){double RoRT = C_R(c,tp) * UNIVERSAL_GAS_CONSTANT * C_T(c,tp);double p_h2 = RoRT * yi[IP_H2][IS_H2]/mw[IP_H2][IS_H2] / 101325.;double p_ch4 = RoRT * yi[IP_CH4][IS_CH4]/mw[IP_CH4][IS_CH4] / 101325.;y_carbon = yi[IP_SOOT][IS_SOOT];mol_weight = mw[IP_SOOT][IS_SOOT];if(C_VOF(c, ts) >= eps_s_small){if (rp_ke)rr_turb = Turbulent_rr(c, t, hr, yi);*rr = rr_h2_gasif(c, t, ts, tp, p_h2, p_ch4, y_carbon, mol_weight, &direction); /* mol/(cm^3 .s) */if( direction > 0.0) /* positive value implies 1/2 C(s) + H2 ---> 1/2 CH4 */*rr = 0.0;else /* negative value implies 1/2 CH4 ---> H2 + (0.5)*1/25 Soot */{*rr = abs(*rr);*rr = MIN(*rr, rr_turb);}}}}double rr_h2_gasif(cell_t c, Thread *t, Thread *ts, Thread *tp, double p_h2, double p_ch4, double y_carbon, double mol_weight, double* direction){double rate = 0.0, prod;double T_g = MIN((MAX(TMIN,C_T(c,tp))), TMAX);double p_h2_star = pow ((p_ch4/(exp(-13.43 + 10999/T_g))), 0.5);prod = y_carbon*C_R(c,ts)*1.e-3/mol_weight * C_VOF(c,ts); /*1e-3 is to convert density from kg/m^3 to g/cm^3 */if(MGAS_Gasif){*direction = p_h2-p_h2_star;if(*direction < 0.0) /* this implies reverse H2 gasification */prod = y_carbon*(C_R(c,tp)*1e-03)/mol_weight*C_VOF(c,tp); /*1e-3 is to convert density from kg/m^3 to g/cm^3 */rate = exp( -7.087 - 8078/T_g )* prod * *direction ; /* mol/cm^3.s */}if(PCCL_Gasif){*direction = p_h2;rate = A_h2_gasification*exp(-E_h2_gasification/Rgas/T_g)*Annealing_h2_gasification * prod * *direction; /* mol/cm^3.s */}rate *= 1000.; /* kmol/(m^3 .s) */return rate;}DEFINE_HET_RXN_RATE(coal_combustion,c,t,hr,mw,yi,rr,rr_t){Thread **pt = THREAD_SUB_THREADS(t);Thread *tp = pt[0]; /* gas phase */int index_phase = Get_Phase_Index(hr);Thread *ts = pt[index_phase]; /* solid phase */double mol_weight, y_carbon, y_ash;*rr = 0.0;/* Set the phase and species indices. Ash species index is initialized to zero, with all other indices.Ash species index is used as a flag to execute SetSpeciesIndex only once. This is done by the firstreaction, defined in the heterogeneous reaction panel in FLUENT GUI.*/if(IS_ASH == 0)SetSpeciesIndex();if( C_YI(c,tp,IS_O2) >= spe_small){SolidFuel_Reactant(c, t, hr, &y_carbon, &mol_weight);y_ash = yi[index_phase][IS_ASH];*rr = rr_combustion(c, t, ts, tp, yi[IP_O2][IS_O2], y_ash, y_carbon); /* mol/(cm^3 .s) */*rr *= 1000.; /* kmol/(m^3 .s) */}}double rr_combustion(cell_t c, Thread *t, Thread *ts, Thread *tp, double yi_O2, double y_ash,double y_carbon){double rd, k_f, k_r, factor, k_a, rate = 0.0, vrel;double Pt = MAX(0.1, (op_pres+C_P(c,t))/101325);double gas_constant = 82.06; /* atm.cm^3/mol.K */double T = C_T(c,tp), T_s = C_T(c,ts), D_p = C_PHASE_DIAMETER(c,ts)*100.;double p_o2 = C_R(c,tp)*UNIVERSAL_GAS_CONSTANT* T *yi_O2/mw[IP_O2][IS_O2] / 101325.; /* atm */if(fc_ar > 0.){if (y_ash > 0.){rd = pow( (y_carbon * ash_ar/100.)/(y_ash * fc_ar/100.), (1./3.) );rd = MIN(1., rd);}elserd = 1.;}elserd = 0.;double diff = MAX((4.26 * pow((T/1800.),1.75)/Pt), 1.e-10); /* cm^2/s */double Sc1o3 = pow(C_MU_L(c,tp)/(C_R(c,tp) * diff * 1.e-4), 1./3.);#if RP_2Dvrel = pow(( (C_U(c,tp)-C_U(c,ts))*(C_U(c,tp)-C_U(c,ts)) +(C_V(c,tp)-C_V(c,ts))*(C_V(c,tp)-C_V(c,ts))), 0.5);#endif#if RP_3Dvrel = pow(( (C_U(c,tp)-C_U(c,ts))*(C_U(c,tp)-C_U(c,ts)) +(C_V(c,tp)-C_V(c,ts))*(C_V(c,tp)-C_V(c,ts)) +(C_W(c,tp)-C_W(c,ts))*(C_W(c,tp)-C_W(c,ts)) ), 0.5);#endifdouble Re = C_VOF(c,tp) * D_p/100. * vrel * C_R(c,tp)/(C_MU_L(c,tp)+SMALL_S);double N_sherwood = (7. - 10. * C_VOF(c,tp) + 5. * C_VOF(c,tp) * C_VOF(c,tp) )*(1. + 0.7 * pow(Re, 0.2) * Sc1o3) +(1.33 - 2.4 * C_VOF(c,tp) + 1.2 * C_VOF(c,tp) * C_VOF(c,tp)) *pow(Re, 0.7) * Sc1o3;if ( rd <= 0. || C_VOF(c, ts) <= 0. ){rate = 0.;}else{k_f = diff * N_sherwood / (D_p * gas_constant/mw[IP_O2][IS_O2] * T ); /* g/(atm.cm^2.s) */k_r = A_c_combustion * exp( -E_c_combustion/Rgas/T_s ) * rd * rd;if ( rd >= 1.){rate = 1. / (1./k_f + 1./k_r);}else{k_a = 2. * rd * diff * f_ep_a / (D_p * (1.-rd) * gas_constant/mw[IP_O2][IS_O2] * T_s );rate = 1. / (1./k_f + 1./k_r + 1./k_a);}factor = y_carbon / (y_carbon + 1.e-6);rate *= p_o2 * 6. * C_VOF(c,ts) * factor / (D_p * 32.); /* mol/(cm^3 .s) */ }return rate;}#if !RP_NODE || !PARALLELvoid volatile_mass_fractions(){read_c3m_data();/* pan2 : Oct 2012 ... added CX_Messages for debugging */CX_Message("PCCL_Devol = %d\n",PCCL_Devol);CX_Message("MGAS_Devol = %d\n",MGAS_Devol);CX_Message("CPD_Devol = %d\n",CPD_Devol);CX_Message("FGDVC_Devol = %d\n",FGDVC_Devol);CX_Message("HPTR_Devol = %d\n",HPTR_Devol);CX_Message("MGAS_Moisture = %d\n",MGAS_Moisture);CX_Message("PCCL_Moisture = %d\n",PCCL_Moisture);CX_Message("MGAS_TarCracking = %d\n",MGAS_TarCracking);CX_Message("PCCL_2nd_Pyro = %d\n",PCCL_2nd_Pyro);CX_Message("MGAS_gasif = %d\n",MGAS_Gasif);CX_Message("PCCL_gasif = %d\n",PCCL_Gasif);CX_Message("PCCL_TarCracking = %d\n",PCCL_TarCracking);CX_Message("MGAS_WGS = %d\n",MGAS_WGS);CX_Message("PCCL_soot_gasif = %d\n",PCCL_soot_gasif);CX_Message("MGAS_char_combustion = %d\n",MGAS_char_combustion);CX_Message("PCCL_char_combustion = %d\n",PCCL_char_combustion);CX_Message("PCCL_soot_oxidation = %d\n",PCCL_soot_oxidation);CX_Message("TAR_oxidation = %d\n",TAR_oxidation);CX_Message("MGAS_gas_phase_oxidation = %d\n",MGAS_gas_phase_oxidation);CX_Message("fc_ar = %f\n",fc_ar);CX_Message("vm_ar = %f\n",vm_ar);CX_Message("ash_ar = %f\n",ash_ar);CX_Message("moist_ar = %f\n",moist_ar);CX_Message("a1_devolatilization = %f\n",A1_devolatilization);CX_Message("e1_devolatilization = %f\n",E1_devolatilization);CX_Message("a2_devolatilization = %f\n",A2_devolatilization);CX_Message("e2_devolatilization = %f\n",E2_devolatilization);CX_Message("a_tar_cracking = %f\n",A_tar_cracking);CX_Message("e_tar_cracking = %f\n",E_tar_cracking);CX_Message("A_steam_gasification = %f\n",A_steam_gasification);CX_Message("E_steam_gasification = %f\n",E_steam_gasification);CX_Message("K_steam_gasification = %f\n",K_steam_gasification);CX_Message("N_steam_gasification = %f\n",N_steam_gasification);CX_Message("Annealing_steam_gasification = %f\n",Annealing_steam_gasification);CX_Message("A_co2_gasification = %f\n",A_co2_gasification);CX_Message("E_co2_gasification = %f\n",E_co2_gasification);CX_Message("K_co2_gasification = %f\n",K_co2_gasification);CX_Message("N_co2_gasification = %f\n",N_co2_gasification);CX_Message("Annealing_co2_gasification = %f\n",Annealing_co2_gasification);CX_Message("A_h2_gasification = %f\n",A_h2_gasification);CX_Message("E_h2_gasification = %f\n",E_h2_gasification);CX_Message("N_h2_gasification = %f\n",N_h2_gasification);CX_Message("Annealing_h2_gasification = %f\n",Annealing_h2_gasification);CX_Message("A_soot_steam_gasification = %f\n",A_soot_steam_gasification);CX_Message("E_soot_steam_gasification = %f\n",E_soot_steam_gasification);CX_Message("K_soot_steam_gasification = %f\n",K_soot_steam_gasification);CX_Message("N_soot_steam_gasification = %f\n",N_soot_steam_gasification);CX_Message("Annealing_soot_steam_gasification= %f\n",Annealing_soot_steam_gasification);CX_Message("A_soot_co2_gasification = %f\n",A_soot_co2_gasification);CX_Message("E_soot_co2_gasification = %f\n",E_soot_co2_gasification);CX_Message("K_soot_co2_gasification = %f\n",K_soot_co2_gasification);CX_Message("N_soot_co2_gasification = %f\n",N_soot_co2_gasification);CX_Message("Annealing_soot_co2_gasification= %f\n",Annealing_soot_co2_gasification);CX_Message("A_soot_h2_gasification = %f\n",A_soot_h2_gasification);CX_Message("E_soot_h2_gasification = %f\n",E_soot_h2_gasification);CX_Message("N_soot_h2_gasification = %f\n",N_soot_h2_gasification);CX_Message("Annealing_soot_h2_gasification = %f\n",Annealing_soot_h2_gasification);CX_Message("A_soot_combustion = %f\n",A_Soot_Combustion);CX_Message("E_soot_combustion = %f\n",E_Soot_Combustion);CX_Message("A_c_combustion = %f\n",A_c_combustion);CX_Message("E_c_combustion = %f\n",E_c_combustion);CX_Message("N_c_combustion = %f\n",N_c_combustion);CX_Message("Annealing_c_combustion = %f\n",Annealing_c_combustion);CX_Message("A_moisture_release = %f\n",A_moisture_release);CX_Message("E_moisture_release = %f\n",E_moisture_release);CX_Message("wg3 = %f\n",wg3);CX_Message("moisture_flux = %f\n",Moisture_Flux);/* f_ep_a is used in shrinking core model, in coal combustion model */double ep_a = 0.25 + 0.75*(1-ash_ar/100.);f_ep_a = pow(ep_a,2.5);}/* pan c3m start */static void SetBooleanValue(char * var , char * svalue) /* pan : Oct 2012 ... new function */ {cxboolean value;if ( strcmp(svalue,"true") == 0) /* pan2 : Oct 2012 : correction */value = TRUE;elsevalue = FALSE;if (strcmp(var,"pccl_devol") == 0) PCCL_Devol = value;if (strcmp(var,"mgas_devol") == 0) MGAS_Devol = value;if (strcmp(var,"cpd_devol") == 0) CPD_Devol = value;if (strcmp(var,"fgdvc_devol") == 0) FGDVC_Devol = value;if (strcmp(var,"hptr_devol") == 0) HPTR_Devol = value;if (strcmp(var,"mgas_moisture") == 0) MGAS_Moisture = value;if (strcmp(var,"pccl_moisture") == 0) PCCL_Moisture = value;if (strcmp(var,"mgas_tarcracking") == 0) MGAS_TarCracking = value;if (strcmp(var,"pccl_2nd_pyro") == 0) PCCL_2nd_Pyro = value;if (strcmp(var,"mgas_gasif") == 0) MGAS_Gasif = value;if (strcmp(var,"pccl_gasif") == 0) PCCL_Gasif = value;if (strcmp(var,"pccl_tarcracking") == 0) PCCL_TarCracking = value;if (strcmp(var,"mgas_wgs") == 0) MGAS_WGS = value;if (strcmp(var,"pccl_soot_gasif") == 0) PCCL_soot_gasif = value;if (strcmp(var,"mgas_char_combustion") == 0) MGAS_char_combustion = value;if (strcmp(var,"pccl_char_combustion") == 0) PCCL_char_combustion = value;if (strcmp(var,"pccl_soot_oxidation") == 0) PCCL_soot_oxidation = value;if (strcmp(var,"tar_oxidation") == 0) TAR_oxidation = value;if (strcmp(var,"mgas_gas_phase_oxidation") == 0) MGAS_gas_phase_oxidation = value;}static void SetValue(char * var , char * svalue) /* pan : oct 2012 ... replace entire function */{char * pEnd;double value = strtod(svalue,&pEnd);if (strcmp(var,"fc_ar") == 0) fc_ar = value;if (strcmp(var,"vm_ar") == 0) vm_ar = value;if (strcmp(var,"ash_ar") == 0) ash_ar = value;if (strcmp(var,"moist_ar") == 0) moist_ar = value;/* NOTE:If the two stage devolatilization model is not used, A2 and E2 have to be set to zero*/if (strcmp(var,"a1_devolatilization") == 0) A1_devolatilization = value;if (strcmp(var,"e1_devolatilization") == 0) E1_devolatilization = value;if (strcmp(var,"a2_devolatilization") == 0) A2_devolatilization = value;if (strcmp(var,"e2_devolatilization") == 0) E2_devolatilization = value;/* if (strcmp(var,"c3m_aem") == 0) c3m_aem = value; should not be used pan : Oct 2012 */if (strcmp(var,"a_tar_cracking") == 0) A_tar_cracking = value;if (strcmp(var,"e_tar_cracking") == 0) E_tar_cracking = value;if (strcmp(var,"a_steam_gasification") == 0) A_steam_gasification = value;if (strcmp(var,"e_steam_gasification") == 0) E_steam_gasification = value;if (strcmp(var,"annealing_steam_gasification") == 0) Annealing_steam_gasification = value;if (strcmp(var,"k_steam_gasification") == 0) K_steam_gasification = value;if (strcmp(var,"n_steam_gasification") == 0) N_steam_gasification = value;if (strcmp(var,"a_co2_gasification") == 0) A_co2_gasification = value;if (strcmp(var,"e_co2_gasification") == 0) E_co2_gasification = value;if (strcmp(var,"annealing_co2_gasification") == 0) Annealing_co2_gasification = value;if (strcmp(var,"k_co2_gasification") == 0) K_co2_gasification = value;if (strcmp(var,"n_co2_gasification") == 0) N_co2_gasification = value;if (strcmp(var,"a_h2_gasification") == 0) A_h2_gasification = value;if (strcmp(var,"e_h2_gasification") == 0) E_h2_gasification = value;if (strcmp(var,"annealing_h2_gasification") == 0) Annealing_h2_gasification = value;if (strcmp(var,"n_h2_gasification") == 0) N_h2_gasification = value;if (strcmp(var,"a_soot_steam_gasification") == 0) A_soot_steam_gasification = value;if (strcmp(var,"e_soot_steam_gasification") == 0) E_soot_steam_gasification = value;if (strcmp(var,"annealing_soot_steam_gasification") == 0) Annealing_soot_steam_gasification = value;if (strcmp(var,"k_soot_steam_gasification") == 0) K_soot_steam_gasification = value;if (strcmp(var,"n_soot_steam_gasification") == 0) N_soot_steam_gasification = value;if (strcmp(var,"a_soot_co2_gasification") == 0) A_soot_co2_gasification = value;if (strcmp(var,"e_soot_co2_gasification") == 0) E_soot_co2_gasification = value;if (strcmp(var,"annealing_soot_co2_gasification") == 0) Annealing_soot_co2_gasification = value;if (strcmp(var,"k_soot_co2_gasification") == 0) K_soot_co2_gasification = value;if (strcmp(var,"n_soot_co2_gasification") == 0) N_soot_co2_gasification= value;if (strcmp(var,"a_soot_h2_gasification") == 0) A_soot_h2_gasification = value;if (strcmp(var,"e_soot_h2_gasification") == 0) E_soot_h2_gasification = value;if (strcmp(var,"annealing_soot_h2_gasification") == 0) Annealing_soot_h2_gasification = value;if (strcmp(var,"n_soot_h2_gasification") == 0) N_soot_h2_gasification = value;if (strcmp(var,"a_soot_combustion") == 0) A_Soot_Combustion = value;if (strcmp(var,"e_soot_combustion") == 0) E_Soot_Combustion = value;if (strcmp(var,"a_c_combustion") == 0) A_c_combustion = value;if (strcmp(var,"e_c_combustion") == 0) E_c_combustion = value;if (strcmp(var,"n_c_combustion") == 0) N_c_combustion = value;if (strcmp(var,"annealing_c_combustion") == 0) Annealing_c_combustion = value;if (strcmp(var,"a_moisture_release") == 0) A_moisture_release = value;if (strcmp(var,"e_moisture_release") == 0) E_moisture_release = value;if (strcmp(var,"c3m_wg3") == 0) wg3 = value;if (strcmp(var,"moisture_flux") == 0) Moisture_Flux = value;/* if (strcmp(var,"c3m_ae5") == 0) c3m_ae5 = value; pan : Oct 2012 ... should not be used */}void read_c3m_data(){CX_Message("start of read_c3m_data \n");FILE * pFile;char line[80];char field1[80];char field2[80];char *pch;int i , field_index;pFile = fopen("fluent_c3m_udf.inp","r");if (pFile != NULL){CX_Message("fopen OK \n");char * pEnd = line;while (pEnd != NULL){/* CX_Message("read a line - start \n"); */pEnd = fgets(line,80,pFile);if (pEnd != NULL){for (i=0; i<80; ++i) line[i] = tolower(line[i]);pch = strtok(line," ,\t\n=");field_index = 0;while (pch != NULL){if (field_index == 0){strcpy(field1,pch);field_index = 1;/* CX_Message("token %s \n",field1); */}else{strcpy(field2,pch);field_index = 0;SetValue(field1,field2);SetBooleanValue(field1,field2); /* pan : Oct 2012 *//*CX_Message("\n\n");CX_Message("token %s \n",field1);CX_Message("token %s \n",field2);*/}pch = strtok(NULL," ,\t\n=");}}}fclose(pFile);}。
