From 7671a89e142450fb223ae01b4b9f23ddf70a8f01 Mon Sep 17 00:00:00 2001 From: Nicolas Kruse Date: Fri, 6 Jun 2025 13:15:34 +0200 Subject: [PATCH] Examples added, example readme and example unit test updated --- docs/source/render_examples.py | 2 +- examples/README.md | 2 +- examples/carbon_activity.md | 68 +++++++++++++++++++++++++++ examples/methane_mixtures.md | 60 +++++++++++++++++++++++ examples/soec_methane.md | 39 ++++++++++----- examples/sulfur_oxygen_equalibrium.md | 52 ++++++++++++++++++++ tests/md_to_code.py | 5 +- tests/test_doc_examples.py | 4 +- 8 files changed, 215 insertions(+), 17 deletions(-) create mode 100644 examples/carbon_activity.md create mode 100644 examples/methane_mixtures.md create mode 100644 examples/sulfur_oxygen_equalibrium.md diff --git a/docs/source/render_examples.py b/docs/source/render_examples.py index ce33166..d38ad0b 100644 --- a/docs/source/render_examples.py +++ b/docs/source/render_examples.py @@ -16,7 +16,7 @@ def run_cmd(command: list[str]): assert (not result.stderr or any('RuntimeWarning: ' in line for line in result.stderr.splitlines()) or - any('[NbConvertApp]' in line for line in result.stderr.splitlines())), 'ERROR: ' + result.stderr + any('[NbConvertApp]' in line and 'error' not in line.lower() for line in result.stderr.splitlines())), 'ERROR: ' + result.stderr def run_rendering(input_path: str, output_directory: str): diff --git a/examples/README.md b/examples/README.md index a77d071..ba40e2f 100644 --- a/examples/README.md +++ b/examples/README.md @@ -21,5 +21,5 @@ jupyter nbconvert --to markdown docs/source/_autogenerated/soec_methane.ipynb -- A new example Markdown file can be created from a Jupyter Notebook running the following command: ``` bash -jupyter nbconvert --to new_example.ipynb --NbConvertApp.use_output_suffix=False --ClearOutputPreprocessor.enabled=True --output-dir examples/ --output new_example.md +jupyter nbconvert --to markdown new_example.ipynb --NbConvertApp.use_output_suffix=False --ClearOutputPreprocessor.enabled=True --output-dir examples/ --output new_example.md ``` \ No newline at end of file diff --git a/examples/carbon_activity.md b/examples/carbon_activity.md new file mode 100644 index 0000000..4e8c67e --- /dev/null +++ b/examples/carbon_activity.md @@ -0,0 +1,68 @@ +# Carbon Activity + +This example shows the equilibrium calculation for solid carbon. + + +```python +import gaspype as gp +import numpy as np +import matplotlib.pyplot as plt +``` + +Setting temperatures and pressure: + + +```python +t_range = np.array([600, 700, 800, 900, 1100, 1500]) # °C + +p = 1e5 # Pa + +fs = gp.fluid_system(['H2', 'H2O', 'CO2', 'CO', 'CH4']) +``` + +Equilibrium calculation for methane steam mixtures: + + +```python +ratio = np.linspace(0.01, 1.5, num=128) + +fl = gp.fluid({'CH4': 1}, fs) + ratio * gp.fluid({'H2O': 1}, fs) +``` + +gaspype.carbon_activity supports currently only 0D fluids therefore we build this helper function: + + +```python +def partial_c_activity(fl: gp.fluid, t: float, p: float): + fls = fl.array_composition.shape + + eq_fl = gp.equilibrium(fl, t, p) + + ret = np.zeros(fls[0]) + for i in range(fls[0]): + ret[i] = gp.carbon_activity(gp.fluid(eq_fl.array_composition[i,:], fs), t, p) + + return ret +``` + +Now we use the helper function to calculate the carbon activitx for all +compositions in equilibrium_h2o times all temperatures in t_range: + + +```python +carbon_activity = np.vstack([partial_c_activity(fl, tc + 273.15, p) for tc in t_range]) +``` + +Plot carbon activities, a activity of > 1 means there is thermodynamically the formation of sold carbon favored. + + +```python +fig, ax = plt.subplots(figsize=(6, 4), dpi=120) +ax.set_xlabel("CO2/CH4") +ax.set_ylabel("carbon activity") +ax.set_ylim(1e-1, 1e3) +ax.set_yscale('log') +ax.plot(ratio, carbon_activity.T) +ax.hlines(1, np.min(ratio), np.max(ratio), colors='k', linestyles='dashed') +ax.legend([f'{tc} °C' for tc in t_range]) +``` diff --git a/examples/methane_mixtures.md b/examples/methane_mixtures.md new file mode 100644 index 0000000..a39268d --- /dev/null +++ b/examples/methane_mixtures.md @@ -0,0 +1,60 @@ +# Methane Mixtures + +This example shows equilibria of methane mixed with steam and CO2 + + +```python +import gaspype as gp +import numpy as np +import matplotlib.pyplot as plt +``` + +Setting temperature and pressure: + + +```python +t = 900 + 273.15 +p = 1e5 + +fs = gp.fluid_system(['H2', 'H2O', 'CO2', 'CO', 'CH4', 'O2']) +``` + +Equilibrium calculation for methane steam mixtures: + + +```python +ratio = np.linspace(0.01, 1.5, num=64) + +el = gp.fluid({'CH4': 1}, fs) + ratio * gp.fluid({'H2O': 1}, fs) +equilibrium_h2o = gp.equilibrium(el, t, p) +``` + + +```python +fig, ax = plt.subplots(figsize=(6, 4), dpi=120) +ax.set_xlabel("H2O/CH4") +ax.set_ylabel("molar fraction") +ax.set_ylim(0, 1.1) +#ax.set_xlim(0, 100) +ax.plot(ratio, equilibrium_h2o.get_x()) +ax.legend(fs.active_species) +``` + +Equilibrium calculation for methane CO2 mixtures: + + +```python +el = gp.fluid({'CH4': 1}, fs) + ratio * gp.fluid({'H2O': 1}, fs) +equilibrium_co2 = gp.equilibrium(el, t, p) +``` + + +```python +fig, ax = plt.subplots(figsize=(6, 4), dpi=120) +ax.set_xlabel("CO2/CH4") +ax.set_ylabel("molar fraction") +ax.set_ylim(0, 1.1) +#ax.set_xlim(0, 100) +ax.plot(ratio, equilibrium_co2.get_x()) +ax.legend(fs.active_species) +``` diff --git a/examples/soec_methane.md b/examples/soec_methane.md index d7d0e7e..a2814f4 100644 --- a/examples/soec_methane.md +++ b/examples/soec_methane.md @@ -1,4 +1,8 @@ -# SOEC example +# SOEC with Methane + +This example shows a 1D isothermal SOEC (Solid oxide electrolyzer cell) model. + +The operating parameters chosen here are not necessary realistic ```python import gaspype as gp @@ -7,7 +11,8 @@ import numpy as np import matplotlib.pyplot as plt ``` - +Calculation of the local equilibrium compositions on the fuel and air +side in counter flow along the fuel flow direction: ```python fuel_utilization = 0.90 air_utilization = 0.5 @@ -28,9 +33,8 @@ fuel_side = gp.equilibrium(feed_fuel + perm_oxygen, t, p) air_side = gp.equilibrium(feed_air - perm_oxygen, t, p) ``` - +Plot compositions of the fuel and air side: ```python -#Plot compositions on fuel and air side fig, ax = plt.subplots() ax.set_xlabel("Conversion") ax.set_ylabel("Molar fraction") @@ -44,13 +48,14 @@ ax.plot(conversion, air_side.get_x(), '-') ax.legend(air_side.species) ``` +Calculation of the oxygen partial pressures: ```python o2_fuel_side = gp.oxygen_partial_pressure(fuel_side, t, p) o2_air_side = air_side.get_x('O2') * p ``` +Plot oxygen partial pressures: ```python -#Plot oxygen partial pressure fig, ax = plt.subplots() ax.set_xlabel("Conversion") ax.set_ylabel("Oxygen partial pressure / Pa") @@ -59,13 +64,14 @@ ax.plot(conversion, np.stack([o2_fuel_side, o2_air_side], axis=1), '-') ax.legend(['o2_fuel_side', 'o2_air_side']) ``` +Calculation of the local nernst potential between fuel and air side: ```python z_O2 = 4 nernst_voltage = R*t / (z_O2*F) * np.log(o2_air_side/o2_fuel_side) ``` +#Plot nernst potential: ```python -#Plot voltage potential fig, ax = plt.subplots() ax.set_xlabel("Conversion") ax.set_ylabel("Voltage / V") @@ -73,24 +79,33 @@ ax.plot(conversion, nernst_voltage, '-') print(np.min(nernst_voltage)) ``` +The model uses between each node a constant conversion. Because +current density depends strongly on the position along the cell +the constant conversion does not relate to a constant distance. + +![Alt text](../../media/soc_inverted.svg) + +To calculate the local current density (**node_current**) as well +as the total cell current (**terminal_current**) the (relative) +physical distance between the nodes (**dz**) must be calculated: ```python cell_voltage = 0.77 #V ASR = 0.2 #Ohm*cm² -node_current = (nernst_voltage - cell_voltage) / ASR #mA/cm² (Current density at each node) +node_current = (nernst_voltage - cell_voltage) / ASR # mA/cm² (Current density at each node) -current = (node_current[1:] + node_current[:-1]) / 2 #mA/cm² (Average current density between the nodes) +current = (node_current[1:] + node_current[:-1]) / 2 # mA/cm² (Average current density between the nodes) -dz = 1/current / np.sum(1/current) #Relative distance between each node +dz = 1/current / np.sum(1/current) # Relative distance between each node -terminal_current = np.sum(current * dz) #mA/cm² (Total cell current per cell area) +terminal_current = np.sum(current * dz) # mA/cm² (Total cell current per cell area) print(f'Terminal current: {terminal_current:.2f} A/cm²') ``` +Plot the local current density: ```python -#Plot current density -z_position = np.concatenate([[0], np.cumsum(dz)]) #Relative position of each node +z_position = np.concatenate([[0], np.cumsum(dz)]) # Relative position of each node fig, ax = plt.subplots() ax.set_xlabel("Relative cell position") diff --git a/examples/sulfur_oxygen_equalibrium.md b/examples/sulfur_oxygen_equalibrium.md new file mode 100644 index 0000000..4831ec3 --- /dev/null +++ b/examples/sulfur_oxygen_equalibrium.md @@ -0,0 +1,52 @@ +# Sulfur Oxygen Equilibrium + +This example shows equilibrium calculations for sulfur/oxygen mixtures. + +```python +import gaspype as gp +import numpy as np +import matplotlib.pyplot as plt +``` + +List possible sulfur/oxygen species: + +```python +gp.species(element_names = 'S, O') +``` + +Or more specific by using regular expressions: + + +```python +gp.species('S?[2-3]?O?[2-5]?', use_regex=True) +``` + +Calculation of the molar equilibrium fractions for sulfur and oxygen depending on the oxygen to sulfur ratio: + + +```python +fs = gp.fluid_system(['S2', 'S2O', 'SO2', 'SO3', 'O2']) + +oxygen_ratio = np.linspace(0.5, 3, num=128) +el = gp.elements({'S': 1}, fs) + oxygen_ratio * gp.elements({'O': 1}, fs) + +composition = gp.equilibrium(el, 800+273.15, 1e4) + +plt.plot(oxygen_ratio, composition.get_x()) +plt.legend(composition.species) +``` + +Calculation of the molar equilibrium fractions for sulfur and oxygen depending on temperature in °C: + + +```python +fs = gp.fluid_system(['S2', 'S2O', 'SO2', 'SO3', 'O2']) + +el = gp.elements({'S': 1, 'O':2.5}, fs) + +t_range = np.linspace(500, 1300, num=32) +composition = gp.equilibrium(el, t_range+273.15, 1e4) + +plt.plot(t_range, composition.get_x()) +plt.legend(composition.species) +``` diff --git a/tests/md_to_code.py b/tests/md_to_code.py index bb009e7..2d9c7dc 100644 --- a/tests/md_to_code.py +++ b/tests/md_to_code.py @@ -93,7 +93,10 @@ def segments_to_test(segments: Iterable[markdown_segment], script_language: str if segment.code_block: if segment.language == script_language: lines = [line for line in segment.text.splitlines() if line.strip()] - ret_block_flag = lines[-1] if not re.match(r'^[^(]*=', lines[-1]) and not lines[-1].startswith('import ') else None + ret_block_flag = lines[-1] if (not re.match(r'^[^(]*=', lines[-1]) and + not lines[-1].startswith('import ') and + not lines[-1].startswith('from ') and + not lines[-1].startswith(' ')) else None # print('Last line: ', ret_block_flag, '-----------', lines[-1]) yield '' diff --git a/tests/test_doc_examples.py b/tests/test_doc_examples.py index 0b498de..61222a4 100644 --- a/tests/test_doc_examples.py +++ b/tests/test_doc_examples.py @@ -11,7 +11,7 @@ def test_readme(): def test_example_code(): - filter = 'docs/source/examples/*.md' + filter = 'examples/*.md' files = glob(filter) for path in files: @@ -19,7 +19,7 @@ def test_example_code(): if not file_name.lower() == 'readme': print(f"> Test Example {file_name} ...") md_to_code.convert_to('test', path, f'tests/autogenerated_{file_name}.py') - mod = importlib.import_module(file_name) + mod = importlib.import_module(f'autogenerated_{file_name}') mod.run_test()