Add eich profile

documentation/source/physics-models/plasma_scrape_off_layer.md

@@ -25,7 +25,43 @@ $$
 A_{\parallel,u} = 2\pi\lambda_{\text{q,u}}R_{\text{u}}\frac{B_{\text{p,u}}}{B_{\text{Tot,u}}}
 $$
 
+---------------
 
+## Upstream radial decay | `calculate_outboard_midplane_near_sol_radial_profile()`
+
+The radial decay length $\lambda_{\text{q}}$ of the scrape-off layer (SOL) at the outer midplane of a tokamak is defined as the e-folding distance over which plasma heat and particle fluxes decay exponentially outside the last closed flux surface. Therefore the decay of the total heat flux outside the separatrix towards the vessel walls can be modelled as [^eich_2011] [^eich_2013]:
+
+$$
+q_{\text{u}}(r) = q_{\parallel,\text{u}}e^{\frac{-r}{\lambda_{\text{q}}}}
+$$
+
+where $r = R - R_{\text{sep}}$, $R_{\text{sep}}$ being the major radius of the separatrix, $\lambda_{\text{q}}$ the [power decay length](#power-decay-lengths) and $q_{\parallel}$ the [upstream energy flux density](#upstream-radial-decay--calculate_outboard_midplane_near_sol_radial_profile)
+
+----------------
+
+## Eich parallel flux at target | `calculate_eich_target_heat_flux_profile()`
+
+The Eich formula (often called the standard SOL heat flux profile) is the primary mathematical model used to describe the distribution of heat target loads on tokamak divertor plates. It convolutionally connects the physics of the plasma edge at the outer midplane with the geometric projection of the heat hitting the divertor surface [^eich_2011] [^eich_2013].
+
+Heat transport into the private flux region is modeled by convolving the power profile $q_{\text{u}}(r)$ with a Gaussian function of width $S$ known as the [spreading parameter](#spreading-parameter). 
+
+$$
+q_{\parallel,t}(s) = \frac{q_0}{2}\times \exp\left[\left(\frac{S}{2\lambda_{\text{q}}f_x}\right)^2- \frac{s-s_0}{\lambda_q f_x}\right] \times \operatorname{erfc}\left(\frac{S}{2\lambda_{\text{q}}f_x}- \frac{s-s_0}{S}\right) + q_{\text{BG}}
+$$
+
+where $s$ is the coordinate along the divertor target, $s_0$ is the strike-point location on the target, $\operatorname{erfc}$ is the complementary error function, $q_{\text{BG}}$ is the background heat flux, $\lambda_{\text{q}}$ is the [power decay length](#power-decay-lengths), $f_x$ is the effective flux expansion in the region,
+
+A compact equivalent form is:
+
+$$
+q_{\parallel,t}(\overline{s}) = \frac{q_0}{2}\times \exp\left[\left(\frac{S}{2\lambda_{\text{q}}f_x}\right)^2- \frac{\overline{s}}{\lambda_q f_x}\right] \times \operatorname{erfc}\left(\frac{S}{2\lambda_{\text{q}}f_x}- \frac{\overline{s}}{S}\right) + q_{\text{BG}}
+$$
+
+The connection to upstream midplane coordinates is usually:
+
+$$
+\overline{s} = f_x(R-R_{\text{sep}})
+$$
 
 ------------------
 
@@ -83,11 +119,42 @@ The $R^2$ value for this fit is 0.55
 
 --------
 
+## Spreading Parameter
+
+The scrape-off layer (SOL) spreading parameter $S$ represents a Gaussian width that quantifies additional perpendicular heat spreading in the divertor leg. It works alongside the upstream heat flux decay length $\lambda_{q}$ to determine total target heat loads on the divertor.
+
+Unlike $\lambda_{q}$, which is governed by robust upstream parallel and perpendicular transport physics at the plasma midplane, $S$ is inherently a "local" divertor parameter. Deriving a single, absolute multi-machine formula for $S$ is incredibly difficult due to several overlapping regional variables:
+
+- Divertor Geometry: The path length from the X-point to the target tile heavily impacts how much the heat spreads radially.
+
+- Plasma Recycling Regimes: Low-recycling, high-recycling, and detached plasma conditions completely alter the cross-field diffusion rates.
+
+- Localized Radiation: Impurity seeding and neutral gas interactions dissipate power unevenly along the divertor leg, altering the effective Gaussian profile width.
+
+-----------
+
+### Scarabosio 2015 | `calculate_scarabosio2015_power_spreading_factor()`
+
+The H-mode SOL spreading factor, $S$ is given in $\text{m}$ by[^scarabosio_2015]:
+
+$$
+S = (0.12(\pm0.07)\times 10^{-3}) P_{\text{sep}}^{0.21(\pm0.11)}R_0^{0.71(\pm0.5)}B_{\text{p}}(a)^{-0.82(\pm0.27)}n_{\text{sep}}^{0.71(\pm0.5)}
+$$
+
+- This was fitted from ASDEX Upgrade and JET outer target data
+- The $R^2$ value of the regression fit was 0.65
+
+------------
+
 [^eich_2013]: T. Eich et al., “Scaling of the tokamak near the scrape-off layer H-mode power width and implications for ITER,” Nuclear Fusion, vol. 53, no. 9 p. 093031, Aug. 2013, doi: 10.1088/0029-5515/53/9/093031.
 
 [^mast_2014]: A. J. Thornton and A. Kirk, “Scaling of the scrape-off layer width during inter-ELM H modes on MAST as measured by infrared thermography,”
 Plasma Physics and Controlled Fusion, vol. 56, no. 5, p. 055008, Apr. 2014, doi: 10.1088/0741-3335/56/5/055008.
 
 [^stangeby_boundary]: P. C. Stangeby, “The Plasma Boundary of Magnetic Fusion Devices,” Jan. 2000, doi: 10.1201/9780367801489.
 
-[^henderson_step]: S. S. Henderson et al., “An overview of the STEP divertor design and the simple models driving the plasma exhaust scenario,” Nuclear Fusion, vol. 65, no. 1, pp. 016033–016033, Nov. 2024, doi: 10.1088/1741-4326/ad93e7.
+[^henderson_step]: S. S. Henderson et al., “An overview of the STEP divertor design and the simple models driving the plasma exhaust scenario,” Nuclear Fusion, vol. 65, no. 1, pp. 016033–016033, Nov. 2024, doi: 10.1088/1741-4326/ad93e7.
+
+[^scarabosio_2015]: A. Scarabosio et al., “Scaling of the divertor power spreading (S-factor) in open and closed divertor operation in JET and ASDEX Upgrade,” Journal of Nuclear Materials, vol. 463, pp. 49-54, Aug. 2015, doi: 10.1016/j.jnucmat.2014.11.076.
+
+[^eich_2011]: T. Eich, B. Sieglin, A. Scarabosio, W. Fundamenski, R. J. Goldston, and A. Herrmann, “Inter-ELM Power Decay Length for JET and ASDEX Upgrade: Measurement and Comparison with Heuristic Drift-Based Model,” Physical Review Letters, vol. 107, no. 21, Nov. 2011, doi: https://doi.org/10.1103/PhysRevLett.107.215001

process/core/io/plot/summary.py

@@ -84,6 +84,7 @@
     PlasmaShapeModelType,
 )
 from process.models.physics.profiles import PlasmaProfileShapeType
+from process.models.physics.scrape_off_layer import ScrapeOffLayer
 from process.models.pulse import PulseTimings
 from process.models.superconductors import SuperconductorModel
 from process.models.tfcoil.base import (
@@ -9349,6 +9350,160 @@ def make_bbox_props(power: float) -> dict[str, Any]:
     axis.get_yaxis().set_ticks([])
 
 
+def plot_midplane_near_sol_radial_profile(
+    axis: plt.Axes, mfile: MFile, scan: int, colour_scheme: int
+):
+    """Function to plot the radial profile of the near SOL at the midplane."""
+    rmajor = mfile.get("rmajor", scan=scan)
+    rminor = mfile.get("rminor", scan=scan)
+    len_sol_outboard_power_decay = mfile.get("len_sol_outboard_power_decay", scan=scan)
+    r = np.linspace(
+        (rmajor + rminor), (rmajor + rminor) + (7 * len_sol_outboard_power_decay), 100
+    )
+
+    radial_profile = (
+        ScrapeOffLayer().calculate_outboard_midplane_near_sol_radial_profile(
+            rmajor=rmajor,
+            rminor=rminor,
+            len_plasma_sol_power_decay=len_sol_outboard_power_decay,
+            pflux_plasma_outboard_sol_parallel_mw=mfile.get(
+                "pflux_plasma_outboard_sol_parallel_mw", scan=scan
+            ),
+            r=r,
+        )
+    )
+
+    x_ref = rmajor + rminor + len_sol_outboard_power_decay
+    plasma_boundary = rmajor + rminor
+
+    axis.axvspan(
+        0,
+        plasma_boundary,
+        color=PLASMA_COLOUR[colour_scheme - 1],
+        alpha=0.35,
+        label="Plasma",
+    )
+    axis.plot(r, radial_profile, label=r"$q_{||}$ profile")
+    axis.axvline(
+        x_ref,
+        color="black",
+        linestyle="--",
+        linewidth=1,
+        label=r"$\lambda_{q,\mathrm{out}}$",
+    )
+    axis.grid()
+    axis.legend()
+    axis.set_xlim(
+        (rmajor + rminor - (3 * len_sol_outboard_power_decay)),
+        (rmajor + rminor) + (7 * len_sol_outboard_power_decay),
+    )
+    axis.set_title(r"Midplane Near SOL Radial Profile")
+    axis.minorticks_on()
+    axis.tick_params(axis="x", labelbottom=False)
+    axis.set_ylabel(r"$q_{||}$ [MW/m$^2$]")
+
+
+def plot_div_lower_outboard_eich_target_profile(axis: plt.Axes, mfile: MFile, scan: int):
+    """Function to plot the Eich target profile at the lower outboard divertor."""
+    rmajor = mfile.get("rmajor", scan=scan)
+    rminor = mfile.get("rminor", scan=scan)
+    len_plasma_sol_power_decay = mfile.get(
+        "len_plasma_sol_eich13_power_decay", scan=scan
+    )
+    len_div_outboard_lower_power_spreading = mfile.get(
+        "len_div_outboard_lower_power_spreading", scan=scan
+    )
+    f_b_flux_expansion = mfile.get("f_b_div_outboard_lower_flux_expansion", scan=scan)
+    r = np.linspace(
+        (rmajor + rminor)
+        - ((f_b_flux_expansion / 2) * len_plasma_sol_power_decay) * f_b_flux_expansion,
+        (rmajor + rminor) + (3 * len_plasma_sol_power_decay) * f_b_flux_expansion,
+        200,
+    )
+
+    pflux_target_profile = ScrapeOffLayer().calculate_eich_target_heat_flux_profile(
+        rmajor=rmajor,
+        rminor=rminor,
+        pflux_plasma_sol_parallel_mw=mfile.get(
+            "pflux_plasma_outboard_sol_parallel_mw", scan=scan
+        ),
+        len_plasma_sol_power_decay=mfile.get("len_sol_outboard_power_decay", scan=scan),
+        f_b_div_flux_expansion=f_b_flux_expansion,
+        len_plasma_sol_power_spreading=len_div_outboard_lower_power_spreading,
+        pflux_target_background_heat_flux_mw=0.0,
+        r=r,
+    )
+    peak_idx = np.argmax(pflux_target_profile)
+    peak_r = r[peak_idx]
+    peak_q = pflux_target_profile[peak_idx]
+
+    axis.plot(r, pflux_target_profile)
+    axis.axvline(peak_r, color="black", linestyle="--", linewidth=1)
+    axis.axhline(peak_q, color="black", linestyle="--", linewidth=1)
+    axis.text(
+        0.8,
+        0.9,
+        f"$f_x$ = {f_b_flux_expansion:.2f}\n$S$ = {len_div_outboard_lower_power_spreading * 1e3:.3f} mm",
+        transform=axis.transAxes,
+        ha="left",
+        va="top",
+        bbox={"boxstyle": "round", "facecolor": "white", "alpha": 1.0},
+    )
+    axis.grid()
+    axis.minorticks_on()
+    axis.set_title(r"Lower Outboard Eich Target Parallel Heat Flux Profile")
+    axis.set_xlabel("Radial Position [m]")
+    axis.set_xlim(r[0], r[-1])
+    axis.set_ylabel(r"$q_{||,t}$ [MW/m$^2$]")
+
+
+def plot_sol_power_flux_profiles(axis: plt.Axes, mfile: MFile, scan: int, colour_scheme):
+    """Plot separatrix power split fractions as a bar chart."""
+    plot_plasma(axis=axis, mfile=mfile, scan=scan, colour_scheme=colour_scheme)
+    rmajor, rminor, kappa = mfile.get_variables(
+        "rmajor",
+        "rminor",
+        "kappa",
+        scan=scan,
+    )
+    len_sol_outboard_power_decay = mfile.get("len_sol_outboard_power_decay", scan=scan)
+    a_plasma_outboard_sol_parallel = mfile.get(
+        "a_plasma_outboard_sol_parallel", scan=scan
+    )
+    pflux_plasma_outboard_sol_parallel_mw = mfile.get(
+        "pflux_plasma_outboard_sol_parallel_mw", scan=scan
+    )
+    plasma_scale = max(rminor, abs(kappa * rminor), 1e-6)
+    scale_factor = min(max(plasma_scale / 2.0, 0.7), 1.0)
+    text_fontsize = 9 * scale_factor
+
+    outboard_pos = (rmajor + rminor, 0.0)
+
+    axis.text(
+        *outboard_pos,
+        f"$\\lambda_q = {len_sol_outboard_power_decay * 1e3:.3f}$ mm\n"
+        f"$A_{{||}} = {a_plasma_outboard_sol_parallel:.4f}$ m$^2$\n"
+        f"$q_{{||}} = {pflux_plasma_outboard_sol_parallel_mw:,.2f}$ MW/m$^2$",
+        fontsize=text_fontsize,
+        verticalalignment="center",
+        horizontalalignment="center",
+        bbox={
+            "boxstyle": f"round,pad={0.3 * scale_factor:.3f}",
+            "alpha": 1.0,
+            "linewidth": 2 * scale_factor,
+            "edgecolor": "black",
+        },
+        zorder=101,
+    )
+
+    axis.spines["top"].set_visible(False)
+    axis.spines["right"].set_visible(False)
+    axis.spines["bottom"].set_visible(False)
+    axis.spines["left"].set_visible(False)
+    axis.get_xaxis().set_ticks([])
+    axis.get_yaxis().set_ticks([])
+
+
 def plot_h_threshold_comparison(axis: plt.Axes, mfile: MFile, scan: int, u_seed=None):
     """Function to plot a scatter box plot of L-H threshold power comparisons.
 
@@ -16731,6 +16886,20 @@ def _add_page(name: str | None = None):
         pages["plasma_exhaust"].add_subplot(122), m_file, scan, colour_scheme
     )
 
+    plot_sol_power_flux_profiles(
+        _add_page("sol_powerfluxes").add_subplot(121), m_file, scan, colour_scheme
+    )
+
+    ax_midplane_near_sol = pages["sol_powerfluxes"].add_subplot(336)
+    plot_midplane_near_sol_radial_profile(
+        ax_midplane_near_sol, m_file, scan, colour_scheme
+    )
+
+    ax_div_lower_outboard = pages["sol_powerfluxes"].add_subplot(
+        339, sharex=ax_midplane_near_sol
+    )
+    plot_div_lower_outboard_eich_target_profile(ax_div_lower_outboard, m_file, scan)
+
     plot_debye_length_profile(
         _add_page("microscopic_quantities").add_subplot(232), m_file, scan
     )

process/data_structure/physics_variables.py

@@ -1761,6 +1761,15 @@ class PhysicsData:
     - =3 MAST 2014 scaling 2
     """
 
+    f_b_div_outboard_lower_flux_expansion: float = 5.0
+    """Outboard lower divertor flux expansion factor for the divertor targets (fₓ)"""
+
+    len_div_outboard_lower_scarabosio15_power_spreading: float = 0.0
+    """Scarabosio 2015 H-mode power spreading length/factor in the scrape-off layer scaling (S) [m]"""
+
+    len_div_outboard_lower_power_spreading: float = 0.0
+    """Power spreading length/factor at the outboard lower divertor target (S) [m]"""
+
     dt_power_density_plasma: float = 0.0
     sigmav_dt_average: float = 0.0
     dhe3_power_density: float = 0.0