Skip to content

Commit

Permalink
Site updated: 2024-10-09 22:28:20
Browse files Browse the repository at this point in the history
  • Loading branch information
@grabtheNS
grabtheNS committed Oct 9, 2024
1 parent 9d28a38 commit 8cb9a3d
Show file tree
Hide file tree
Showing 4 changed files with 25 additions and 44 deletions.
63 changes: 24 additions & 39 deletions 2024/10/07/在OpenFOAM中构建SE双方程模型/index.html
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,10 +26,12 @@
<meta property="og:site_name" content="jiaqianzhang">
<meta property="og:description" content="SE双方程模型SE双方程模型&#x3D;单分布模型，它是一种简化的模型假设，假定碳烟的粒径分布为单分布。具体的模型假定可以参考SootLib: Monodispersed (byu.edu)。 值得说明的是，目前的SootLib仍然有很多的BUG，目前通过我的测试的搭配： SootModelSelect{ nucleationModel LL; growthModel LL;">
<meta property="og:locale" content="zh_CN">
<meta property="og:image" content="http://example.com/images/%E8%BE%90%E5%B0%84%E5%90%B8%E6%94%B6%E7%B3%BB%E6%95%B0%E6%8B%9F%E5%90%88%E6%9B%B2%E7%BA%BF.png">
<meta property="article:published_time" content="2024-10-07T14:57:59.000Z">
<meta property="article:modified_time" content="2024-10-08T08:36:48.597Z">
<meta property="article:modified_time" content="2024-10-09T14:26:12.412Z">
<meta property="article:author" content="jiaqianzhang">
<meta name="twitter:card" content="summary_large_image">
<meta name="twitter:image" content="http://example.com/images/%E8%BE%90%E5%B0%84%E5%90%B8%E6%94%B6%E7%B3%BB%E6%95%B0%E6%8B%9F%E5%90%88%E6%9B%B2%E7%BA%BF.png">



Expand Down Expand Up @@ -235,7 +237,7 @@
<span class="post-meta mr-2">
<i class="iconfont icon-chart"></i>

739 字
1.1k 字

</span>

Expand All @@ -246,7 +248,7 @@



7 分钟
10 分钟

</span>

