GeezerVisonThe World Wide Web, restored to 800 pixels of glory. |
| Home|About|Gaming cabinet|Archives|How it works|Starter sites |
GeezerVison is showing github.com — MFlowCode.github.io/documentation/equations.html at master · MFlowCode/MFlowCode.github.io Expand file tree/ Copy pathequations.htmlMore file actions 1239 lines (1237 loc) · 91.3 KB / Copy pathequations.htmlFile metadata and controls1239 lines (1237 loc) · 91.3 KB 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703 704 705 706 707 708 709 710 711 712 713 714 715 716 717 718 719 720 721 722 723 724 725 726 727 728 729 730 731 732 733 734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750 751 752 753 754 755 756 757 758 759 760 761 762 763 764 765 766 767 768 769 770 771 772 773 774 775 776 777 778 779 780 781 782 783 784 785 786 787 788 789 790 791 792 793 794 795 796 797 798 799 800 801 802 803 804 805 806 807 808 809 810 811 812 813 814 815 816 817 818 819 820 821 822 823 824 825 826 827 828 829 830 831 832 833 834 835 836 837 838 839 840 841 842 843 844 845 846 847 848 849 850 851 852 853 854 855 856 857 858 859 860 861 862 863 864 865 866 867 868 869 870 871 872 873 874 875 876 877 878 879 880 881 882 883 884 885 886 887 888 889 890 891 892 893 894 895 896 897 898 899 900 901 902 903 904 905 906 907 908 909 910 911 912 913 914 915 916 917 918 919 920 921 922 923 924 925 926 927 928 929 930 931 932 933 934 935 936 937 938 939 940 941 942 943 944 945 946 947 948 949 950 951 952 953 954 955 956 957 958 959 960 961 962 963 964 965 966 967 968 969 970 971 972 973 974 975 976 977 978 979 980 981 982 983 984 985 986 987 988 989 990 991 992 993 994 995 996 997 998 999 1000 <!-- HTML header for doxygen 1.9.1--> <!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN" "https://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd"> <html xmlns="http://www.w3.org/1999/xhtml"> <head> <!-- Google tag (gtag.js) --> <script async src="https://www.googletagmanager.com/gtag/js?id=G-SY496B9L99"></script> <script> window.dataLayer = window.dataLayer || []; function gtag(){dataLayer.push(arguments);} gtag('js', new Date()); gtag('config', 'G-SY496B9L99'); </script> <meta http-equiv="Content-Type" content="text/xhtml;charset=UTF-8"/> <meta http-equiv="X-UA-Compatible" content="IE=9"/> <meta name="generator" content="Doxygen 1.16.1"/> <meta name="viewport" content="width=device-width, initial-scale=1"/> <title>MFC: Equations</title> <meta name="description" content="Equations — MFC documentation. Open-source exascale multiphase flow solver." /> <meta name="keywords" content="exascale, fluid dynamics, cfd, computational fluid dynamics, compressible, hpc, bryngelson, colonius, subgrid, multiphase, frontier, summit, el capitan, aurora, amd gpu, gpu, nvidia"/> <link href="tabs.css" rel="stylesheet" type="text/css"/> <script type="text/javascript" src="jquery.js"></script> <script type="text/javascript" src="dynsections.js"></script> <link href="navtree.css" rel="stylesheet" type="text/css"/> <script type="text/javascript" src="navtreedata.js"></script> <script type="text/javascript" src="navtree.js"></script> <script type="text/javascript" src="cookie.js"></script> <link href="search/search.css" rel="stylesheet" type="text/css"/> <script type="text/javascript" src="search/searchdata.js"></script> <script type="text/javascript" src="search/search.js"></script> <script type="text/x-mathjax-config"> MathJax.Hub.Config({ extensions: ["tex2jax.js", "TeX/AMSmath.js", "TeX/AMSsymbols.js"], jax: ["input/TeX","output/HTML-CSS"], }); // This file is set as MATHJAX_CODEFILE in the Doxyfile. It configures how // MathJax renders expressions in Markdown so that it is consistent with GitHub. MathJax.Hub.Config({ extensions: ["tex2jax.js"], jax: ["input/TeX", "output/HTML-CSS"], tex2jax: { inlineMath: [ ['$', '$'], ["\\(","\\)"] ], displayMath: [ ['$$','$$'], ["\\[","\\]"] ], processEscapes: true, ignoreClass: "line" // Ignore code blocks: https://web.archive.org/web/20120430100225/http://www.mathjax.org/docs/1.1/options/tex2jax.html }, "HTML-CSS": { fonts: ["TeX"] } }); </script> <script type="text/javascript" async="async" src="https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.5/MathJax.js"></script> <link href="doxygen.css" rel="stylesheet" type="text/css" /> <link rel="shortcut icon" href="icon.ico" type="image/x-icon" /> <link href="doxygen-awesome.css" rel="stylesheet" type="text/css"/> <link href="doxygen-awesome-sidebar-only.css" rel="stylesheet" type="text/css"/> <link href="custom.css" rel="stylesheet" type="text/css"/> </head> <body> <div id="top"><!-- do not remove this div, it is closed by doxygen! --> <div id="titlearea"> <table cellspacing="0" cellpadding="0"> <tbody> <tr style="height: 56px;"> <td id="projectlogo"><img alt="Logo" src="icon.ico"/></td> <td id="projectalign" style="padding-left: 0.5em;"> <div id="projectname">MFC </div> <div id="projectbrief">Exascale flow solver</div> </td> </tr> </tbody> </table> </div> <!-- Cross-navigation injected into sidebar via script below --> <script> document.addEventListener('DOMContentLoaded', function() { var nav = document.createElement('div'); nav.id = 'mfc-nav'; var items = [ ['../documentation/index.html', 'documentation', 'User Guide'], ['../api/index.html', 'api', 'API Documentation'] ]; var path = window.location.pathname; var apiPaths = ['/api/', '/pre_process/', '/simulation/', '/post_process/']; for (var i = 0; i < items.length; i++) { var a = document.createElement('a'); a.href = items[i][0]; a.textContent = items[i][2]; if (items[i][1] === 'api') { for (var j = 0; j < apiPaths.length; j++) { if (path.indexOf(apiPaths[j]) !== -1) { a.className = 'active'; break; } } } else { if (path.indexOf('/' + items[i][1] + '/') !== -1) a.className = 'active'; } nav.appendChild(a); } var sideNav = document.getElementById('side-nav'); if (sideNav) sideNav.insertBefore(nav, sideNav.firstChild); }); </script> <!-- end header part --> <!-- Generated by Doxygen 1.16.1 --> <script type="text/javascript"> var searchBox = new SearchBox("searchBox", "search/",'.html'); </script> <script type="text/javascript"> $(function() { codefold.init(); }); </script> <script type="text/javascript" src="menudata.js"></script> <script type="text/javascript" src="menu.js"></script> <script type="text/javascript"> $(function() { initMenu('',true,false,'search.php','Search',true); $(function() { init_search(); }); }); </script> <div id="main-nav"></div> </div><!-- top --> <div id="side-nav" class="ui-resizable side-nav-resizable"> <div id="nav-tree"> <div id="nav-tree-contents"> <div id="nav-sync" class="sync"></div> </div> </div> <div id="splitbar" style="-moz-user-select:none;" class="ui-resizable-handle"> </div> </div> <script type="text/javascript"> $(function(){initNavTree('equations.html','',''); }); </script> <div id="container"> <div id="doc-content"> <!-- window showing the filter options --> <div id="MSearchSelectWindow" onmouseover="return searchBox.OnSearchSelectShow()" onmouseout="return searchBox.OnSearchSelectHide()" onkeydown="return searchBox.OnSearchSelectKey(event)"> </div> <!-- iframe showing the search results (closed by default) --> <div id="MSearchResultsWindow"> <div id="MSearchResults"> <div class="SRPage"> <div id="SRIndex"> <div id="SRResults"></div> <div class="SRStatus" id="Loading">Loading...</div> <div class="SRStatus" id="Searching">Searching...</div> <div class="SRStatus" id="NoMatches">No Matches</div> </div> </div> </div> </div> <div><div class="header"> <div class="headertitle"><div class="title">Equations </div></div> </div><!--header--> <div class="contents"> <div class="textblock"><h1 class="doxsection"><a class="anchor" id="autotoc_md144"></a> MFC: Comprehensive Equations Reference</h1> <p>This document catalogs every equation solved by MFC, organized by physical model. Each section notes the input parameter(s) that activate the corresponding physics module and cross-references the relevant source files.</p> <p>The models and algorithms described here are detailed in Wilfong et al. <a class="el" href="citelist.html#CITEREF_wilfong26">[68]</a> (MFC 5.0) and Bryngelson et al. <a class="el" href="citelist.html#CITEREF_bryngelson21">[8]</a>. Foundational references for each model are cited inline; see the <a class="el" href="citelist.html" title="Bibliography">Bibliography</a> for full details.</p> <p>For parameter details and allowed values, see <a class="el" href="case.html" title="Case Files">Case Files</a> and the <a class="el" href="parameters.html" title="Case Parameters">Case Parameters</a> reference.</p> <hr /> <h2 class="doxsection"><a class="anchor" id="autotoc_md146"></a> 1. Overview</h2> <p>MFC solves the compressible Navier-Stokes equations (or Euler equations when viscosity is off) in a finite volume framework. The general semi-discrete form is:</p> <p class="formulaDsp"> \[\frac{\partial \mathbf{q}}{\partial t} + \nabla \cdot \mathbf{F}(\mathbf{q}) + \mathbf{h}(\mathbf{q})\,\nabla \cdot \mathbf{u} = \mathbf{s}(\mathbf{q})\] </p> <p>where:</p><ul> <li>\(\mathbf{q}\) is the conservative variable vector,</li> <li>\(\mathbf{F}\) is the flux tensor,</li> <li>\(\mathbf{h}(\mathbf{q})\,\nabla \cdot \mathbf{u}\) contains non-conservative terms (volume fraction advection),</li> <li>\(\mathbf{s}(\mathbf{q})\) is the source vector (bubbles, body forces, chemistry, etc.).</li> </ul> <p>The parameter <span class="tt">model_eqns</span> (1, 2, or 3) selects the governing equation set.</p> <p><b>Key source files:</b> <span class="tt">src/simulation/m_rhs.fpp</span> (RHS evaluation), <span class="tt">src/common/m_eos.fpp</span> (equations of state), <span class="tt">src/common/m_variables_conversion.fpp</span> (variable conversion and mixture rules).</p> <hr /> <h2 class="doxsection"><a class="anchor" id="sec-units-dimensions"></a> 1b. Units, Dimensions, and Non-Dimensionalization</h2> <h3 class="doxsection"><a class="anchor" id="sec-dimensional-handling"></a> General Users: Dimensional Handling</h3> <h4 class="doxsection"><a class="anchor" id="sec-dimensions-in-out"></a> Dimensions In = Dimensions Out</h4> <p>The main flow solver (Navier-Stokes equations, Riemann solvers, viscous stress, body forces, surface tension, etc.) is <b>unit-agnostic</b>: whatever units the user provides for the initial and boundary conditions, the solver preserves them throughout the computation. If the user inputs SI units, the outputs are in SI units. If the user inputs CGS, the outputs are in CGS. No internal non-dimensionalization is performed by the flow solver.</p> <p>This means that for simulations <b>without</b> sub-grid bubble models, the user can work in any consistent unit system without additional effort.</p> <h4 class="doxsection"><a class="anchor" id="sec-stored-forms"></a> Stored Parameter Conventions</h4> <p>Several EOS and transport parameters use <b>transformed stored forms</b> that differ from the standard physical values. This is the most common source of input errors:</p> <table class="markdownTable"> <tr class="markdownTableHead"> <th class="markdownTableHeadNone">Parameter </th><th class="markdownTableHeadNone">Physical quantity </th><th class="markdownTableHeadNone">What MFC expects (stored form) </th></tr> <tr class="markdownTableRowOdd"> <td class="markdownTableBodyNone"><span class="tt">fluid_pp(i)%gamma</span> </td><td class="markdownTableBodyNone">Heat capacity ratio \(\gamma\) </td><td class="markdownTableBodyNone">\(\Gamma = \frac{1}{\gamma - 1}\) </td></tr> <tr class="markdownTableRowEven"> <td class="markdownTableBodyNone"><span class="tt">fluid_pp(i)%pi_inf</span> </td><td class="markdownTableBodyNone">Stiffness pressure \(\pi_\infty\) [Pa] </td><td class="markdownTableBodyNone">\(\Pi_\infty = \frac{\gamma\,\pi_\infty}{\gamma - 1}\) [Pa] </td></tr> <tr class="markdownTableRowOdd"> <td class="markdownTableBodyNone"><span class="tt">fluid_pp(i)%Re(1)</span> </td><td class="markdownTableBodyNone">Dynamic viscosity \(\mu\) </td><td class="markdownTableBodyNone">\(1/\mu\) (inverse viscosity) </td></tr> <tr class="markdownTableRowEven"> <td class="markdownTableBodyNone"><span class="tt">fluid_pp(i)%Re(2)</span> </td><td class="markdownTableBodyNone">Bulk viscosity \(\mu_b\) </td><td class="markdownTableBodyNone">\(1/\mu_b\) (inverse bulk viscosity) </td></tr> </table> <p>These transformations arise because MFC internally solves the energy equation using the transformed variables \(\Gamma\) and \(\Pi_\infty\) (see Section 3.1), and the viscous stress is computed by dividing by <span class="tt">Re</span> rather than multiplying by \(\mu\).</p> <p><b>Common mistake:</b> setting <span class="tt">fluid_pp(1)%gamma = 1.4</span> for air. The correct value is <span class="tt">1.0 / (1.4 - 1.0) = 2.5</span>. Setting <span class="tt">gamma = 1.4</span> corresponds to a physical \(\gamma \approx 1.71\), which is not a standard gas.</p> <h4 class="doxsection"><a class="anchor" id="sec-material-values"></a> Common Material Values</h4> <p>Pre-computed stored-form values for common fluids (SI units):</p> <table class="markdownTable"> <tr class="markdownTableHead"> <th class="markdownTableHeadNone">Material </th><th class="markdownTableHeadNone">\(\gamma\) </th><th class="markdownTableHeadNone">\(\pi_\infty\) [Pa] </th><th class="markdownTableHeadNone"><span class="tt">gamma</span> (stored) </th><th class="markdownTableHeadNone"><span class="tt">pi_inf</span> (stored) [Pa] </th></tr> <tr class="markdownTableRowOdd"> <td class="markdownTableBodyNone">Air </td><td class="markdownTableBodyNone">1.4 </td><td class="markdownTableBodyNone">0 </td><td class="markdownTableBodyNone">2.5 </td><td class="markdownTableBodyNone">0 </td></tr> <tr class="markdownTableRowEven"> <td class="markdownTableBodyNone">Helium </td><td class="markdownTableBodyNone">5/3 </td><td class="markdownTableBodyNone">0 </td><td class="markdownTableBodyNone">1.5 </td><td class="markdownTableBodyNone">0 </td></tr> <tr class="markdownTableRowOdd"> <td class="markdownTableBodyNone">Water (Tait) </td><td class="markdownTableBodyNone">4.4 </td><td class="markdownTableBodyNone">6.0e8 </td><td class="markdownTableBodyNone">0.2941 </td><td class="markdownTableBodyNone">7.76e8 </td></tr> <tr class="markdownTableRowEven"> <td class="markdownTableBodyNone">Water (Le Metayer et al. <a class="el" href="citelist.html#CITEREF_lemetayer04">[30]</a>) </td><td class="markdownTableBodyNone">6.12 </td><td class="markdownTableBodyNone">3.43e8 </td><td class="markdownTableBodyNone">0.1953 </td><td class="markdownTableBodyNone">4.10e8 </td></tr> </table> <p>Example for an air-water simulation:</p> <div class="fragment"><div class="line"><span class="comment"># Air (fluid 1)</span></div> <div class="line">gam_a = 1.4</div> <div class="line"><span class="stringliteral">"fluid_pp(1)%gamma"</span>: 1.0 / (gam_a - 1.0), <span class="comment"># = 2.5</span></div> <div class="line"><span class="stringliteral">"fluid_pp(1)%pi_inf"</span>: 0.0,</div> <div class="line"> </div> <div class="line"><span class="comment"># Water (fluid 2)</span></div> <div class="line">gam_w = 4.4</div> <div class="line">pi_w = 6.0e8 <span class="comment"># Pa</span></div> <div class="line"><span class="stringliteral">"fluid_pp(2)%gamma"</span>: 1.0 / (gam_w - 1.0), <span class="comment"># ≈ 0.294</span></div> <div class="line"><span class="stringliteral">"fluid_pp(2)%pi_inf"</span>: gam_w * pi_w / (gam_w - 1.0), <span class="comment"># ≈ 7.76e8</span></div> </div><!-- fragment --><p>For viscous cases, provide the <b>reciprocal</b> of the dynamic viscosity:</p> <div class="fragment"><div class="line">mu = 1.002e-3 <span class="comment"># water viscosity [Pa·s]</span></div> <div class="line"><span class="stringliteral">"fluid_pp(1)%Re(1)"</span>: 1.0 / mu, <span class="comment"># ≈ 998</span></div> </div><!-- fragment --><h4 class="doxsection"><a class="anchor" id="sec-unit-consistency"></a> Unit Consistency</h4> <p>The solver does not check or convert units. All inputs must use the <b>same consistent unit system</b> (e.g., all SI or all CGS). Mixing units — for example, pressures in atmospheres with densities in kg/m³ — will produce silently incorrect results.</p> <h3 class="doxsection"><a class="anchor" id="sec-bubble-nondim"></a> Bubble Users: Non-Dimensional Framework</h3> <h4 class="doxsection"><a class="anchor" id="sec-nondim-bubble-dynamics"></a> Non-Dimensional Bubble Dynamics</h4> <p>The sub-grid bubble models (<span class="tt">bubbles_euler = .true.</span> or <span class="tt">bubbles_lagrange = .true.</span>) solve the bubble wall dynamics in <b>non-dimensional form</b>. The bubble wall pressure equation as implemented is:</p> <p class="formulaDsp"> \[p_{bw} = \left(\text{Ca} + \frac{2}{\text{We}_b\,R_0}\right)\left(\frac{R_0}{R}\right)^{3\gamma} - \text{Ca} - \frac{4\,\text{Re}_{\text{inv}}\,\dot{R}}{R} - \frac{2}{R\,\text{We}_b}\] </p> <p>Here \(R\) and \(R_0\) are non-dimensional radii (scaled by \(x_0\)), and \(\dot{R}\) is a non-dimensional wall speed (scaled by \(u_0\)); the entire bubble ODE is solved in non-dimensional variables.</p> <p>The dimensionless groups are:</p> <table class="markdownTable"> <tr class="markdownTableHead"> <th class="markdownTableHeadNone">Dimensionless group </th><th class="markdownTableHeadNone">Definition </th><th class="markdownTableHeadNone">Code variable </th><th class="markdownTableHeadNone">Computed from </th></tr> <tr class="markdownTableRowOdd"> <td class="markdownTableBodyNone">\(\text{Ca}\) (Cavitation number) </td><td class="markdownTableBodyNone">\(p_{0,\text{ref}} - p_v\) </td><td class="markdownTableBodyNone"><span class="tt">Ca</span> </td><td class="markdownTableBodyNone"><span class="tt">bub_pp%p0ref - bub_pp%pv</span> </td></tr> <tr class="markdownTableRowEven"> <td class="markdownTableBodyNone">\(\text{Eu}\) (Euler number) </td><td class="markdownTableBodyNone">\(p_{0,\text{ref}}\) </td><td class="markdownTableBodyNone"><span class="tt">Eu</span> </td><td class="markdownTableBodyNone"><span class="tt">bub_pp%p0ref</span> </td></tr> <tr class="markdownTableRowOdd"> <td class="markdownTableBodyNone">\(\text{We}_b\) (bubble Weber number) </td><td class="markdownTableBodyNone">\(1/\sigma\) </td><td class="markdownTableBodyNone"><span class="tt">Web</span> </td><td class="markdownTableBodyNone"><span class="tt">1 / bub_pp%ss</span> </td></tr> <tr class="markdownTableRowEven"> <td class="markdownTableBodyNone">\(\text{Re}_{\text{inv}}\) (inverse bubble Reynolds number) </td><td class="markdownTableBodyNone">\(\mu_l\) </td><td class="markdownTableBodyNone"><span class="tt">Re_inv</span> </td><td class="markdownTableBodyNone"><span class="tt">bub_pp%mu_l</span> </td></tr> </table> <p>Because the bubble equations use these dimensionless numbers directly, all <span class="tt">bub_pp%%</span> inputs are interpreted by the code as <b>already non-dimensional</b>. The code does <b>not</b> non-dimensionalize bubble quantities internally. Therefore, when bubbles are enabled, the simulation must be run in a <b>fully non-dimensional</b> form: <b>all</b> inputs — flow ICs/BCs, EOS parameters, domain lengths, <span class="tt">dt</span>, and <span class="tt">bub_pp%%</span> values — must be scaled with the same \((x_0, p_0, \rho_0, u_0, t_0, T_0)\) reference quantities, or the coupled solution will be physically incorrect.</p> <h4 class="doxsection"><a class="anchor" id="sec-reference-scales"></a> Reference Scales</h4> <p>When using bubble models, the user must choose reference scales and non-dimensionalize <b>all</b> inputs (flow and bubble) consistently. The standard convention used in the MFC examples is:</p> <table class="markdownTable"> <tr class="markdownTableHead"> <th class="markdownTableHeadNone">Reference quantity </th><th class="markdownTableHeadNone">Symbol </th><th class="markdownTableHeadNone">Typical choice </th></tr> <tr class="markdownTableRowOdd"> <td class="markdownTableBodyNone">Length </td><td class="markdownTableBodyNone">\(x_0\) </td><td class="markdownTableBodyNone">\(R_{0,\text{ref}}\) (reference bubble radius) </td></tr> <tr class="markdownTableRowEven"> <td class="markdownTableBodyNone">Pressure </td><td class="markdownTableBodyNone">\(p_0\) </td><td class="markdownTableBodyNone">\(p_{0,\text{ref}}\) (reference bubble pressure) </td></tr> <tr class="markdownTableRowOdd"> <td class="markdownTableBodyNone">Density </td><td class="markdownTableBodyNone">\(\rho_0\) </td><td class="markdownTableBodyNone">\(\rho_{0,\text{ref}}\) (reference liquid density) </td></tr> <tr class="markdownTableRowEven"> <td class="markdownTableBodyNone">Velocity </td><td class="markdownTableBodyNone">\(u_0\) </td><td class="markdownTableBodyNone">\(\sqrt{p_0 / \rho_0}\) (derived) </td></tr> <tr class="markdownTableRowOdd"> <td class="markdownTableBodyNone">Time </td><td class="markdownTableBodyNone">\(t_0\) </td><td class="markdownTableBodyNone">\(x_0 / u_0\) (derived) </td></tr> <tr class="markdownTableRowEven"> <td class="markdownTableBodyNone">Temperature </td><td class="markdownTableBodyNone">\(T_0\) </td><td class="markdownTableBodyNone">\(T_{0,\text{ref}}\) (reference temperature) </td></tr> </table> <h4 class="doxsection"><a class="anchor" id="sec-nondim-inputs"></a> Non-Dimensionalization of Input Parameters</h4> <p>The following table lists every <span class="tt">bub_pp%%</span> parameter and its required non-dimensionalization:</p> <table class="markdownTable"> <tr class="markdownTableHead"> <th class="markdownTableHeadNone">Parameter </th><th class="markdownTableHeadNone">Physical meaning </th><th class="markdownTableHeadNone">Non-dimensional form </th></tr> <tr class="markdownTableRowOdd"> <td class="markdownTableBodyNone"><span class="tt">bub_pp%R0ref</span> </td><td class="markdownTableBodyNone">Reference bubble radius </td><td class="markdownTableBodyNone">\(R_{0,\text{ref}} / x_0\) </td></tr> <tr class="markdownTableRowEven"> <td class="markdownTableBodyNone"><span class="tt">bub_pp%p0ref</span> </td><td class="markdownTableBodyNone">Reference bubble pressure </td><td class="markdownTableBodyNone">\(p_{0,\text{ref}} / p_0\) </td></tr> <tr class="markdownTableRowOdd"> <td class="markdownTableBodyNone"><span class="tt">bub_pp%rho0ref</span> </td><td class="markdownTableBodyNone">Reference liquid density </td><td class="markdownTableBodyNone">\(\rho_{0,\text{ref}} / \rho_0\) </td></tr> <tr class="markdownTableRowEven"> <td class="markdownTableBodyNone"><span class="tt">bub_pp%T0ref</span> </td><td class="markdownTableBodyNone">Reference temperature </td><td class="markdownTableBodyNone">\(T_{0,\text{ref}} / T_0\) (typically 1) </td></tr> <tr class="markdownTableRowOdd"> <td class="markdownTableBodyNone"><span class="tt">bub_pp%ss</span> </td><td class="markdownTableBodyNone">Surface tension \(\sigma\) </td><td class="markdownTableBodyNone">\(\sigma / (\rho_0\,x_0\,u_0^2)\) </td></tr> <tr class="markdownTableRowEven"> <td class="markdownTableBodyNone"><span class="tt">bub_pp%pv</span> </td><td class="markdownTableBodyNone">Vapor pressure </td><td class="markdownTableBodyNone">\(p_v / p_0\) </td></tr> <tr class="markdownTableRowOdd"> <td class="markdownTableBodyNone"><span class="tt">bub_pp%mu_l</span> </td><td class="markdownTableBodyNone">Liquid dynamic viscosity </td><td class="markdownTableBodyNone">\(\mu_l / (\rho_0\,x_0\,u_0)\) </td></tr> <tr class="markdownTableRowEven"> <td class="markdownTableBodyNone"><span class="tt">bub_pp%mu_v</span> </td><td class="markdownTableBodyNone">Vapor dynamic viscosity </td><td class="markdownTableBodyNone">\(\mu_v / (\rho_0\,x_0\,u_0)\) </td></tr> <tr class="markdownTableRowOdd"> <td class="markdownTableBodyNone"><span class="tt">bub_pp%mu_g</span> </td><td class="markdownTableBodyNone">Gas dynamic viscosity </td><td class="markdownTableBodyNone">\(\mu_g / (\rho_0\,x_0\,u_0)\) </td></tr> <tr class="markdownTableRowEven"> <td class="markdownTableBodyNone"><span class="tt">bub_pp%vd</span> </td><td class="markdownTableBodyNone">Vapor diffusivity </td><td class="markdownTableBodyNone">\(D / (x_0\,u_0)\) </td></tr> <tr class="markdownTableRowOdd"> <td class="markdownTableBodyNone"><span class="tt">bub_pp%k_v</span> </td><td class="markdownTableBodyNone">Vapor thermal conductivity </td><td class="markdownTableBodyNone">\(k_v\,T_0 / (x_0\,\rho_0\,u_0^3)\) </td></tr> <tr class="markdownTableRowEven"> <td class="markdownTableBodyNone"><span class="tt">bub_pp%k_g</span> </td><td class="markdownTableBodyNone">Gas thermal conductivity </td><td class="markdownTableBodyNone">\(k_g\,T_0 / (x_0\,\rho_0\,u_0^3)\) </td></tr> <tr class="markdownTableRowOdd"> <td class="markdownTableBodyNone"><span class="tt">bub_pp%cp_v</span> </td><td class="markdownTableBodyNone">Vapor specific heat </td><td class="markdownTableBodyNone">\(c_{p,v}\,T_0 / u_0^2\) </td></tr> <tr class="markdownTableRowEven"> <td class="markdownTableBodyNone"><span class="tt">bub_pp%cp_g</span> </td><td class="markdownTableBodyNone">Gas specific heat </td><td class="markdownTableBodyNone">\(c_{p,g}\,T_0 / u_0^2\) </td></tr> <tr class="markdownTableRowOdd"> <td class="markdownTableBodyNone"><span class="tt">bub_pp%R_v</span> </td><td class="markdownTableBodyNone">Vapor gas constant </td><td class="markdownTableBodyNone">\(R_v\,T_0 / u_0^2\) </td></tr> <tr class="markdownTableRowEven"> <td class="markdownTableBodyNone"><span class="tt">bub_pp%R_g</span> </td><td class="markdownTableBodyNone">Gas gas constant </td><td class="markdownTableBodyNone">\(R_g\,T_0 / u_0^2\) </td></tr> <tr class="markdownTableRowOdd"> <td class="markdownTableBodyNone"><span class="tt">bub_pp%gam_v</span> </td><td class="markdownTableBodyNone">Vapor heat capacity ratio </td><td class="markdownTableBodyNone">Already dimensionless (no scaling) </td></tr> <tr class="markdownTableRowEven"> <td class="markdownTableBodyNone"><span class="tt">bub_pp%gam_g</span> </td><td class="markdownTableBodyNone">Gas heat capacity ratio </td><td class="markdownTableBodyNone">Already dimensionless (no scaling) </td></tr> <tr class="markdownTableRowOdd"> <td class="markdownTableBodyNone"><span class="tt">bub_pp%M_v</span> </td><td class="markdownTableBodyNone">Vapor molar mass </td><td class="markdownTableBodyNone">Consistent units; only ratios are used (no scaling needed) </td></tr> <tr class="markdownTableRowEven"> <td class="markdownTableBodyNone"><span class="tt">bub_pp%M_g</span> </td><td class="markdownTableBodyNone">Gas molar mass </td><td class="markdownTableBodyNone">Consistent units; only ratios are used (no scaling needed) </td></tr> </table> <p>When the reference scales match the bubble reference values (e.g., \(x_0 = R_{0,\text{ref}}\), \(p_0 = p_{0,\text{ref}}\), \(\rho_0 = \rho_{0,\text{ref}}\)), the reference parameters simplify to unity: <span class="tt">bub_pp%R0ref = 1</span>, <span class="tt">bub_pp%p0ref = 1</span>, <span class="tt">bub_pp%rho0ref = 1</span>.</p> <h4 class="doxsection"><a class="anchor" id="sec-flow-params-bubbles"></a> Flow Parameters with Bubbles</h4> <p>When bubbles are enabled, the flow-level parameters must also be non-dimensionalized with the same reference scales:</p> <table class="markdownTable"> <tr class="markdownTableHead"> <th class="markdownTableHeadNone">Parameter </th><th class="markdownTableHeadNone">Non-dimensional form </th></tr> <tr class="markdownTableRowOdd"> <td class="markdownTableBodyNone"><span class="tt">x_domain%beg</span>, <span class="tt">x_domain%end</span> </td><td class="markdownTableBodyNone">Domain bounds divided by \(x_0\) </td></tr> <tr class="markdownTableRowEven"> <td class="markdownTableBodyNone"><span class="tt">patch_icpp(i)%pres</span> </td><td class="markdownTableBodyNone">Pressure divided by \(p_0\) </td></tr> <tr class="markdownTableRowOdd"> <td class="markdownTableBodyNone"><span class="tt">patch_icpp(i)%alpha_rho(j)</span> </td><td class="markdownTableBodyNone">Partial density divided by \(\rho_0\) </td></tr> <tr class="markdownTableRowEven"> <td class="markdownTableBodyNone"><span class="tt">patch_icpp(i)%vel(j)</span> </td><td class="markdownTableBodyNone">Velocity divided by \(u_0\) </td></tr> <tr class="markdownTableRowOdd"> <td class="markdownTableBodyNone"><span class="tt">fluid_pp(i)%gamma</span> </td><td class="markdownTableBodyNone">\(1/(\gamma_i - 1)\) (dimensionless, same as without bubbles) </td></tr> <tr class="markdownTableRowEven"> <td class="markdownTableBodyNone"><span class="tt">fluid_pp(i)%pi_inf</span> </td><td class="markdownTableBodyNone">\(\gamma_i\,\pi_{\infty,i} / [(\gamma_i - 1)\,p_0]\) (scaled by reference pressure) </td></tr> <tr class="markdownTableRowOdd"> <td class="markdownTableBodyNone"><span class="tt">fluid_pp(i)%Re(1)</span> </td><td class="markdownTableBodyNone">\(\rho_0\,x_0\,u_0 / \mu_i\) (Reynolds number, inverse viscosity) </td></tr> <tr class="markdownTableRowEven"> <td class="markdownTableBodyNone"><span class="tt">dt</span> </td><td class="markdownTableBodyNone">Time step divided by \(t_0\) </td></tr> </table> <h4 class="doxsection"><a class="anchor" id="sec-two-viscosities"></a> Two Different Viscosity Parameters</h4> <p>MFC has two conceptually distinct viscosity-related parameters that serve different physical roles:</p> <ol type="1"> <li><b><span class="tt">fluid_pp(i)%Re(1)</span></b> — Used for the <b>macroscopic flow viscous stress tensor</b> (Navier-Stokes equations). This is \(1/\mu\) in dimensional simulations, or \(\rho_0 x_0 u_0 / \mu\) (a Reynolds number) when non-dimensionalized. It appears as a <b>divisor</b> in the viscous stress computation: <p class="formulaDsp"> \[\tau_{ij} \propto \frac{\nabla u}{\text{Re}}\] </p> Stored in the physical_parameters derived type (<span class="tt">src/common/m_derived_types.fpp</span>).</li> <li><b><span class="tt">bub_pp%mu_l</span></b> — Used for <b>microscale bubble wall viscous damping</b> (Rayleigh-Plesset / Keller-Miksis equations). This is the non-dimensional liquid viscosity \(\mu_l / (\rho_0 x_0 u_0)\). It appears as a <b>multiplier</b> in the bubble wall pressure: <p class="formulaDsp"> \[p_{bw} \ni -\frac{4\,\text{Re}_{\text{inv}}\,\dot{R}}{R}\] </p> Stored in the subgrid_bubble_physical_parameters derived type (<span class="tt">src/common/m_derived_types.fpp</span>).</li> </ol> <p>These two parameters represent viscous effects at fundamentally different scales — bulk flow dissipation vs. single-bubble-wall damping — and are stored in separate derived types with separate code paths. They are <b>not</b> interchangeable: <span class="tt">fluid_pp%Re(1)</span> is an inverse viscosity while <span class="tt">bub_pp%mu_l</span> is a viscosity (non-dimensionalized).</p> <h3 class="doxsection"><a class="anchor" id="sec-nondim-example"></a> Worked Examples</h3> <h4 class="doxsection"><a class="anchor" id="sec-bubble-example"></a> Example: Non-Dimensionalizing a Bubble Case</h4> <p>A typical bubble case setup in <span class="tt">case.py</span> follows this pattern:</p> <div class="fragment"><div class="line"><span class="keyword">import</span> math</div> <div class="line"> </div> <div class="line"><span class="comment"># Physical properties (SI units)</span></div> <div class="line">rho_l = 1.0e03 <span class="comment"># liquid density [kg/m³]</span></div> <div class="line">mu_l = 1.002e-03 <span class="comment"># liquid viscosity [kg/(m·s)]</span></div> <div class="line">ss = 0.07275 <span class="comment"># surface tension [kg/s²]</span></div> <div class="line">pv = 2.3388e03 <span class="comment"># vapor pressure [Pa]</span></div> <div class="line">gam_l = 7.15 <span class="comment"># liquid stiffened gas gamma</span></div> <div class="line">pi_inf = 306.0e06 <span class="comment"># liquid stiffened gas pi_inf [Pa]</span></div> <div class="line"> </div> <div class="line"><span class="comment"># Bubble reference values (SI)</span></div> <div class="line">R0ref = 10.0e-06 <span class="comment"># reference bubble radius [m]</span></div> <div class="line">p0ref = 112.9e03 <span class="comment"># reference bubble pressure [Pa]</span></div> <div class="line">rho0ref = rho_l <span class="comment"># reference density [kg/m³]</span></div> <div class="line"> </div> <div class="line"><span class="comment"># Derived reference scales</span></div> <div class="line">x0 = R0ref</div> <div class="line">p0 = p0ref</div> <div class="line">rho0 = rho0ref</div> <div class="line">u0 = math.sqrt(p0 / rho0)</div> <div class="line">t0 = x0 / u0</div> <div class="line"> </div> <div class="line"><span class="comment"># Non-dimensional inputs</span></div> <div class="line">params = {</div> <div class="line"> <span class="stringliteral">"bub_pp%R0ref"</span>: R0ref / x0, <span class="comment"># = 1.0</span></div> <div class="line"> <span class="stringliteral">"bub_pp%p0ref"</span>: p0ref / p0, <span class="comment"># = 1.0</span></div> <div class="line"> <span class="stringliteral">"bub_pp%rho0ref"</span>: rho0ref / rho0, <span class="comment"># = 1.0</span></div> <div class="line"> <span class="stringliteral">"bub_pp%ss"</span>: ss / (rho0 * x0 * u0**2), <span class="comment"># surface tension</span></div> <div class="line"> <span class="stringliteral">"bub_pp%pv"</span>: pv / p0, <span class="comment"># vapor pressure</span></div> <div class="line"> <span class="stringliteral">"bub_pp%mu_l"</span>: mu_l / (rho0 * x0 * u0), <span class="comment"># liquid viscosity</span></div> <div class="line"> </div> <div class="line"> <span class="stringliteral">"fluid_pp(1)%gamma"</span>: 1.0 / (gam_l - 1.0),</div> <div class="line"> <span class="stringliteral">"fluid_pp(1)%pi_inf"</span>: gam_l * (pi_inf / p0) / (gam_l - 1.0),</div> <div class="line"> <span class="stringliteral">"fluid_pp(1)%Re(1)"</span>: rho0 * x0 * u0 / mu_l, <span class="comment"># flow Re (inverse!)</span></div> <div class="line">}</div> </div><!-- fragment --><p>Note the inverse relationship: <span class="tt">fluid_pp%Re(1) = 1 / bub_pp%mu_l</span> when both use the same reference scales and the same physical viscosity. This is expected — they encode the same physical viscosity but in reciprocal forms for their respective equations.</p> <hr /> <h2 class="doxsection"><a class="anchor" id="autotoc_md149"></a> 2. Governing PDEs</h2> <h3 class="doxsection"><a class="anchor" id="autotoc_md150"></a> 2.1 Five-Equation Model (<span class="tt">model_eqns = 2</span>)</h3> <p>The primary workhorse model (Allaire et al. <a class="el" href="citelist.html#CITEREF_allaire02">[1]</a>; Wilfong et al. <a class="el" href="citelist.html#CITEREF_wilfong26">[68]</a> Sec. 2.1). The state vector is:</p> <p class="formulaDsp"> \[\mathbf{q} = \bigl(\alpha_1 \rho_1,\;\alpha_2 \rho_2,\;\ldots,\;\rho u_1,\;\rho u_2,\;\rho u_3,\;\rho E,\;\alpha_1,\;\alpha_2,\;\ldots\bigr)^T\] </p> <p><b>Continuity</b> (one per component):</p> <p class="formulaDsp"> \[\frac{\partial (\alpha_i \rho_i)}{\partial t} + \nabla \cdot (\alpha_i \rho_i\,\mathbf{u}) = 0\] </p> <p><b>Momentum:</b></p> <p class="formulaDsp"> \[\frac{\partial (\rho \mathbf{u})}{\partial t} + \nabla \cdot \bigl(\rho\,\mathbf{u} \otimes \mathbf{u} + p\,\mathbf{I} - \boldsymbol{\tau}^v\bigr) = 0\] </p> <p><b>Energy:</b></p> <p class="formulaDsp"> \[\frac{\partial (\rho E)}{\partial t} + \nabla \cdot \bigl[(\rho E + p)\,\mathbf{u} - \boldsymbol{\tau}^v \cdot \mathbf{u}\bigr] = 0\] </p> <p><b>Volume fraction advection:</b></p> <p class="formulaDsp"> \[\frac{\partial \alpha_i}{\partial t} + \mathbf{u} \cdot \nabla \alpha_i = K\,\nabla \cdot \mathbf{u}\] </p> <p>where the \(K\) term enforces interface conditions via the Wood sound speed:</p> <p class="formulaDsp"> \[K = \frac{\rho_2 c_2^2 - \rho_1 c_1^2}{\displaystyle\frac{\rho_1 c_1^2}{\alpha_1} + \displaystyle\frac{\rho_2 c_2^2}{\alpha_2}}\] </p> <p>Setting <span class="tt">alt_soundspeed = .true.</span> enables the \(K\) correction (Kapila et al. <a class="el" href="citelist.html#CITEREF_kapila01">[27]</a>, with Wood sound speed). Setting <span class="tt">alt_soundspeed = .false.</span> uses the Allaire variant without the \(K\) correction, which is conservative but does not strictly obey the second law of thermodynamics.</p> <p><b>Mixture rules:</b></p> <p class="formulaDsp"> \[1 = \sum_i \alpha_i, \qquad \rho = \sum_i \alpha_i \rho_i, \qquad \rho e = \sum_i \alpha_i \rho_i e_i\] </p> <h3 class="doxsection"><a class="anchor" id="autotoc_md151"></a> 2.2 Six-Equation Model (<span class="tt">model_eqns = 3</span>)</h3> <p>Allows pressure disequilibrium between phases (Saurel et al. <a class="el" href="citelist.html#CITEREF_saurel09">[53]</a>; Wilfong et al. <a class="el" href="citelist.html#CITEREF_wilfong26">[68]</a> Sec. 2.1).</p> <p><b>Continuity and momentum:</b> Same as the five-equation model.</p> <p><b>Separate phasic internal energy:</b></p> <p class="formulaDsp"> \[\frac{\partial (\alpha_i \rho_i e_i)}{\partial t} + \nabla \cdot (\alpha_i \rho_i e_i\,\mathbf{u}) + \alpha_i p_i\,\nabla \cdot \mathbf{u} = -\mu\,p_I\,(p_2 - p_1) - \alpha_i\,\boldsymbol{\tau}_i^v : \nabla \mathbf{u}\] </p> <p><b>Volume fraction:</b></p> <p class="formulaDsp"> \[\frac{\partial \alpha_1}{\partial t} + \mathbf{u} \cdot \nabla \alpha_1 = \mu\,(p_1 - p_2)\] </p> <p><b>Interfacial pressure:</b></p> <p class="formulaDsp"> \[p_I = \frac{z_2\,p_1 + z_1\,p_2}{z_1 + z_2}, \qquad z_i = \rho_i\,c_i\] </p> <p>Infinite pressure relaxation is applied at each Runge-Kutta stage to drive toward pressure equilibrium.</p> <p><b>Mixture speed of sound:</b></p> <p class="formulaDsp"> \[c^2 = \sum_k Y_k\,c_k^2\] </p> <p>With phase change (<span class="tt">relax = .true.</span>), additional source terms appear in the phasic energy and volume fraction equations:</p><ul> <li><b>Pressure relaxation:</b> \(\mu\,\delta p\) where \(\delta p = p_1 - p_2\)</li> <li><b>Thermal transfer:</b> \(Q = \theta\,(T_2 - T_1)\)</li> <li><b>Mass transfer:</b> \(\dot{m} = \nu\,(g_2 - g_1)\) (Gibbs free energy difference)</li> </ul> <p>See Section 8 (Phase Change) below for details.</p> <h3 class="doxsection"><a class="anchor" id="autotoc_md152"></a> 2.3 Other Model Variants</h3> <ul> <li><span class="tt">model_eqns = 1</span>: <b>Gamma/pi_inf model</b> — simplified single-fluid formulation using mixture \(\gamma\) and \(\pi_\infty\) directly without tracking individual volume fractions (Johnsen <a class="el" href="citelist.html#CITEREF_johnsen08">[26]</a>).</li> </ul> <hr /> <h2 class="doxsection"><a class="anchor" id="autotoc_md154"></a> 3. Equations of State</h2> <h3 class="doxsection"><a class="anchor" id="autotoc_md155"></a> 3.1 Stiffened Gas EOS (Menikoff and Plohr <a class="el" href="citelist.html#CITEREF_menikoff89">[40]</a>; Le Metayer et al. <a class="el" href="citelist.html#CITEREF_lemetayer04">[30]</a>; Wilfong et al. <a class="el" href="citelist.html#CITEREF_wilfong26">[68]</a> Sec. 2.2)</h3> <p>The primary closure for each phase:</p> <p class="formulaDsp"> \[p_k = (\gamma_k - 1)\,\rho_k\,e_k - \gamma_k\,\pi_{\infty,k}\] </p> <p>Equivalently:</p> <p class="formulaDsp"> \[e_k = \frac{p_k + \gamma_k\,\pi_{\infty,k}}{(\gamma_k - 1)\,\rho_k}\] </p> <p><b>Total energy relation:</b></p> <p class="formulaDsp"> \[\rho E = \Gamma\,p + \Pi_\infty + \frac{1}{2}\rho\,|\mathbf{u}|^2 + q_v\] </p> <p>where MFC internally tracks the transformed thermodynamic quantities:</p> <p class="formulaDsp"> \[\Gamma_k = \frac{1}{\gamma_k - 1}, \qquad \Pi_{\infty,k} = \frac{\gamma_k\,\pi_{\infty,k}}{\gamma_k - 1}\] </p> <p>and the mixture rules are arithmetic averages of these transformed quantities:</p> <p class="formulaDsp"> \[\Gamma = \sum_i \frac{\alpha_i}{\gamma_i - 1}, \qquad \Pi_\infty = \sum_i \frac{\alpha_i\,\gamma_i\,\pi_{\infty,i}}{\gamma_i - 1}, \qquad q_v = \sum_i \alpha_i\,\rho_i\,q_{v,i}\] </p> <p>The pressure is recovered from the total energy as:</p> <p class="formulaDsp"> \[p = \frac{\rho E - \frac{1}{2}\rho\,|\mathbf{u}|^2 - \Pi_\infty - q_v}{\Gamma}\] </p> <p><b>Phasic speed of sound:</b></p> <p class="formulaDsp"> \[c_k = \sqrt{\frac{\gamma_k\,(p + \pi_{\infty,k})}{\rho_k}}\] </p> <p><b>Wood mixture sound speed:</b></p> <p class="formulaDsp"> \[\frac{1}{\rho\,c^2} = \sum_k \frac{\alpha_k}{\rho_k\,c_k^2}\] </p> <p>Input parameters per fluid: <span class="tt">gamma</span> ( \(\Gamma_k = 1/(\gamma_k - 1)\)), <span class="tt">pi_inf</span> ( \(\Pi_{\infty,k} = \gamma_k\,\pi_{\infty,k}/(\gamma_k - 1)\)), <span class="tt">cv</span> ( \(c_{v,k}\)), <span class="tt">qv</span> ( \(q_{v,k}\)), <span class="tt">qvp</span> ( \(q'_{v,k}\)). Note that <span class="tt">gamma</span> and <span class="tt">pi_inf</span> are stored in transformed form, not as the raw physical values (see Section 1b).</p> <h3 class="doxsection"><a class="anchor" id="autotoc_md156"></a> 3.2 Ideal Gas EOS (Chemistry, <span class="tt">chemistry = .true.</span>)</h3> <p>For reacting gas mixtures:</p> <p class="formulaDsp"> \[p = \frac{\rho\,R_u\,T}{W}, \qquad W = \left(\sum_m \frac{Y_m}{W_m}\right)^{-1}\] </p> <p>Temperature is obtained from the internal energy by Newton iteration:</p> <p class="formulaDsp"> \[e_g - \sum_m e_m(T)\,Y_m = 0\] </p> <p><b>Species internal energy from enthalpy:</b></p> <p class="formulaDsp"> \[e_m(T) = \frac{\hat{h}_m(T) - R_u\,T}{W_m}\] </p> <p><b>NASA polynomial enthalpies</b> (McBride et al. <a class="el" href="citelist.html#CITEREF_mcbride93">[38]</a>):</p> <p class="formulaDsp"> \[\frac{\hat{h}_m}{R_u\,T} = \frac{C_0}{T} + \sum_{r=1}^{5} \frac{C_r\,T^{r-1}}{r}\] </p> <hr /> <h2 class="doxsection"><a class="anchor" id="autotoc_md158"></a> 4. Viscous Stress Tensor (<span class="tt">viscous = .true.</span>)</h2> <p><b>Newtonian viscous stress (no bulk viscosity by default):</b></p> <p class="formulaDsp"> \[\boldsymbol{\tau}^v = 2\,\eta\left(\mathbf{D} - \frac{1}{3}\,\text{tr}(\mathbf{D})\,\mathbf{I}\right)\] </p> <p>where the strain rate tensor is:</p> <p class="formulaDsp"> \[\mathbf{D} = \frac{1}{2}\bigl(\nabla \mathbf{u} + (\nabla \mathbf{u})^T\bigr)\] </p> <p><b>With bulk viscosity:</b></p> <p class="formulaDsp"> \[\tau_{ij} = \mu\left(\frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i}\right) + \left(\zeta - \frac{2\mu}{3}\right)\delta_{ij}\,\frac{\partial u_k}{\partial x_k}\] </p> <p><b>Cartesian components:</b></p> <p class="formulaDsp"> \[\tau_{xx} = \mu\left(2\,\frac{\partial u}{\partial x} - \frac{2}{3}\nabla\cdot\mathbf{u}\right), \qquad \tau_{xy} = \mu\left(\frac{\partial u}{\partial y} + \frac{\partial v}{\partial x}\right)\] </p> <p>and similarly for all other components. Cylindrical coordinate formulations include additional \(1/r\) terms.</p> <p><b>Viscosity averaging:</b></p> <p class="formulaDsp"> \[\frac{1}{\text{Re}_\text{mix}} = \sum_j \frac{\alpha_j}{\text{Re}_j}\] </p> <p>Input parameters: <span class="tt">Re_inv</span> (shear and volume Reynolds numbers per fluid).</p> <p><b>Non-Newtonian viscosity (<span class="tt">fluid_pp(i)%non_newtonian = .true.</span>):</b></p> <p>Per-fluid Herschel-Bulkley rheology with Papanastasiou regularization (C. <a class="el" href="citelist.html#CITEREF_papanastasiou87">[44]</a>) replaces the constant viscosity with a shear-rate-dependent effective viscosity, recomputed from the local strain-rate tensor at every time step:</p> <p class="formulaDsp"> \[\mu_{\rm eff}(\dot\gamma) = \frac{\tau_0}{\dot\gamma}\left(1 - e^{-m\,\dot\gamma}\right) + K\,\dot\gamma^{\,n-1}, \qquad \dot\gamma = \sqrt{2\,\mathbf{D}:\mathbf{D}}\] </p> <p>This covers power-law shear-thinning/thickening ( \(\tau_0 = 0\)), Bingham plastic ( \(n = 1,\ \tau_0 > 0\)), and general Herschel-Bulkley yield-stress fluids. The mixture rule above applies with \(1/\text{Re}_j = \mu_{{\rm eff},j}\) evaluated at the local shear rate; Newtonian and non-Newtonian fluids can be mixed. Supported with the HLL and HLLC Riemann solvers, including immersed boundaries. See <a class="el" href="case.html#sec-non-newtonian" title="19. Non-Newtonian (Herschel-Bulkley) Viscosity">the case documentation</a> for parameters, constraints, and validated example cases.</p> <hr /> <h2 class="doxsection"><a class="anchor" id="autotoc_md160"></a> 4a. Fourier Heat Conduction (<span class="tt">fluid_pp(i)%k_therm</span>)</h2> <p><b>Source:</b> <span class="tt">src/simulation/m_conduction.fpp</span></p> <p>Setting <span class="tt">fluid_pp(i)%k_therm > 0</span> on any fluid adds a Fourier conduction term to the energy equation:</p> <p class="formulaDsp"> \[\frac{\partial(\rho E)}{\partial t} + \nabla\cdot\bigl[(\rho E + p)\,\mathbf{u}\bigr] = \cdots + \nabla\cdot(k\,\nabla T)\] </p> <p>using a single thermal-equilibrium mixture temperature \(T\) (the same value carried in <span class="tt">q_T_sf</span> and written by <span class="tt">T_wrt</span>) and a volume-fraction-weighted mixture conductivity:</p> <p class="formulaDsp"> \[k = \sum_i \alpha_i\,k_i, \qquad T = \frac{(\Gamma+1)\,p + \Pi_\infty}{\sum_i \alpha_i\rho_i\,c_{v,i}\,n_i}\] </p> <p>The flux is direction-split and face-centered: a two-point difference of \(T\) across each face, exact for this term since it has no cross-derivatives (an axis-cell correction applies in cylindrical coordinates; see the limitation below). The closure is implemented for the stiffened-gas and ideal-gas equations of state only, and for <span class="tt">model_eqns = 2</span> (5-equation) or <span class="tt">model_eqns = 3</span> (6-equation) — both carry the volume fractions \(\alpha_i\) that weight \(k\), which <span class="tt">model_eqns = 1</span> (gamma law) does not. Heat conduction is independent of <span class="tt">viscous</span>: it can be active in an otherwise inviscid simulation. It is not supported together with <span class="tt">igr</span> or <span class="tt">chemistry</span>; see <a class="el" href="case.html#sec-fluid-materials" title="5. Fluid Material's">5. Fluid Material's</a> in the case documentation for the full set of input constraints.</p> <p>A thermal diffusion CFL limit (<span class="tt">TCFL</span>) is added to the adaptive time-step candidates alongside <span class="tt">ICFL</span>/<span class="tt">VCFL</span>/<span class="tt">CCFL</span>.</p> <p><b>Known limitation:</b> in cylindrical coordinates, the cell adjacent to the axis carries a non-converging -3% error in the conduction term. It is inherited from the two-point face-gradient every MFC diffusive flux uses, applied across the coordinate singularity at the axis — the same pattern as the two-point gradient in the chemistry diffusion flux (<span class="tt">src/common/m_chemistry.fpp</span>).</p> <hr /> <h2 class="doxsection"><a class="anchor" id="autotoc_md162"></a> 5. Cylindrical Coordinates (<span class="tt">cyl_coord = .true.</span>) (Wilfong et al. <a class="el" href="citelist.html#CITEREF_wilfong26">[68]</a> Sec. 2.3)</h2> <p>Additional geometric source terms appear with \(1/r\) factors in the continuity, momentum, and energy equations. Key modifications:</p> <ul> <li><b>Radial momentum:</b> extra \(p/r\) and \(\tau_{\theta\theta}/r\) terms</li> <li><p class="startli"><b>Viscous stress:</b> \(\tau_{yy}\) includes \(v/r\) corrections:</p> <p class="formulaDsp"> \[\tau_{yy} = \mu\left(\frac{4}{3}\frac{\partial v}{\partial r} - \frac{2}{3}\frac{\partial u}{\partial x} - \frac{2}{3}\frac{v}{r}\right)\] </p> </li> <li><b>Axis singularity:</b> axis placed at cell boundary with spectral filtering in the azimuthal direction</li> </ul> <hr /> <h2 class="doxsection"><a class="anchor" id="autotoc_md164"></a> 6. Sub-Grid Bubble Dynamics (Wilfong et al. <a class="el" href="citelist.html#CITEREF_wilfong26">[68]</a> Sec. 4.1)</h2> <h3 class="doxsection"><a class="anchor" id="autotoc_md165"></a> 6.1 Euler-Euler Bubbles (<span class="tt">bubbles_euler = .true.</span>)</h3> <p><b>Source:</b> <span class="tt">src/simulation/m_bubbles_EE.fpp</span>, <span class="tt">src/simulation/m_bubbles.fpp</span></p> <h4 class="doxsection"><a class="anchor" id="autotoc_md166"></a> 6.1.1 Method of Classes (Commander and Prosperetti <a class="el" href="citelist.html#CITEREF_commander89">[13]</a>; Ando et al. <a class="el" href="citelist.html#CITEREF_ando11">[2]</a>)</h4> <p><b>Modified mixture pressure:</b></p> <p class="formulaDsp"> \[p = (1 - \alpha)\,p_l + \alpha\left(\frac{R^3\,p_{bw}}{\bar{R}^3} + \frac{\rho\,R^3\,\dot{R}^2}{\bar{R}^3}\right)\] </p> <p><b>Modified stiffened gas for the liquid phase:</b></p> <p class="formulaDsp"> \[\Gamma_l\,p_l + \Pi_{\infty,l} = \frac{1}{1 - \alpha}\left(E - \frac{1}{2}\rho\,|\mathbf{u}|^2\right)\] </p> <p><b>Bubble wall pressure (polytropic):</b></p> <p class="formulaDsp"> \[p_{bw} = \left(p_0 + \frac{2\sigma}{R_0}\right)\left(\frac{R_0}{R}\right)^{3\gamma} - \frac{4\mu\,\dot{R}}{R} - \frac{2\sigma}{R}\] </p> <p><b>Void fraction transport:</b></p> <p class="formulaDsp"> \[\frac{\partial \alpha}{\partial t} + \mathbf{u} \cdot \nabla \alpha = \frac{3\,\alpha\,\bar{R}^2\,\dot{R}}{\bar{R}^3}\] </p> <p><b>Number density conservation:</b></p> <p class="formulaDsp"> \[\frac{\partial n_\text{bub}}{\partial t} + \nabla \cdot (n_\text{bub}\,\mathbf{u}) = 0\] </p> <p>where \(n = \frac{3}{4\pi}\,\frac{\alpha}{\bar{R}^3}\).</p> <p><b>Polydispersity</b> (<span class="tt">polydisperse = .true.</span>): Log-normal PDF discretized into \(N_\text{bin}\) equilibrium radii with standard deviation <span class="tt">poly_sigma</span>, integrated via Simpson's rule.</p> <h4 class="doxsection"><a class="anchor" id="autotoc_md167"></a> 6.1.2 Rayleigh-Plesset (<span class="tt">bubble_model = 3</span>) (Lord Rayleigh <a class="el" href="citelist.html#CITEREF_rayleigh17">[32]</a>; Plesset <a class="el" href="citelist.html#CITEREF_plesset49">[47]</a>)</h4> <p class="formulaDsp"> \[R\,\ddot{R} + \frac{3}{2}\,\dot{R}^2 = \frac{p_{bw} - p_\infty}{\rho_l}\] </p> <h4 class="doxsection"><a class="anchor" id="autotoc_md168"></a> 6.1.3 Keller-Miksis (<span class="tt">bubble_model = 2</span>) (Keller and Miksis <a class="el" href="citelist.html#CITEREF_keller80">[28]</a>)</h4> <p class="formulaDsp"> \[R\,\ddot{R}\left(1 - \frac{\dot{R}}{c}\right) + \frac{3}{2}\,\dot{R}^2\left(1 - \frac{\dot{R}}{3c}\right) = \frac{p_{bw} - p_\infty}{\rho_l}\left(1 + \frac{\dot{R}}{c}\right) + \frac{R\,\dot{p}_{bw}}{\rho_l\,c}\] </p> <h4 class="doxsection"><a class="anchor" id="autotoc_md169"></a> 6.1.4 Gilmore (<span class="tt">bubble_model = 1</span>) (Gilmore <a class="el" href="citelist.html#CITEREF_gilmore52">[19]</a>)</h4> <p>Enthalpy-based formulation with compressibility corrections via the Tait EOS:</p> <p class="formulaDsp"> \[R\,\ddot{R}\left(1 - \frac{\dot{R}}{C}\right) + \frac{3}{2}\,\dot{R}^2\left(1 - \frac{\dot{R}}{3C}\right) = H\left(1 + \frac{\dot{R}}{C}\right) + \frac{R\,\dot{H}}{C}\left(1 - \frac{\dot{R}}{C}\right)\] </p> <p>where the enthalpy difference is:</p> <p class="formulaDsp"> \[H = \frac{n_\text{tait}(1 + B)}{n_\text{tait} - 1}\left[\left(\frac{p_{bw}}{1+B} + 1\right)^{(n_\text{tait}-1)/n_\text{tait}} - \left(\frac{p_\infty}{1+B} + 1\right)^{(n_\text{tait}-1)/n_\text{tait}}\right]\] </p> <p>and the local liquid sound speed:</p> <p class="formulaDsp"> \[C = \sqrt{n_\text{tait}(1+B)\left(\frac{p_\infty}{1+B} + 1\right)^{(n_\text{tait}-1)/n_\text{tait}} + (n_\text{tait} - 1)\,H}\] </p> <h4 class="doxsection"><a class="anchor" id="autotoc_md170"></a> 6.1.5 Non-Polytropic Thermal Model (<span class="tt">polytropic = .false.</span>) (Preston et al. <a class="el" href="citelist.html#CITEREF_preston07">[48]</a>)</h4> <p><b>Internal bubble pressure ODE:</b></p> <p class="formulaDsp"> \[\dot{p}_b = \frac{3\gamma_b}{R}\left(-\dot{R}\,p_b + R_v\,T_{bw}\,\dot{m}_v + \frac{\gamma_b - 1}{\gamma_b}\,k_{bw}\left.\frac{\partial T}{\partial r}\right|_R\right)\] </p> <p><b>Vapor mass flux:</b></p> <p class="formulaDsp"> \[\dot{m}_v = \frac{D\,\rho_{bw}}{1 - \chi_{vw}}\left.\frac{\partial \chi_v}{\partial r}\right|_R\] </p> <h4 class="doxsection"><a class="anchor" id="autotoc_md171"></a> 6.1.6 QBMM Moment Transport (<span class="tt">qbmm = .true.</span>) (Bryngelson et al. <a class="el" href="citelist.html#CITEREF_bryngelson20">[7]</a>)</h4> <p><b>Population balance equation:</b></p> <p class="formulaDsp"> \[\frac{\partial f}{\partial t} + \frac{\partial (f\,\dot{R})}{\partial R} + \frac{\partial (f\,\ddot{R})}{\partial \dot{R}} = 0\] </p> <p><b>Moment transport:</b></p> <p class="formulaDsp"> \[\frac{\partial (n_\text{bub}\,\mu_i)}{\partial t} + \nabla \cdot (n_\text{bub}\,\mu_i\,\mathbf{u}) = n_\text{bub}\,\dot{\mu}_i\] </p> <p>where moments \(\mu_{i_1,i_2} = \int R^{i_1}\,\dot{R}^{i_2}\,f\,dR\,d\dot{R}\).</p> <p><b>CHyQMOM inversion</b> recovers 4 quadrature nodes \((w_j, R_j, \dot{R}_j)\) from 6 moments via:</p> <p class="formulaDsp"> \[\bar{u} = \frac{\mu_{10}}{\mu_{00}}, \quad \bar{v} = \frac{\mu_{01}}{\mu_{00}}, \quad c_{20} = \frac{\mu_{20}}{\mu_{00}} - \bar{u}^2, \quad c_{11} = \frac{\mu_{11}}{\mu_{00}} - \bar{u}\bar{v}, \quad c_{02} = \frac{\mu_{02}}{\mu_{00}} - \bar{v}^2\] </p> <h3 class="doxsection"><a class="anchor" id="autotoc_md172"></a> 6.2 Euler-Lagrange Bubbles (<span class="tt">bubbles_lagrange = .true.</span>) (Maeda and Colonius <a class="el" href="citelist.html#CITEREF_maeda18">[36]</a>)</h3> <p><b>Source:</b> <span class="tt">src/simulation/m_bubbles_EL.fpp</span></p> <p>Volume-averaged carrier flow equations with bubble source terms:</p> <p><b>Continuity:</b></p> <p class="formulaDsp"> \[\frac{\partial \rho_l}{\partial t} + \nabla \cdot (\rho_l\,\mathbf{u}_l) = \frac{\rho_l}{1 - \alpha}\left[\frac{\partial \alpha}{\partial t} + \mathbf{u}_l \cdot \nabla \alpha\right]\] </p> <p><b>Momentum:</b></p> <p class="formulaDsp"> \[\frac{\partial (\rho_l\,\mathbf{u}_l)}{\partial t} + \nabla \cdot (\rho_l\,\mathbf{u}_l \otimes \mathbf{u}_l + p\,\mathbf{I} - \boldsymbol{\tau}_l) = \frac{\rho_l\,\mathbf{u}_l}{1 - \alpha}\left[\frac{\partial \alpha}{\partial t} + \mathbf{u}_l \cdot \nabla \alpha\right] - \frac{\alpha}{1 - \alpha}\,\nabla \cdot (p\,\mathbf{I} - \boldsymbol{\tau}_l)\] </p> <p><b>Energy:</b></p> <p class="formulaDsp"> \[\frac{\partial E_l}{\partial t} + \nabla \cdot \bigl[(E_l + p)\,\mathbf{u}_l - \boldsymbol{\tau}_l \cdot \mathbf{u}_l\bigr] = \frac{E_l}{1 - \alpha}\left[\frac{\partial \alpha}{\partial t} + \mathbf{u}_l \cdot \nabla \alpha\right] - \frac{\alpha}{1 - \alpha}\,\nabla \cdot (p\,\mathbf{u}_l - \boldsymbol{\tau}_l \cdot \mathbf{u}_l)\] </p> <p>The left-hand side is the standard conservation law for the liquid phase; the right-hand side source terms capture the effect of the bubbles on the host liquid.</p> <p><b>Void fraction via regularization kernel:</b></p> <p class="formulaDsp"> \[\alpha(\mathbf{x}) = \sum_n V_n\,\delta_\sigma(\mathbf{x} - \mathbf{x}_n)\] </p> <p>where \(\delta_\sigma\) is a Gaussian kernel:</p> <p class="formulaDsp"> \[\delta_\sigma(\mathbf{r}) = \frac{1}{(2\pi\sigma^2)^{3/2}}\exp\!\left(-\frac{|\mathbf{r}|^2}{2\sigma^2}\right)\] </p> <p>with \(\sigma = \varepsilon_b \max(\Delta x^{1/3}_\text{cell},\;R_\text{bubble})\).</p> <p>Each bubble is tracked individually with Keller-Miksis dynamics and 4th-order adaptive Runge-Kutta time integration.</p> <p><b>Translational motion (<span class="tt">vel_model > 0</span>):</b></p> <p>Bubbles may also translate through the carrier flow. Let \(\mathbf{x}_b\) be the bubble position, \(\mathbf{u}_b\) its velocity, \(a\) its radius, and \(\mathbf{u}_l(\mathbf{x}_b)\) the carrier velocity interpolated to the bubble location (a Lagrange polynomial of order <span class="tt">fd_order</span>, which must be set when <span class="tt">vel_model > 0</span>).</p> <ul> <li><b>Tracer bubbles (<span class="tt">vel_model = 1</span>)</b> follow the local carrier velocity: <p class="formulaDsp"> \[\frac{d\mathbf{x}_b}{dt} = \mathbf{u}_l(\mathbf{x}_b).\] </p> </li> <li><b>Newton's second law (<span class="tt">vel_model = 2</span>)</b> integrates the bubble momentum: <p class="formulaDsp"> \[m_b\,\frac{d\mathbf{u}_b}{dt} = \mathbf{F}_D + \mathbf{F}_p + \mathbf{F}_g, \qquad \frac{d\mathbf{x}_b}{dt} = \mathbf{u}_b,\] </p> with bubble mass \(m_b\) and forces acting on the slip velocity \(\mathbf{u}_\text{rel} = \mathbf{u}_b - \mathbf{u}_l\):</li> </ul> <table class="markdownTable"> <tr class="markdownTableHead"> <th class="markdownTableHeadNone">Force </th><th class="markdownTableHeadNone"><span class="tt">case.py</span> control </th><th class="markdownTableHeadNone">Non-dimensional form </th></tr> <tr class="markdownTableRowOdd"> <td class="markdownTableBodyNone">Drag \(\mathbf{F}_D\) </td><td class="markdownTableBodyNone"><span class="tt">drag_model</span> </td><td class="markdownTableBodyNone">\(-\,c_D\,\pi\,a\,\mathbf{u}_\text{rel} / \text{Re}\), with \(c_D = 4\) free-slip (Hadamard <a class="el" href="citelist.html#CITEREF_hadamard1911">[22]</a>; Rybczynski <a class="el" href="citelist.html#CITEREF_rybczynski1911">[50]</a>), \(c_D = 6\) no-slip Stokes (Stokes <a class="el" href="citelist.html#CITEREF_stokes1851">[57]</a>), \(c_D = 12\) Levich (Levich <a class="el" href="citelist.html#CITEREF_levich1962">[31]</a>) </td></tr> <tr class="markdownTableRowEven"> <td class="markdownTableBodyNone">Pressure \(\mathbf{F}_p\) </td><td class="markdownTableBodyNone"><span class="tt">pressure_force</span> </td><td class="markdownTableBodyNone">\(-V_b\,\nabla p\), with bubble volume \(V_b = \frac{4}{3}\pi a^3\) </td></tr> <tr class="markdownTableRowOdd"> <td class="markdownTableBodyNone">Gravity \(\mathbf{F}_g\) </td><td class="markdownTableBodyNone"><span class="tt">gravity_force</span> </td><td class="markdownTableBodyNone">\(m_b\,\mathbf{g}\), with \(\mathbf{g}\) the body-force acceleration </td></tr> </table> <p>Here \(\text{Re}\) is the mixture Reynolds number — the same <span class="tt">fluid_pp%Re(1)</span> that scales the viscous stress tensor (<a class="el" href="#sec-two-viscosities" title="Two Different Viscosity Parameters">Two Different Viscosity Parameters</a>) — so the drag scales with the liquid viscosity. The three drag models increase in magnitude free-slip \(<\) no-slip \(<\) Levich; see Magnaudet and Eames <a class="el" href="citelist.html#CITEREF_magnaudet2000">[37]</a> for a review of these bubble-drag regimes.</p> <hr /> <h2 class="doxsection"><a class="anchor" id="autotoc_md174"></a> 7. Fluid-Structure Interaction</h2> <h3 class="doxsection"><a class="anchor" id="autotoc_md175"></a> 7.1 Hypoelastic Model (<span class="tt">hypoelasticity = .true.</span>) (Rodriguez and Johnsen <a class="el" href="citelist.html#CITEREF_rodriguez19">[49]</a>; Wilfong et al. <a class="el" href="citelist.html#CITEREF_wilfong26">[68]</a> Sec. 4.1.6)</h3> <p><b>Source:</b> <span class="tt">src/simulation/m_hypoelastic.fpp</span></p> <p><b>Cauchy stress decomposition:</b></p> <p class="formulaDsp"> \[\sigma_{ij} = -p\,\delta_{ij} + \tau_{ij}^{(v)} + \tau_{ij}^{(e)}\] </p> <p><b>Elastic energy contribution to total energy:</b></p> <p class="formulaDsp"> \[E = e + \frac{|\mathbf{u}|^2}{2} + \frac{\boldsymbol{\tau}^e : \boldsymbol{\tau}^e}{4\,\rho\,G}\] </p> <p><b>Elastic stress evolution:</b></p> <p class="formulaDsp"> \[\frac{\partial (\rho\,\boldsymbol{\tau}^e)}{\partial t} + \nabla \cdot (\rho\,\boldsymbol{\tau}^e \otimes \mathbf{u}) = \mathbf{S}^e\] </p> <p><b>Source term:</b></p> <p class="formulaDsp"> \[\mathbf{S}^e = \rho\bigl(\mathbf{l} \cdot \boldsymbol{\tau}^e + \boldsymbol{\tau}^e \cdot \mathbf{l}^T - \boldsymbol{\tau}^e\,\text{tr}(\mathbf{D}) + 2G\,\mathbf{D}^d\bigr)\] </p> <p>where \(\mathbf{l} = \nabla \mathbf{u}\) is the velocity gradient and \(\mathbf{D}^d = \mathbf{D} - \frac{1}{3}\text{tr}(\mathbf{D})\,\mathbf{I}\) is the deviatoric strain rate.</p> <p><b>Lie objective temporal derivative (Kelvin-Voigt):</b></p> <p class="formulaDsp"> \[\hat{\boldsymbol{\tau}}^e = \frac{D\boldsymbol{\tau}^e}{Dt} - \mathbf{l} \cdot \boldsymbol{\tau}^e - \boldsymbol{\tau}^e \cdot \mathbf{l}^T + \boldsymbol{\tau}^e\,\text{tr}(\mathbf{D}) = 2G\,\mathbf{D}^d\] </p> <p>This adds 6 additional transport equations in 3D (symmetric stress tensor: \(\tau_{xx}^e, \tau_{xy}^e, \tau_{yy}^e, \tau_{xz}^e, \tau_{yz}^e, \tau_{zz}^e\)).</p> <h3 class="doxsection"><a class="anchor" id="autotoc_md176"></a> 7.2 Continuum Damage (<span class="tt">cont_damage = .true.</span>) (Cao et al. <a class="el" href="citelist.html#CITEREF_cao19">[10]</a>; Spratt <a class="el" href="citelist.html#CITEREF_spratt24">[56]</a> Sec. 4.1.2)</h3> <p><b>Source:</b> <span class="tt">src/simulation/m_hypoelastic.fpp</span></p> <p>A scalar damage field \(D \in [0,1]\) is transported with the damageable-solid partial mass \(m_s = \sum_{i:\,G_i > 0} \alpha_i \rho_i\):</p> <p class="formulaDsp"> \[\frac{\partial (m_s D)}{\partial t} + \nabla \cdot (m_s D\, \mathbf{u}) = m_s\,\dot{D}\] </p> <p>Damage grows when the maximum principal Cauchy stress \(\sigma_1 = \lambda_{\max}(-p\mathbf{I} + \boldsymbol{\tau}^e)\) exceeds \(\tau^*\):</p> <p class="formulaDsp"> \[\dot{D} = \bigl(\bar{\alpha}\,\max(\sigma_1 - \tau^*,\, 0)\bigr)^{s}\] </p> <p>The damaged shear modulus is</p> <p class="formulaDsp"> \[G = G_0(1-D),\] </p> <p>and is used in the elastic stress evolution and HLL/HLLC wave speeds; elastic energy uses the undamaged modulus \(G_0\).</p> <h2 class="doxsection"><a class="anchor" id="autotoc_md177"></a> 8. Phase Change (<span class="tt">relax = .true.</span>) (Wilfong et al. <a class="el" href="citelist.html#CITEREF_wilfong26">[68]</a> Sec. 4.1.3)</h2> <p><b>Source:</b> <span class="tt">src/common/m_phase_change.fpp</span></p> <h3 class="doxsection"><a class="anchor" id="autotoc_md178"></a> 8.1 pT-Relaxation (<span class="tt">relax_model = 5</span>) (Saurel et al. <a class="el" href="citelist.html#CITEREF_saurel08">[52]</a>)</h3> <p>\(N\)-fluid pressure-temperature equilibrium. The equilibrium condition is:</p> <p class="formulaDsp"> \[f(p) = \sum_i \alpha_i - 1 = 0\] </p> <p><b>Temperature from energy conservation:</b></p> <p class="formulaDsp"> \[T = \frac{\rho e + p - \sum_i (\alpha_i \rho_i)\,q_{v,i}}{\sum_i (\alpha_i \rho_i)\,c_{v,i}\,\gamma_i}\] </p> <p><b>Newton residual:</b></p> <p class="formulaDsp"> \[g(p) = \sum_i \frac{(\gamma_i - 1)\,(\alpha_i \rho_i)\,c_{v,i}}{(p + \pi_{\infty,i})} \cdot \frac{\rho e + p - \sum_j (\alpha_j \rho_j)\,q_{v,j}}{\sum_j (\alpha_j \rho_j)\,c_{v,j}\,\gamma_j}\] </p> <p>Solved via Newton's method for the equilibrium pressure.</p> <h3 class="doxsection"><a class="anchor" id="autotoc_md179"></a> 8.2 pTg-Relaxation (<span class="tt">relax_model = 6</span>) (Zein et al. <a class="el" href="citelist.html#CITEREF_zein10">[69]</a>)</h3> <p>Two coupled equations for \((\alpha_1 \rho_1,\;p)\):</p> <p><b>Gibbs free energy equilibrium (Clausius-Clapeyron):</b></p> <p class="formulaDsp"> \[F_1 = T\left[(c_{v,l}\gamma_l - c_{v,v}\gamma_v)(1 - \ln T) - (q'_l - q'_v) + c_{v,l}(\gamma_l - 1)\ln(p + \pi_{\infty,l}) - c_{v,v}(\gamma_v - 1)\ln(p + \pi_{\infty,v})\right] + q_{v,l} - q_{v,v} = 0\] </p> <p><b>Energy conservation constraint:</b></p> <p class="formulaDsp"> \[F_2 = \rho e + p + m_l\,(q_{v,v} - q_{v,l}) - m_T\,q_{v,v} - m_{qD} + \frac{m_l\,(c_{v,v}\,\gamma_v - c_{v,l}\,\gamma_l) - m_T\,c_{v,v}\,\gamma_v - m_{cpD}}{m_l\left(\frac{c_{v,l}\,(\gamma_l - 1)}{p + \pi_{\infty,l}} - \frac{c_{v,v}\,(\gamma_v - 1)}{p + \pi_{\infty,v}}\right) + \frac{m_T\,c_{v,v}\,(\gamma_v - 1)}{p + \pi_{\infty,v}} + m_{cvgp}} = 0\] </p> <p>where \(m_T\) is the total mass, \(m_l = \alpha_l \rho_l\) is the liquid partial density, and \(m_{qD}\), \(m_{cpD}\), \(m_{cvgp}\) are auxiliary thermodynamic sums over additional fluids (beyond the phase-changing pair).</p> <p>Solved via 2D Newton-Raphson.</p> <hr /> <h2 class="doxsection"><a class="anchor" id="autotoc_md181"></a> 9. Chemistry and Combustion (<span class="tt">chemistry = .true.</span>) (Wilfong et al. <a class="el" href="citelist.html#CITEREF_wilfong26">[68]</a> Sec. 4.1.7)</h2> <p><b>Source:</b> <span class="tt">src/common/m_chemistry.fpp</span></p> <p><b>Species transport:</b></p> <p class="formulaDsp"> \[\frac{\partial (\rho_g\,Y_m)}{\partial t} + \frac{\partial (\rho_g\,u_i\,Y_m)}{\partial x_i} = W_m\,\dot{\omega}_m\] </p> <p><b>Net production rate:</b></p> <p class="formulaDsp"> \[\dot{\omega}_m = \sum_n (\nu''_{mn} - \nu'_{mn})\,\mathcal{R}_n\] </p> <p><b>Reaction rate (law of mass action):</b></p> <p class="formulaDsp"> \[\mathcal{R}_n = k_n(T)\left[\prod_j \left(\frac{\rho_g\,Y_j}{W_j}\right)^{\nu'_{jn}} - \frac{1}{K_n}\prod_k \left(\frac{\rho_g\,Y_k}{W_k}\right)^{\nu''_{kn}}\right]\] </p> <p><b>Arrhenius rate:</b></p> <p class="formulaDsp"> \[k_n(T) = A_n\,T^{b_n}\exp\!\left(-\frac{T_{a,n}}{T}\right)\] </p> <p><b>Molecular diffusion</b> (<span class="tt">transport_model</span>):</p><ul> <li><b>Mixture-average:</b> Species-specific diffusion coefficients \(D_m^\text{mix}\), mass flux: \(\dot{m}_k = \rho\,D_k^\text{mix}\,(W_k / W_\text{mix})\,\partial X_k / \partial x\)</li> <li><b>Unity Lewis number:</b> \(D_m = \lambda / (\rho\,c_p)\)</li> </ul> <p>Enthalpy flux with diffusion:</p> <p class="formulaDsp"> \[q_\text{diff} = \lambda\,\frac{\partial T}{\partial x} + \sum_k h_k\,\dot{m}_k\] </p> <p>Reaction mechanisms are compiled into Fortran by MFC's own thermochemistry generator, using Cantera to load mechanism and transport data. The generator derives from the Fortran path of Pyrometheus (Cisneros-Garibay et al. <a class="el" href="citelist.html#CITEREF_cisneros26">[12]</a>), with its MIT attribution retained. The generated routines support CPU, OpenACC, and OpenMP execution; MFC does not require Pyrometheus, JAX, or automatic differentiation. See <a class="el" href="thermochemistry.html" title="Thermochemistry implementation">Thermochemistry implementation</a>.</p> <hr /> <h2 class="doxsection"><a class="anchor" id="autotoc_md183"></a> 10. Surface Tension (<span class="tt">surface_tension = .true.</span>) (Schmidmayer et al. <a class="el" href="citelist.html#CITEREF_schmidmayer17">[54]</a>; Wilfong et al. <a class="el" href="citelist.html#CITEREF_wilfong26">[68]</a> Sec. 4.1.8)</h2> <p><b>Source:</b> <span class="tt">src/simulation/m_surface_tension.fpp</span></p> <p><b>Color function advection:</b></p> <p class="formulaDsp"> \[\frac{\partial c}{\partial t} + \mathbf{u} \cdot \nabla c = 0\] </p> <p><b>Capillary stress tensor (CSF model):</b></p> <p class="formulaDsp"> \[\boldsymbol{\Omega} = -\sigma\left(\|\nabla c\|\,\mathbf{I} - \frac{\nabla c \otimes \nabla c}{\|\nabla c\|}\right)\] </p> <p>In component form, with \(\hat{w}_i = (\partial c / \partial x_i) / \|\nabla c\|\):</p> <p class="formulaDsp"> \[\Omega_{xx} = -\sigma\,(\hat{w}_y^2 + \hat{w}_z^2)\,\|\nabla c\|, \qquad \Omega_{xy} = \sigma\,\hat{w}_x\,\hat{w}_y\,\|\nabla c\|\] </p> <p>The capillary stress divergence is added to the momentum and energy equations. The total energy equation becomes:</p> <p class="formulaDsp"> \[\frac{\partial (\rho E + \varepsilon_0)}{\partial t} + \nabla \cdot \bigl[(\rho E + \varepsilon_0 + p)\,\mathbf{u} + (\boldsymbol{\Omega} - \boldsymbol{\tau}^v) \cdot \mathbf{u}\bigr] = 0\] </p> <p><b>Capillary mixture energy:</b></p> <p class="formulaDsp"> \[\varepsilon_0 = \sigma\,\|\nabla c\|\] </p> <hr /> <h2 class="doxsection"><a class="anchor" id="autotoc_md185"></a> 11. Magnetohydrodynamics</h2> <h3 class="doxsection"><a class="anchor" id="autotoc_md186"></a> 11.1 Ideal MHD (<span class="tt">mhd = .true.</span>) (Wilfong et al. <a class="el" href="citelist.html#CITEREF_wilfong26">[68]</a> Sec. 4.1.9)</h3> <p><b>Continuity:</b></p> <p class="formulaDsp"> \[\frac{\partial \rho}{\partial t} + \nabla \cdot (\rho\,\mathbf{u}) = 0\] </p> <p><b>Momentum:</b></p> <p class="formulaDsp"> \[\frac{\partial (\rho\,\mathbf{u})}{\partial t} + \nabla \cdot \left[\rho\,\mathbf{u} \otimes \mathbf{u} + \left(p + \frac{|\mathbf{B}|^2}{2}\right)\mathbf{I} - \mathbf{B} \otimes \mathbf{B}\right] = 0\] </p> <p><b>Energy:</b></p> <p class="formulaDsp"> \[\frac{\partial \mathcal{E}}{\partial t} + \nabla \cdot \left[\left(\mathcal{E} + p + \frac{|\mathbf{B}|^2}{2}\right)\mathbf{u} - (\mathbf{u} \cdot \mathbf{B})\,\mathbf{B}\right] = 0\] </p> <p><b>Induction:</b></p> <p class="formulaDsp"> \[\frac{\partial \mathbf{B}}{\partial t} + \nabla \cdot (\mathbf{u} \otimes \mathbf{B} - \mathbf{B} \otimes \mathbf{u}) = 0\] </p> <p><b>Total energy:</b></p> <p class="formulaDsp"> \[\mathcal{E} = \rho\,e + \frac{1}{2}\rho\,|\mathbf{u}|^2 + \frac{|\mathbf{B}|^2}{2}\] </p> <p><b>Fast magnetosonic speed:</b></p> <p class="formulaDsp"> \[c_f = \sqrt{\frac{1}{2}\left(c_s^2 + v_A^2 + \sqrt{(c_s^2 + v_A^2)^2 - 4\,c_s^2\,v_A^2\cos^2\theta}\right)}\] </p> <p><b>Alfven speed:</b></p> <p class="formulaDsp"> \[v_A = \sqrt{\frac{|\mathbf{B}|^2}{\rho}}\] </p> <p>Uses the HLLD Riemann solver (<span class="tt">riemann_solver = 4</span>). Hyperbolic divergence cleaning (<span class="tt">hyper_cleaning = .true.</span>) via the GLM method (Dedner et al. <a class="el" href="citelist.html#CITEREF_dedner02">[16]</a>).</p> <h3 class="doxsection"><a class="anchor" id="autotoc_md187"></a> 11.2 Relativistic MHD (<span class="tt">relativity = .true.</span>) (Wilfong et al. <a class="el" href="citelist.html#CITEREF_wilfong26">[68]</a> Sec. 4.1.10)</h3> <p><b>Conserved variables:</b></p> <p class="formulaDsp"> \[\mathbf{U} = (D,\;\mathbf{m},\;\tau,\;\mathbf{B})^T\] </p> <p>where:</p> <p class="formulaDsp"> \[D = \Gamma\,\rho, \qquad \mathbf{m} = \Gamma^2\rho h\,\mathbf{u} + |\mathbf{B}|^2\mathbf{u} - (\mathbf{u} \cdot \mathbf{B})\,\mathbf{B}\] </p> <p class="formulaDsp"> \[\tau = \Gamma^2\rho h - p + \frac{|\mathbf{B}|^2}{2} + \frac{|\mathbf{u}|^2|\mathbf{B}|^2 - (\mathbf{B} \cdot \mathbf{u})^2}{2} - \Gamma\,\rho\] </p> <p>Primitive recovery uses Newton-Raphson on the nonlinear conserved-to-primitive relation.</p> <hr /> <h2 class="doxsection"><a class="anchor" id="autotoc_md189"></a> 12. Information Geometric Regularization (<span class="tt">igr = .true.</span>) (Wilfong et al. <a class="el" href="citelist.html#CITEREF_wilfong25a">[67]</a>)</h2> <p><b>Source:</b> <span class="tt">src/simulation/m_igr.fpp</span></p> <p><b>Modified momentum with entropic pressure</b> \(\Sigma\)**:**</p> <p class="formulaDsp"> \[\frac{\partial (\rho\,\mathbf{u})}{\partial t} + \nabla \cdot \bigl[\rho\,\mathbf{u} \otimes \mathbf{u} + (p + \Sigma)\,\mathbf{I} - \boldsymbol{\tau}\bigr] = 0\] </p> <p><b>Elliptic PDE for</b> \(\Sigma\)**:**</p> <p class="formulaDsp"> \[\alpha\left[\text{tr}(\nabla \mathbf{u})^2 + \text{tr}^2(\nabla \mathbf{u})\right] = \frac{\Sigma}{\rho} - \alpha\,\nabla \cdot \left(\frac{\nabla \Sigma}{\rho}\right)\] </p> <p>where \(\alpha \sim \Delta x^2\) (regularization strength proportional to mesh spacing squared):</p> <p class="formulaDsp"> \[\alpha_\text{IGR} = \alpha_\text{factor} \cdot \max(\Delta x,\;\Delta y,\;\Delta z)^2\] </p> <p><b>RHS strain-rate source (3D):</b></p> <p class="formulaDsp"> \[\text{RHS} = \alpha\left[2\left(\frac{\partial u}{\partial y}\frac{\partial v}{\partial x} + \frac{\partial u}{\partial z}\frac{\partial w}{\partial x} + \frac{\partial v}{\partial z}\frac{\partial w}{\partial y}\right) + \left(\frac{\partial u}{\partial x}\right)^2 + \left(\frac{\partial v}{\partial y}\right)^2 + \left(\frac{\partial w}{\partial z}\right)^2 + (\nabla \cdot \mathbf{u})^2\right]\] </p> <p><b>Iterative solver:</b> Jacobi (<span class="tt">igr_iter_solver = 1</span>) or Gauss-Seidel (<span class="tt">igr_iter_solver = 2</span>), up to <span class="tt">num_igr_iters</span> iterations (default 5).</p> <p>Uses Lax-Friedrichs flux (replaces WENO + Riemann solver).</p> <hr /> <h2 class="doxsection"><a class="anchor" id="autotoc_md191"></a> 13. Body Forces (<span class="tt">bf_x</span>, <span class="tt">bf_y</span>, <span class="tt">bf_z</span>)</h2> <p><b>Source:</b> <span class="tt">src/simulation/m_body_forces.fpp</span></p> <p><b>Time-dependent acceleration:</b></p> <p class="formulaDsp"> \[a_i(t) = g_i + k_i\sin(\omega_i\,t - \phi_i)\] </p> <p><b>Momentum source:</b></p> <p class="formulaDsp"> \[\frac{\partial (\rho\,u_i)}{\partial t}\bigg|_\text{bf} = \rho\,a_i(t)\] </p> <p><b>Energy source:</b></p> <p class="formulaDsp"> \[\frac{\partial E}{\partial t}\bigg|_\text{bf} = \rho\,\mathbf{u} \cdot \mathbf{a}(t)\] </p> <p><b>Mean pressure gradient</b> (<span class="tt">mean_dpdx</span>, <span class="tt">mean_mflux</span>), \(f_i = -\partial\bar p/\partial x_i\):</p> <p class="formulaDsp"> \[\frac{\partial (\rho\,u_i)}{\partial t}\bigg|_\text{mpg} = f_i, \qquad \frac{\partial E}{\partial t}\bigg|_\text{mpg} = \mathbf{u}\cdot\mathbf{f}\] </p> <p>With a target bulk mass flux \(Q_i^\ast\), \(f_i^{n+1} = f_i^n + (G/\Delta t)(Q_i^\ast - 2Q_i^n + Q_i^{n-1})\).</p> <p><b>Mean gas temperature source</b> (<span class="tt">mean_qdot</span>, <span class="tt">mean_temp</span>):</p> <p class="formulaDsp"> \[\frac{\partial E}{\partial t}\bigg|_\text{mts} = \dot q, \qquad \dot q^{n+1} = \dot q^n + \frac{G\,\overline{\rho c_v}}{\Delta t}\left(T_b^\ast - 2T_b^n + T_b^{n-1}\right),\] </p> <p>with \(T_b\) the mass-averaged fluid temperature.</p> <hr /> <h2 class="doxsection"><a class="anchor" id="autotoc_md193"></a> 14. Acoustic Sources (<span class="tt">acoustic_source = .true.</span>)</h2> <p><b>Source:</b> <span class="tt">src/simulation/m_acoustic_src.fpp</span></p> <p>Source terms added to the RHS of the governing equations.</p> <p>Formulation follows Maeda and Colonius <a class="el" href="citelist.html#CITEREF_maeda17">[35]</a>.</p> <p><b>Discrete delta function (spatial support):</b></p> <p class="formulaDsp"> \[\delta_h(r) = \frac{1}{(2\pi\sigma^2)^{d/2}}\exp\!\left(-\frac{r^2}{2\sigma^2}\right)\] </p> <p><b>Forcing form (added to conservative variables):</b></p> <p class="formulaDsp"> \[\mathbf{s}_\text{ac} = \Omega_\Gamma\,f(t)\left[\frac{1}{c_0},\;\cos\theta,\;\sin\theta,\;\frac{c_0^2}{\gamma - 1}\right]^T\] </p> <p><b>Temporal profiles:</b></p><ul> <li><b>Sine</b> (<span class="tt">pulse = 1</span>): \(S(t) = M\sin(\omega(t - t_\text{delay}))\)</li> <li><b>Gaussian</b> (<span class="tt">pulse = 2</span>): \(S(t) = M\exp\!\bigl(-\frac{1}{2}((t - t_\text{delay})/\sigma_t)^2\bigr)\)</li> <li><b>Square</b> (<span class="tt">pulse = 3</span>): \(S(t) = M\,\text{sign}(\sin(\omega(t - t_\text{delay})))\)</li> <li><b>Broadband</b> (<span class="tt">pulse = 4</span>): superposition of multiple frequencies across a bandwidth</li> </ul> <p><b>Spatial supports:</b> planar, spherical transducer, cylindrical transducer, transducer array (arcuate, annular, circular).</p> <hr /> <h2 class="doxsection"><a class="anchor" id="autotoc_md195"></a> 15. Numerical Methods</h2> <h3 class="doxsection"><a class="anchor" id="autotoc_md196"></a> 15.1 Spatial Reconstruction</h3> <h4 class="doxsection"><a class="anchor" id="autotoc_md197"></a> WENO (<span class="tt">weno_order = 3, 5, 7</span>)</h4> <p><b>Source:</b> <span class="tt">src/simulation/m_weno.fpp</span></p> <p>Weighted sum of candidate polynomials at cell interfaces:</p> <p class="formulaDsp"> \[f_{i+1/2} = \sum_r \omega_r\,f_{i+1/2}^{(r)}\] </p> <p><b>WENO-JS</b> (Jiang and Shu <a class="el" href="citelist.html#CITEREF_jiang96">[25]</a>, default):</p> <p class="formulaDsp"> \[\alpha_r = \frac{d_r}{(\beta_r + \varepsilon)^2}, \qquad \omega_r = \frac{\alpha_r}{\sum_s \alpha_s}\] </p> <p>where \(d_r\) are ideal weights, \(\beta_r\) are smoothness indicators, and \(\varepsilon\) is a small regularization parameter (<span class="tt">weno_eps</span>).</p> <p><b>WENO-M</b> (<span class="tt">mapped_weno = .true.</span>): Henrick et al. <a class="el" href="citelist.html#CITEREF_henrick05">[24]</a> mapped weights for improved accuracy at critical points:</p> <p class="formulaDsp"> \[\omega_M^{(r)} = \frac{d_r\bigl(1 + d_r - 3\omega_0^{(r)} + (\omega_0^{(r)})^2\bigr)\,\omega_0^{(r)}}{d_r^2 + \omega_0^{(r)}(1 - 2d_r)}, \qquad \omega^{(r)} = \frac{\omega_M^{(r)}}{\sum_s \omega_M^{(s)}}\] </p> <p><b>WENO-Z</b> (<span class="tt">wenoz = .true.</span>): Borges et al. <a class="el" href="citelist.html#CITEREF_borges08">[6]</a> improved weights with global smoothness measure:</p> <p class="formulaDsp"> \[\alpha_r = d_r\left(1 + \left(\frac{\tau}{\beta_r + \varepsilon}\right)^q\right), \qquad \tau = |\beta_0 - \beta_{k-1}|\] </p> <p>The parameter \(q\) controls the convergence rate at critical points (typically \(q = 1\) for fifth-order reconstruction, as used in MFC).</p> <p><b>TENO</b> (<span class="tt">teno = .true.</span>): Fu et al. <a class="el" href="citelist.html#CITEREF_fu16">[17]</a> targeted ENO with smoothness threshold \(C_T\) (<span class="tt">teno_CT</span>):</p> <p class="formulaDsp"> \[\gamma_r = 1 + \frac{\tau}{\beta_r}, \qquad \xi_r = \frac{\gamma_r}{\sum_s \gamma_s}\] </p> <p>If \(\xi_r < C_T\), set \(\alpha_r = 0\) (stencil excluded).</p> <p>Primitive variable reconstruction is used to avoid spurious oscillations at interfaces.</p> <h4 class="doxsection"><a class="anchor" id="autotoc_md198"></a> MUSCL (<span class="tt">muscl_order = 2</span>)</h4> <p><b>Source:</b> <span class="tt">src/simulation/m_muscl.fpp</span></p> <p class="formulaDsp"> \[q_L(j) = q(j) - \frac{1}{2}\,\phi(r)\,\Delta q, \qquad q_R(j) = q(j) + \frac{1}{2}\,\phi(r)\,\Delta q\] </p> <p><b>Five slope limiters</b> (<span class="tt">muscl_lim</span>):</p><ol type="1"> <li><b>Minmod:</b> \(\phi = \text{sign}(a)\,\min(|a|,\;|b|)\) if \(ab > 0\), else \(0\)</li> |
| You are visitor 01402520 — GeezerVison does not store the pages you look up. Best viewed with Netscape Navigator 4.7 or Internet Explorer 6 at 800×600. No cookies required. |