.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "api/gallery/heat_transfer/phase_change/plot_01_stefan_problem.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code or to run this example in your browser via JupyterLite. .. rst-class:: sphx-glr-example-title .. _sphx_glr_api_gallery_heat_transfer_phase_change_plot_01_stefan_problem.py: Stefan's problem: a moving freezing front (1891) ================================================ In 1891 Josef Stefan asked how fast sea ice thickens in the polar winter. The ice conducts heat from the freezing water below to the cold air above, and every millimetre of new ice releases latent heat at the moving interface, so the boundary is part of the unknown. Franz Neumann had found the exact similarity solution in his lectures of the 1860s: the front advances as :math:`s = 2\lambda\sqrt{\alpha t}` with :math:`\lambda e^{\lambda^2}\operatorname{erf}\lambda = Ste/\sqrt\pi` (:func:`~engineeringkit.heat_transfer.neumann_stefan_front`, :func:`~engineeringkit.heat_transfer.stefan_number`). This example grows a winter of sea ice, checks the energy balance at the front, and compares with Stefan's own quasi-steady estimate. .. GENERATED FROM PYTHON SOURCE LINES 19-25 .. code-block:: Python import matplotlib.pyplot as plt import numpy as np from scipy import special from engineeringkit.heat_transfer import neumann_stefan_front, stefan_number .. GENERATED FROM PYTHON SOURCE LINES 26-32 Sea ice under a -20 °C sky -------------------------- Ice: :math:`k = 2.22` W/(m K), :math:`\rho = 917` kg/m³, :math:`c_p = 2050` J/(kg K), latent heat 334 kJ/kg. The top surface is held 20 K below the freezing point. .. GENERATED FROM PYTHON SOURCE LINES 32-45 .. code-block:: Python k, rho, cp, latent, dT = 2.22, 917.0, 2050.0, 334e3, 20.0 alpha = k / (rho * cp) Ste = stefan_number(cp, dT, latent) day = 86400.0 days = np.array([1, 10, 30, 100]) s = neumann_stefan_front(days * day, alpha, Ste) lam = s[0] / (2 * np.sqrt(alpha * day)) print(f"Ste = {Ste:.4f}, lambda = {lam:.4f}") for d, thickness in zip(days, s): print(f"after {d:3d} days: {thickness * 100:5.1f} cm of ice") assert np.isclose(lam * np.exp(lam**2) * special.erf(lam), Ste / np.sqrt(np.pi)) .. rst-class:: sphx-glr-script-out .. code-block:: none Ste = 0.1228, lambda = 0.2429 after 1 days: 15.5 cm of ice after 10 days: 49.1 cm of ice after 30 days: 85.0 cm of ice after 100 days: 155.2 cm of ice .. GENERATED FROM PYTHON SOURCE LINES 46-52 Latent heat released equals heat conducted ------------------------------------------ At the front, :math:`k\,\partial T/\partial x = \rho h_{sf}\,\mathrm{d}s/\mathrm{d}t`. Inside the ice the temperature is :math:`T_s + \Delta T\operatorname{erf}(x/2\sqrt{\alpha t})/\operatorname{erf}\lambda`. .. GENERATED FROM PYTHON SOURCE LINES 52-68 .. code-block:: Python t, dt = 30 * day, 1.0 front = neumann_stefan_front(t, alpha, Ste) def ice_temperature(x, t): return -dT + dT * special.erf(x / (2 * np.sqrt(alpha * t))) / special.erf(lam) # °C, 0 at the front h = 1e-6 conducted = k * (ice_temperature(front, t) - ice_temperature(front - h, t)) / h released = rho * latent * (neumann_stefan_front(t + dt, alpha, Ste) - neumann_stefan_front(t - dt, alpha, Ste)) / (2 * dt) print(f"at 30 days: conducted {conducted:.3f} W/m², latent heat released {released:.3f} W/m²") assert np.isclose(ice_temperature(front, t), 0.0, atol=1e-12) assert np.isclose(conducted, released, rtol=1e-4) .. rst-class:: sphx-glr-script-out .. code-block:: none at 30 days: conducted 50.216 W/m², latent heat released 50.216 W/m² .. GENERATED FROM PYTHON SOURCE LINES 69-76 Stefan's quasi-steady estimate ------------------------------ Neglecting the heat stored in the ice (a linear profile) gives :math:`s = \sqrt{2k\Delta T t/(\rho h_{sf})} = \sqrt{2\,Ste\,\alpha t}`, Stefan's law of ice growth. It overestimates by a few percent at this Stefan number and becomes exact as :math:`Ste \to 0`. .. GENERATED FROM PYTHON SOURCE LINES 76-99 .. code-block:: Python quasi = np.sqrt(2 * Ste * alpha * days * day) print(f"quasi-steady overestimate: {(quasi / s - 1)[0]:.2%}") for small in [1e-2, 1e-4]: ratio = np.sqrt(2 * small * alpha * day) / neumann_stefan_front(day, alpha, small) assert abs(ratio - 1) < small assert 0 < quasi[0] / s[0] - 1 < 0.05 fig, axes = plt.subplots(1, 2, figsize=(10, 3.8)) winter = np.linspace(0, 180, 200) * day axes[0].plot(winter / day, neumann_stefan_front(winter, alpha, Ste) * 100, label="Neumann (exact)") axes[0].plot(winter / day, np.sqrt(2 * Ste * alpha * winter) * 100, "--", label="Stefan quasi-steady") axes[0].set(xlabel="days", ylabel="ice thickness [cm]", title="Growth as the square root of time") axes[0].legend() for d in [10, 30, 100]: x = np.linspace(0, neumann_stefan_front(d * day, alpha, Ste), 50) axes[1].plot(ice_temperature(x, d * day), x * 100, label=f"{d} days") axes[1].invert_yaxis() axes[1].set(xlabel="ice temperature [°C]", ylabel="depth below the surface [cm]", title="Profiles in the ice") axes[1].legend() fig.tight_layout() plt.show() .. image-sg:: /api/gallery/heat_transfer/phase_change/images/sphx_glr_plot_01_stefan_problem_001.png :alt: Growth as the square root of time, Profiles in the ice :srcset: /api/gallery/heat_transfer/phase_change/images/sphx_glr_plot_01_stefan_problem_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none quasi-steady overestimate: 1.99% .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 0.062 seconds) .. _sphx_glr_download_api_gallery_heat_transfer_phase_change_plot_01_stefan_problem.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: lite-badge .. image:: images/jupyterlite_badge_logo.svg :target: ../../../../lite/lab/index.html?path=api/gallery/heat_transfer/phase_change/plot_01_stefan_problem.ipynb :alt: Launch JupyterLite :width: 150 px .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_01_stefan_problem.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_01_stefan_problem.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_01_stefan_problem.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_