Expand Down Expand Up @@ -301,52 +303,35 @@ <h1 id="SE双方程模型"><a href="#SE双方程模型" class="headerlink" title
<p>值得说明的是，目前的SootLib仍然有很多的BUG，目前通过我的测试的搭配：</p>
<p>SootModelSelect<br>{<br> nucleationModel LL;<br> growthModel LL;<br> oxidationModel OPTG;<br> coagulationModel NONE;<br>}</p>
<h1 id="具体实现耦合框架"><a href="#具体实现耦合框架" class="headerlink" title="具体实现耦合框架"></a>具体实现耦合框架</h1><h2 id="输运方程"><a href="#输运方程" class="headerlink" title="输运方程"></a>输运方程</h2><p>在OpenFOAM中，求解数密度矩的输运方程：</p>
<p> <code>fvScalarMatrix MOM0Eqn</code></p>
<p>​ <code>(</code> </p>
<p>​ <code>fvm::ddt(rho , M0) +fvm::div(phi, M0 ) \+ fvm::div(phiTh, M0)// 添加热泳项</code> </p>
<p>​ <code>==</code> </p>
<p>​ <code>fvm::laplacian(diffusionCoeff * rhoGas , M0 ) + SM0</code> </p>
<p>​ <code>);</code></p>
<p>​ <code>fvScalarMatrix MOM1Eqn</code></p>
<p>​ <code>(</code></p>
<p>​ <code>fvm::ddt(rho , M1) +fvm::div(phi, M1) \+ fvm::div(phiTh, M1) // 添加热泳项</code></p>
<p>​ <code>==</code></p>
<p>​ <code>fvm::laplacian( diffusionCoeff * rhoGas, M1 ) +SM1</code> </p>
<p>​ <code>);</code></p>
<figure class="highlight less"><table><tr><td class="gutter"><pre><span class="line">1</span><br><span class="line">2</span><br><span class="line">3</span><br><span class="line">4</span><br><span class="line">5</span><br><span class="line">6</span><br><span class="line">7</span><br><span class="line">8</span><br><span class="line">9</span><br><span class="line">10</span><br><span class="line">11</span><br><span class="line">12</span><br><span class="line">13</span><br></pre></td><td class="code"><pre><code class="hljs less"> <span class="hljs-selector-tag">fvScalarMatrix</span> <span class="hljs-selector-tag">MOM0Eqn</span><br><br>​ ( <br><br>​ <span class="hljs-attribute">fvm</span>::<span class="hljs-built_in">ddt</span>(rho , M0) +<span class="hljs-attribute">fvm</span>::<span class="hljs-built_in">div</span>(phi, M0 ) \+ <span class="hljs-attribute">fvm</span>::<span class="hljs-built_in">div</span>(phiTh, M0)<span class="hljs-comment">// 添加热泳项 </span><br><br>​ == <span class="hljs-attribute">fvm</span>::<span class="hljs-built_in">laplacian</span>(diffusionCoeff * rhoGas , M0 ) + SM0 );<br><br>​ <span class="hljs-selector-tag">fvScalarMatrix</span> <span class="hljs-selector-tag">MOM1Eqn</span><br><br>​ (<span class="hljs-attribute">fvm</span>::<span class="hljs-built_in">ddt</span>(rho , M1) +<span class="hljs-attribute">fvm</span>::<span class="hljs-built_in">div</span>(phi, M1) \+ <span class="hljs-attribute">fvm</span>::<span class="hljs-built_in">div</span>(phiTh, M1) <span class="hljs-comment">// 添加热泳项</span><br><br>​ ==<span class="hljs-attribute">fvm</span>::<span class="hljs-built_in">laplacian</span>( diffusionCoeff * rhoGas, M1 ) +SM1 );<br></code></pre></td></tr></table></figure>
<p>其中第一项是时间导数项，第二项是对流项，第三项是热泳项，第四项是扩散项，第五项是源项<br>其中，相关文献提到扩散项可以忽略不计，而热泳项在AMREX中并未考虑。具体的建模公式可以参考文献。</p>
<p>这里我是直接进行求解，有文献（NGA求解器）采用分裂算子的方法进行计算。</p>
<p>源项的计算：</p>
<p>​ <code>std::vector&lt;double&gt; yGas&#123;Y_H[cellI], Y_H2[cellI], Y_O[cellI], Y_O2[cellI], Y_OH[cellI],Y_H2O[cellI],Y_CO[cellI], Y_C2H2[cellI] &#125;;</code></p>
<p>​ <code>std::vector&lt;double&gt; Msoot&#123;M0[cellI]* rhoGas[cellI] , M1[cellI]* rhoGas[cellI] , M2[cellI]* rhoGas[cellI], M3[cellI]* rhoGas[cellI]&#125;; // soot moment values [M0, M1, M2, M3]</code></p>
<p>​ <code>S.setState(T[cellI], p[cellI], rhoGas[cellI], muGas[cellI], yGas, yPAH, Msoot, nsoot );</code></p>
<figure class="highlight css"><table><tr><td class="gutter"><pre><span class="line">1</span><br><span class="line">2</span><br><span class="line">3</span><br><span class="line">4</span><br><span class="line">5</span><br></pre></td><td class="code"><pre><code class="hljs css"> `std::vector&lt;double&gt; yGas&#123;Y_H<span class="hljs-selector-attr">[cellI]</span>, Y_H2<span class="hljs-selector-attr">[cellI]</span>, Y_O<span class="hljs-selector-attr">[cellI]</span>, Y_O2<span class="hljs-selector-attr">[cellI]</span>, Y_OH<span class="hljs-selector-attr">[cellI]</span>,Y_H2O<span class="hljs-selector-attr">[cellI]</span>,Y_CO<span class="hljs-selector-attr">[cellI]</span>, Y_C2H2<span class="hljs-selector-attr">[cellI]</span> &#125;;`<br><br>​ `std::vector&lt;double&gt; Msoot&#123;M0<span class="hljs-selector-attr">[cellI]</span>* rhoGas<span class="hljs-selector-attr">[cellI]</span> , M1<span class="hljs-selector-attr">[cellI]</span>* rhoGas<span class="hljs-selector-attr">[cellI]</span> , M2<span class="hljs-selector-attr">[cellI]</span>* rhoGas<span class="hljs-selector-attr">[cellI]</span>, M3<span class="hljs-selector-attr">[cellI]</span>* rhoGas<span class="hljs-selector-attr">[cellI]</span>&#125;; // soot moment values <span class="hljs-selector-attr">[M0, M1, M2, M3]</span>`<br><br>​ `S<span class="hljs-selector-class">.setState</span>(T<span class="hljs-selector-attr">[cellI]</span>, <span class="hljs-selector-tag">p</span><span class="hljs-selector-attr">[cellI]</span>, rhoGas<span class="hljs-selector-attr">[cellI]</span>, muGas<span class="hljs-selector-attr">[cellI]</span>, yGas, yPAH, Msoot, nsoot );<br></code></pre></td></tr></table></figure>
<p>这里大致的求解步骤：<br>给一个小的M1 M0初值，然后与其它所需化学物质温度等变量作为输入进行矩的源项计算，然后进行矩的输运方程计算。</p>
<h2 id="耦合物种消耗"><a href="#耦合物种消耗" class="headerlink" title="耦合物种消耗"></a>耦合物种消耗</h2><p>考虑碳烟带来的物质消耗，以02为例<br><code>if (Yi.name() == &quot;O2&quot;) &#123;</code></p>
<p>​ <code>YiEqn -= O2ConsumptionRate;</code></p>
<h2 id="耦合物种消耗"><a href="#耦合物种消耗" class="headerlink" title="耦合物种消耗"></a>耦合物种消耗</h2><p>考虑碳烟带来的物质消耗，以02为例</p>
<p>在物种的输运方程添加一个源项。</p>
<h2 id="修正质量守恒方程"><a href="#修正质量守恒方程" class="headerlink" title="修正质量守恒方程"></a>修正质量守恒方程</h2><p>同时也要考虑碳烟带来的气相质量损失，修改连续性方程 fvScalarMatrix rhoEqn</p>
<p> <code>(</code></p>
<p>​ <code>fvm::ddt(rho)</code></p>
<p> <code>\+ fvc::div(phi)</code></p>
<p> <code>==</code></p>
<p> <code>totalConsumptionRate</code></p>
<p> <code>);</code><br>其中，totalConsumptionRate为所有化学物质的总消耗</p>
<h2 id="修正过氧化和非物理结果"><a href="#修正过氧化和非物理结果" class="headerlink" title="修正过氧化和非物理结果"></a>修正过氧化和非物理结果</h2><p>​ <code>forAll(mesh.cells(), cellI)</code></p>
<p>​ <code>&#123;</code></p>
<p>​ <code>if( M0[cellI] &lt; smallWeight || M1[cellI] &lt; tolV )</code></p>
<p>​ <code>&#123;</code></p>
<p>​ <code>M1[cellI] = std::max(M1[cellI] , tolV);</code></p>
<p>​ <code>M0[cellI] = M1[cellI] /nuclVol ;</code></p>
<p>​ <code>Info &lt;&lt; &quot;modify in &quot; &lt;&lt; cellI&lt;&lt;endl;</code></p>
<p>​ <code>&#125;</code></p>
<p>​ <code>sootfv[cellI] = M1[cellI] * rho[cellI] /1850 ; //计算碳烟体积分数</code></p>
<p>​ <code>&#125;</code></p>
<p>在结束计算后，由于实际物理意义中矩不可能有负数，而计算过程中由于过氧化带来的负数矩要进行修正，用一个较小的数字进行替代。</p>
<p>其中，totalConsumptionRate为所有化学物质的总消耗`</p>
<h2 id="修正过氧化和非物理结果"><a href="#修正过氧化和非物理结果" class="headerlink" title="修正过氧化和非物理结果"></a>修正过氧化和非物理结果</h2><p>在结束计算后，由于实际物理意义中矩不可能有负数，而计算过程中由于过氧化带来的负数矩要进行修正，用一个较小的数字进行替代。</p>
<h2 id="碳烟ppm转化公式"><a href="#碳烟ppm转化公式" class="headerlink" title="碳烟ppm转化公式"></a>碳烟ppm转化公式</h2><p>求得的M1还不是我们常用的碳烟体积分数，假定碳烟的密度为1850g/cm^3<br>sootfv[cellI] = M1[cellI] * rho[cellI] /1850</p>
<h2 id="压力边界条件设置"><a href="#压力边界条件设置" class="headerlink" title="压力边界条件设置"></a>压力边界条件设置</h2><p>不要设置fixValue边界条件，这样改变燃烧室的压力时候会计算报错。正确的压力边界条件设置：</p>
<p>dimensions [1 -1 -2 0 0 0 0];</p>
<p>internalField uniform 201325;</p>
<p>boundaryField<br>{<br> inletfuel<br> {<br> type zeroGradient;<br> }<br> inletair<br> {<br> type zeroGradient;<br> }<br> outlet<br> {<br> type zeroGradient;<br> }<br> axis<br> {<br> type empty;<br> }<br> leftside<br> {<br> type zeroGradient;<br> }<br> burnerwall<br> {<br> type zeroGradient;<br> }<br> front<br> {<br> type wedge;<br> }<br> back<br> {<br> type wedge;<br> }<br>}</p>
<h2 id="光学薄膜假设辐射项"><a href="#光学薄膜假设辐射项" class="headerlink" title="光学薄膜假设辐射项"></a>光学薄膜假设辐射项</h2><p>理论基础：</p>
<p>全部改为zeroGradient条件</p>
<h2 id="温度边界条件"><a href="#温度边界条件" class="headerlink" title="温度边界条件"></a>温度边界条件</h2><p>杜博建议：将碰嘴改为绝热边界条件可以减少举升</p>
<h2 id="光学薄膜假设辐射项"><a href="#光学薄膜假设辐射项" class="headerlink" title="光学薄膜假设辐射项"></a>光学薄膜假设辐射项</h2><p>理论基础：<br>参考Maciej Chmielewski等人的文章《Planck Mean Absorption Coefficients of H2O, CO2, CO and NO for radiation numerical modeling in combusting flows》</p>
<p>所谓光学薄膜假设实际上就是在计算能量方程的右边添加一个辐射源项</p>
<script type="math/tex; mode=display">
q_{\text {loss }}=\sum_{i=1}^N \frac{x_i}{x_{\text {tot }}} p_{\text {tot }} \kappa_i \sigma\left(T^4-T_b^4\right)</script><p>其中，第一项是需要计算的物种摩尔浓度（无量纲），第二个是总压（atm),第三个是普朗克平均吸收系数（m pa)-1 Tb是环境温度，这里我设置为300K。</p>
<p>碳烟的拟合是通过：</p>
<script type="math/tex; mode=display">
Q_{\text {soot }}=C_{\text {soot }} \cdot f_v \cdot T^5</script><p>其中Csoot为</p>
<script type="math/tex; mode=display">
C_{\text {soot }}=3.334 \cdot 10^{-4} \mathrm{~W} /\left(\mathrm{m}^3 \mathrm{~K}^5\right)</script><p>fv为碳烟的体积分数</p>
<h3 id="普朗克平均吸收系数（m-pa-1-的拟合"><a href="#普朗克平均吸收系数（m-pa-1-的拟合" class="headerlink" title="普朗克平均吸收系数（m pa)-1 的拟合"></a>普朗克平均吸收系数（m pa)-1 的拟合</h3><p>这里通过6次多项式进行拟合，数据来源是PELELMX提供的。</p>
<p><img src="/images/辐射吸收系数拟合曲线.png" srcset="/img/loading.gif" lazyload alt="辐射吸收系数拟合曲线"></p>
<figure class="highlight inform7"><table><tr><td class="gutter"><pre><span class="line">1</span><br><span class="line">2</span><br><span class="line">3</span><br><span class="line">4</span><br><span class="line">5</span><br><span class="line">6</span><br><span class="line">7</span><br><span class="line">8</span><br><span class="line">9</span><br><span class="line">10</span><br><span class="line">11</span><br><span class="line">12</span><br><span class="line">13</span><br><span class="line">14</span><br><span class="line">15</span><br><span class="line">16</span><br><span class="line">17</span><br><span class="line">18</span><br><span class="line">19</span><br><span class="line">20</span><br><span class="line">21</span><br><span class="line">22</span><br><span class="line">23</span><br><span class="line">24</span><br><span class="line">25</span><br><span class="line">26</span><br><span class="line">27</span><br></pre></td><td class="code"><pre><code class="hljs inform7">double calculate_ap(double T, const std::vector&lt;double&gt;&amp; coeffs) &#123;<br> // T 是温度，coeffs 是存储 c0, c1, c2, c3, c4, c5 的系数组<br> double T_inv = T; // <br> double ap = coeffs<span class="hljs-comment">[0]</span> <br> + coeffs<span class="hljs-comment">[1]</span> * T_inv <br> + coeffs<span class="hljs-comment">[2]</span> * std::pow(T_inv, 2) <br> + coeffs<span class="hljs-comment">[3]</span> * std::pow(T_inv, 3) <br> + coeffs<span class="hljs-comment">[4]</span> * std::pow(T_inv, 4) <br> + coeffs<span class="hljs-comment">[5]</span> * std::pow(T_inv, 5)<br> + coeffs<span class="hljs-comment">[6]</span> * std::pow(T_inv, 6);<br> return ap;<br> &#125;<br> if(radiationButton =1)&#123;<br> Qsoot<span class="hljs-comment">[cellI]</span> = Csoot<span class="hljs-comment">[cellI]</span> * pow(T<span class="hljs-comment">[cellI]</span>, 5) * (sootfv<span class="hljs-comment">[cellI]</span>);<br> absorptivity_CO2<span class="hljs-comment">[cellI]</span> = calculate_ap(T<span class="hljs-comment">[cellI]</span>, CO2_coeffs); //注意这里吸收系数算出来单位/(cm*bar)-1<br> absorptivity_CH4<span class="hljs-comment">[cellI]</span> = calculate_ap(T<span class="hljs-comment">[cellI]</span>, CH4_coeffs);<br> absorptivity_CO<span class="hljs-comment">[cellI]</span> = calculate_ap(T<span class="hljs-comment">[cellI]</span>, CO_coeffs);<br> absorptivity_H2O<span class="hljs-comment">[cellI]</span> = calculate_ap(T<span class="hljs-comment">[cellI]</span>, H2O_coeffs);<br> //计算公式：kpl系数*总压（atm)*摩尔分数 *玻尔兹曼常数* T相关 气体的总摩尔质量通过理想气体状态方程进行换算<br> mean_MW<span class="hljs-comment">[cellI]</span> = rho<span class="hljs-comment">[cellI]</span>/ p<span class="hljs-comment">[cellI]</span>* 8.314 * T<span class="hljs-comment">[cellI]</span>;<br> QgasRadLoss<span class="hljs-comment">[cellI]</span> = (p<span class="hljs-comment">[cellI]</span> / 101325) * (<br> (Y_H2O<span class="hljs-comment">[cellI]</span> / 18 * mean_MW<span class="hljs-comment">[cellI]</span> * absorptivity_H2O<span class="hljs-comment">[cellI]</span>) + <br> (Y_CO2<span class="hljs-comment">[cellI]</span> / 44 * mean_MW<span class="hljs-comment">[cellI]</span> * absorptivity_CO2<span class="hljs-comment">[cellI]</span>) + <br> (Y_CO<span class="hljs-comment">[cellI]</span> / 28 * mean_MW<span class="hljs-comment">[cellI]</span> * absorptivity_CO<span class="hljs-comment">[cellI]</span>) + <br> (Y_CH4<span class="hljs-comment">[cellI]</span> / 16 * mean_MW<span class="hljs-comment">[cellI]</span> * absorptivity_CH4<span class="hljs-comment">[cellI]</span>)<br> ) * 100 * (5.67e-8) * (pow4(T<span class="hljs-comment">[cellI]</span>) - pow4(Tb<span class="hljs-comment">[cellI]</span>));&#125;<br><br></code></pre></td></tr></table></figure>


</div>
Expand Down
4 changes: 0 additions & 4 deletions about/index.md
View file

This file was deleted.

Binary file added images/辐射吸收系数拟合曲线.png
View file
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Loading

0 comments on commit 8cb9a3d

Please sign in to comment.