Add wall-only part of Y3 phenolics model #875

tulioricci · 2023-04-13T22:23:31Z

Material response added to "multiphysics":

Y3 model: First PR to include phenolics materials. Adding thermophysical properties and governing equations to evaluate composite material response. It doesn't include the coupling with the fluid part (for now).
Y2 model: Include oxidation model (previously, all in the drivers) in specific files. Tentatively keeping the same structure for solving different but related physics (i.e., Y2 and Y3). In time, the two models will be combined (i.e., both phenolic resin and fibers must be solved simultaneously).

The EOS of gas and wall are included in gas_model.py and wall_model.py. These have been exercised in the McKenna burner driver (which is similar to the prediction driver): https://github.com/illinois-ceesd/drivers_mckenna.

The material properties are added in specific files in the new folder "materials". In time, HARLEM will be added as the prediction-targeted material.

Add new boundary condition to diffusion.py: impose the actual flux into the boundary instead of using numerical flux.

Adding the ablation workshop test case 2.1 as an example (wall-only, no flow).

Tests and more examples will be added with the porous wall coupling (initial work in #900, more PRs to come ).

Questions for the review:

Is the scope and purpose of the PR clear?
- The PR should have a description.
- The PR should have a guide if needed (e.g., an ordering).
Is every top-level method and class documented? Are things that should be documented actually so?
Is the interface understandable? (I.e. can someone figure out what stuff does?) Is it well-defined?
Does the implementation do what the docstring claims?
Is everything that is implemented covered by tests?
Do you see any immediate risks or performance disadvantages with the design? Example: what do interface normals attach to?

Questions for the review (@MTCam):

Is the scope and purpose of the PR clear?
- The PR should have a description.
- The PR should have a guide if needed (e.g., an ordering).
Is every top-level method and class documented? Are things that should be documented actually so?
Is the interface understandable? (I.e. can someone figure out what stuff does?) Is it well-defined?
Does the implementation do what the docstring claims?
Is everything that is implemented covered by tests?
Do you see any immediate risks or performance disadvantages with the design? Example: what do interface normals attach to?

majosm

Added some more thoughts.

examples/ablation-workshop-mpi.py

mirgecom/materials/tacot.py

mirgecom/multiphysics/phenolics.py

mirgecom/wall_model.py

…calls

MTCam · 2023-07-12T16:11:41Z

mirgecom/materials/carbon_fiber.py

+        """Numerical data of thermal conductivity evaluated by PUMA."""
+        # FIXME the fully decomposed conductivity is given by the air, which
+        # is already accounted for in our model.
+        return 0.0988*progress**2 - 0.2751*progress + 0.201


It is my opinion that the constants in this class should be class data items with names. This would allow, for example, one to create a material with this exact model but with changed constants and coefficients.

MTCam · 2023-07-12T16:14:32Z

mirgecom/materials/tacot.py

+        gas_data = np.array([
+            [200.00,  21.996,  1.5119,  1.3334,  -7246.50,  0.086881],
+            [350.00,  21.995,  1.7259,  1.2807,  -7006.30,  0.144380],
+            [500.00,  21.948,  2.2411,  1.2133,  -6715.20,  0.196150],
+            [650.00,  21.418,  4.3012,  1.1440,  -6265.70,  0.243230],
+            [700.00,  20.890,  6.3506,  1.1242,  -6004.60,  0.258610],
+            [750.00,  19.990,  9.7476,  1.1131,  -5607.70,  0.274430],
+            [800.00,  18.644,  14.029,  1.1116,  -5014.40,  0.290920],
+            [850.00,  17.004,  17.437,  1.1171,  -4218.50,  0.307610],
+            [900.00,  15.457,  17.009,  1.1283,  -3335.30,  0.323490],
+            [975.00,  14.119,  8.5576,  1.1620,  -2352.90,  0.344350],
+            [1025.0,  13.854,  4.7840,  1.1992,  -2034.20,  0.356630],
+            [1100.0,  13.763,  3.5092,  1.2240,  -1741.20,  0.373980],
+            [1150.0,  13.737,  3.9008,  1.2087,  -1560.90,  0.385360],
+            [1175.0,  13.706,  4.8067,  1.1899,  -1453.50,  0.391330],
+            [1200.0,  13.639,  6.2353,  1.1737,  -1315.90,  0.397930],
+            [1275.0,  13.256,  8.4790,  1.1633,  -739.700,  0.421190],
+            [1400.0,  12.580,  9.0239,  1.1583,  353.3100,  0.458870],
+            [1525.0,  11.982,  11.516,  1.1377,  1608.400,  0.483230],
+            [1575.0,  11.732,  12.531,  1.1349,  2214.000,  0.487980],
+            [1625.0,  11.495,  11.514,  1.1444,  2826.800,  0.491950],
+            [1700.0,  11.255,  7.3383,  1.1849,  3529.400,  0.502120],
+            [1775.0,  11.139,  5.3118,  1.2195,  3991.000,  0.516020],
+            [1925.0,  11.046,  4.2004,  1.2453,  4681.800,  0.545280],
+            [2000.0,  11.024,  4.0784,  1.2467,  4991.300,  0.559860],
+            [2150.0,  10.995,  4.1688,  1.2382,  5605.400,  0.588820],
+            [2300.0,  10.963,  4.5727,  1.2214,  6257.300,  0.617610],
+            [2450.0,  10.914,  5.3049,  1.2012,  6993.500,  0.646380],
+            [2600.0,  10.832,  6.4546,  1.1815,  7869.600,  0.675410],
+            [2750.0,  10.701,  8.1450,  1.1650,  8956.900,  0.705000],
+            [2900.0,  10.503,  10.524,  1.1528,  10347.00,  0.735570],
+            [3050.0,  10.221,  13.755,  1.1449,  12157.00,  0.767590],
+            [3200.0,  9.8394,  17.957,  1.1408,  14523.00,  0.801520],
+            [3350.0,  9.3574,  22.944,  1.1401,  17584.00,  0.837430],
+        ])


Any chance we'd be amenable to specifying this data at class instantiation time after reading it from a file or hard-coding it only in the driver?

majosm

Just a couple of minor suggestions, otherwise LGTM. 👍

majosm · 2023-07-12T17:39:47Z

doc/operators/composite_materials.rst

+    Note that :mod:`pyrometheus` is used to handle the species properties.
+
+Composite Materials
+^^^^^^^^^^^^^^^^^^^
+
+    This section covers the response of composite materials made of phenolic
+    resin and carbon fibers.
+    The carbon fibers are characterized as a highly porous material,
+    with void fraction $\approx 90\%$. Phenolic resins are added to rigidize
+    the fibers by permeating the material and filling partially the gaps between
+    fibers. As the material is heated up by the flow, the resin pyrolysis, i.e.,
+    it degrades and produces gaseous species.
+
+    The temporal evolution of wall density is solved in
+    order to predict the material degradation. As the
+    :class:`~mirgecom.materials.tacot.Pyrolysis` progresses, the mass of each 
+    $i$ constituents of the resin, denoted by
+    :attr:`~mirgecom.wall_model.PorousFlowDependentVars.material_densities`,
+    is calculated as
+
+    .. math ::
+        \frac{\partial \rho_i}{\partial t} = \dot{\omega}_i \mbox{ .}
+
+    This process creates gaseous species following
+
+    .. math ::
+        \frac{\partial \rho_g}{\partial t} + \nabla \cdot \rho_g \mathbf{u} = 
+            - \sum_i \dot{\omega}_i \mbox{ ,}
+
+    where the
+    :attr:`~mirgecom.fluid.ConservedVars.mass`
+    is $\rho_g$. The source term indicates that all solid resin must become
+    gas in order to satisfy mass conservation. Lastly, the gas velocity
+    $\mathbf{u}$ is obtained by Darcy's law, given by
+
+    .. math::
+        \mathbf{u} = \frac{\mathbf{K}}{\mu \epsilon} \cdot \nabla P \mbox{ .}
+
+    In this equation, $\mathbf{K}$ is the second-order permeability tensor,
+    $\mu$ is the gas viscosity, $\epsilon$ is the void fraction and $P$ is
+    the gas pressure.
+
+    The
+    :attr:`~mirgecom.fluid.ConservedVars.energy`
+    of the bulk material is given by
+
+    .. math::
+        \frac{\partial \rho_b e_b}{\partial t}
+        + \nabla \cdot (\epsilon_{g} \rho_g h_g \mathbf{u})
+        = \nabla \cdot \left( \bar{\boldsymbol{\kappa}} \nabla T \right)
+        + \mu \epsilon_{g}^2 (\bar{\mathbf{K}}^{-1} \cdot \vec{v} ) \cdot \vec{v}
+        \mbox{ .}
+
+    The first term on the RHS is modeled as an effective diffusive transfer
+    using Fourier's law, where the thermal conductivity is given by
+    $\bar{\boldsymbol{\kappa}}$. The second term on the RHS account for the
+    viscous dissipation in the Darcy's flow. The sub-index $b$ is the bulk
+    material consisting of both solid and gas phases, where
+
+    .. math::
+        \rho_b e_b = \epsilon_{g} \rho_g e_g + \epsilon_s \rho_s e_s \mbox{ .}
+
+    The energy of the gas is given by $e_g(T) = h_g(T) - \frac{P}{\rho_g}$,
+    where $h_g$ is its enthalpy.
+
+    From the conserved variables, it is possible to compute the decomposition
+    status, denoted by
+    :attr:`~mirgecom.wall_model.WallDependentVars.tau`.
+    This yields the proportion of virgin (unpyrolyzed material) to char (fully
+    pyrolyzed) and, consequently, the different thermophysicochemical
+    properties of the solid phase. Thus, the instantaneous material properties
+    depend on the current state of the material, as well as the 
+    :attr:`~mirgecom.eos.GasDependentVars.temperature`.
+    It is evaluated using Newton iteration based on both
+    :attr:`~mirgecom.materials.tacot.GasProperties.enthalpy` (tabulated
+    data) or 
+    :attr:`~mirgecom.eos.PyrometheusMixture.get_internal_energy` (Pyrometheus)
+    and
+    :attr:`~mirgecom.wall_model.WallEOS.enthalpy`,
+    as well as their respective derivatives, namely
+    :attr:`~mirgecom.materials.tacot.GasProperties.heat_capacity` and
+    :attr:`~mirgecom.wall_model.WallEOS.heat_capacity`.
+
+    In *MIRGE-Com*, the solid properties are obtained by fitting polynomials
+    to tabulated data for easy evaluation of the properties based on the
+    temperature. The complete list of properties can be find, for instance, in
+    :mod:`~mirgecom.materials.tacot`.
+    Different materials can be incorporated as separated files.
+
+    The
+    :class:`~mirgecom.materials.tacot.GasProperties` are obtained based on
+    tabulated data, assuming chemical equilibrium, and evaluated with splines
+    for interpolation of the entries. However, the data is currently obtained
+    by PICA.
+
+.. important ::
+
+    The current implementation follows the description of [Lachaud_2014]_ 
+    for type 2 code. Additional details, extensive formulation and references
+    are provided in https://github.com/illinois-ceesd/phenolics-notes


Should this go at the top instead of here? In the built docs it ends up towards the bottom, but it seems too important to be down there? Just a thought.

majosm · 2023-07-12T17:46:46Z

mirgecom/multiphysics/phenolics.py

+def initializer(dcoll, gas_model, material_densities, temperature,
+                gas_density=None, pressure=None):


I would either: 1) add an optional dd=DD_VOLUME_ALL argument (even though it's not currently used), or 2) remove the dcoll and pass dim directly instead.

tulioricci added 30 commits March 28, 2023 16:44

The very first of many commits...

7031511

Reorganizing files

154e3a5

Initial draft on dependent vars

2b7f29b

Newton iteration for solid temperature

c08820d

Heat conduction + presc.flux + pyrolysis = OK

73d88cc

Changes to make lazy work

a6ac753

Merge remote-tracking branch 'origin/main' into tulio/phenolics-v0

2493c71

Merge main

da844a5

First attempt at gas data interpolation

fcc7ce2

First attempt at gas data interpolation

e9337eb

Spline + where seems to be working

7b0417b

Merge main

2df138d

Flake8 + pydocstyle + fix tau/progress

e10cc33

Several minor changes

f5a541d

First attempt at Darcy Equation + SSPRK43; add docs; perfume

0cba99b

'Darcy + phenolics' seems to be working just fine

8392189

Improving documentation...

41d69df

Failed attempt to add inviscid RHS

08e4971

Playing around with the inviscid term

bb8d35e

Inviscid RHS seems to be working

6aaed1d

Add operator to its own file

6912f87

Merge main

968bcb2

Merge remote-tracking branch 'origin/main' into tulio/phenolics-v0

c8cb858

Remove aux file from gas_data

8371a21

Remove operator from driver and pass to operator.py

fb31c3d

Fix bug in driver, simplify diff and remove dummy line in ssprk

bcafcb9

Adding missing features to gas.py and phenolics.py

c09f2ed

Minor changes

a2c92a7

Use greater_than

46f3bb5

Debugging phenolics and add 1D file

931031d

tulioricci and others added 10 commits May 25, 2023 11:22

Merge branch 'main' into tulio/phenolics-v0

99983b6

Remove (now) dummy functions to the wall model

c6abae2

Remove (now) dummy functions to the wall model and fix aw.py example

57d6a8a

Merge branch 'main' into tulio/phenolics-v0

356e87d

Merge branch 'main' into tulio/phenolics-v0

d84f35b

Merge main

4f16476

Add missing np to gas_model

82b5d32

Merge branch 'main' into tulio/phenolics-v0

154e114

Merge branch 'main' into tulio/phenolics-v0

7e0331f

Cosmetic changes to gas and wall_model

8e0a910

tulioricci requested review from majosm and MTCam July 4, 2023 01:34

majosm reviewed Jul 11, 2023

View reviewed changes

tulioricci added 7 commits July 11, 2023 12:41

Merge main

cf6b35a

Changes per Matts review

bc52690

Use fancy way to give the path and change PrescDiffFlux

efa834b

flake8; rename wall_density; modify example to use pre/post/advancer …

b9d8385

…calls

Move back the table file and fix wall_density

22e1bb2

Test Matts code to see if the Bprime table works in CI

4861458

Merge branch 'main' into tulio/phenolics-v0

cc0a087

This was referenced Jul 12, 2023

De-fluidize material data for porous media flow #930

Open

Recursion error in ablation example #931

Closed

MTCam reviewed Jul 12, 2023

View reviewed changes

majosm approved these changes Jul 12, 2023

View reviewed changes

tulioricci merged commit 680968d into main Jul 12, 2023

tulioricci deleted the tulio/phenolics-v0 branch July 12, 2023 20:27

tulioricci mentioned this pull request Aug 21, 2023

Reorganize a bit the structure of the wall files #935

Merged

8 tasks

tulioricci mentioned this pull request Mar 11, 2024

Update material functions and fix AW example #1019

Merged

8 tasks

tulioricci mentioned this pull request Mar 26, 2024

Add unified-domain to solve porous flow #1001

Draft

8 tasks

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Add wall-only part of Y3 phenolics model #875

Add wall-only part of Y3 phenolics model #875

tulioricci commented Apr 13, 2023 •

edited by majosm

Loading

majosm left a comment

MTCam Jul 12, 2023

MTCam Jul 12, 2023

majosm left a comment

majosm Jul 12, 2023

majosm Jul 12, 2023

		def initializer(dcoll, gas_model, material_densities, temperature,
		gas_density=None, pressure=None):

Add wall-only part of Y3 phenolics model #875

Add wall-only part of Y3 phenolics model #875

Conversation

tulioricci commented Apr 13, 2023 • edited by majosm Loading

majosm left a comment

Choose a reason for hiding this comment

MTCam Jul 12, 2023

Choose a reason for hiding this comment

MTCam Jul 12, 2023

Choose a reason for hiding this comment

majosm left a comment

Choose a reason for hiding this comment

majosm Jul 12, 2023

Choose a reason for hiding this comment

majosm Jul 12, 2023

Choose a reason for hiding this comment

tulioricci commented Apr 13, 2023 •

edited by majosm

Loading