3258 add divertor and fwblanket lifetimes to st availability model

documentation/proc-pages/eng-models/plant-availability.md

@@ -76,16 +76,24 @@ All availability models in PROCESS require the calculation of the centerpost lif
 
 For superconducting magnets (`i_tf_sup = 1`), the centrepost lifetime is calculated as
 
-$$ t_{\text{CP,life}} = min(f_{\text{TF,max}}/(\phi_{\text{CP,max}}t_{\text{year}}),t_{\text{life}}) $$
+$$ t_{\text{CP,life}} = \min(f_{\text{TF,max}}/(\phi_{\text{CP,max}}t_{\text{year}}),t_{\text{life}}) $$
 
 where $f_{\text{TF,max}}$ is the max fast neutron fluence on the TF coil ($\mathrm{m}^{-2} \mathrm{s}$), $\phi_{\text{CP,max}}$ is the centrepost TF fast neutron flux ($\mathrm{m}^{-2}$ $\mathrm{s}^{-1}$) and $t_{\text{year}}$ is the number of seconds in a year. 
 
 For copper or cryogenic aluminium magnets (`i_tf_sup = 0 or 2`), the centrepost lifetime is
 
-$$ t_{\text{CP,life}} = min(f_{\text{CP, allowable}}/P_{\text{wall}}, t_{\text{life}}) $$
+$$ t_{\text{CP,life}} = \min(f_{\text{CP, allowable}}/P_{\text{wall}}, t_{\text{life}}) $$
 
 where $f_{\text{CP, allowable}}$ is the allowable centrepost neutron fluence and $P_{\text{wall}}$ is the average neutron wall load ($\mathrm{MW} \mathrm{m}^{-2}$).
 
+## Divertor lifetime
+
+The divertor lifetime is calculated as
+
+$$ t_{\text{div, life}} = \max (0, \min(f_{\text{div, allowable}} / P_{\text{div}}, t_{\text{life}})) $$
+
+where $f_{\text{div, allowable}}$ is the allowable divertor heat fluence ($\mathrm{MW}\text{-}\mathrm{yr} \mathrm{m}^{-2}$) and $P_{\text{div}}$ is the heat load to the divertor ($\mathrm{MW} \mathrm{m}^{-2}$).
+
 [^1]: P. J. Knight, *"PROCESS 3020: Plant Availability Model"*, Work File Note
 F/PL/PJK/PROCESS/CODE/<br>
 [^2]: M. Kovari, F. Fox, C. Harrington, R. Kembleton, P. Knight, H. Lux, J. Morris *"PROCESS: a systems code for fusion power plants - Part 2: Engineering"*, Fus. Eng. & Des. 104, 9-20 (2016)

process/availability.py

@@ -77,7 +77,7 @@ def avail(self, output: bool):
 
         # Full power lifetime (in years)
         if ifev.ife != 1:
-            # Caculate DPA per FPY - based on neutronics-derived fusion power relation to DEMO blanket lifetime provided by Matti Coleman
+            # Calculate DPA per FPY - based on neutronics-derived fusion power relation to DEMO blanket lifetime provided by Matti Coleman
             # Detailed and cited in T. Franke 2020, "The EU DEMO equatorial outboard limiter — Design and port integration concept"
             # https://www.sciencedirect.com/science/article/pii/S0920379620301952#bib0075
             # Scaling w.r.t. fusion power drops out a large number of factors relating to neutronics, such as:
@@ -117,7 +117,7 @@ def avail(self, output: bool):
                 dv.hldiv = 1.0e-10
 
             # Divertor lifetime (years)
-            cv.divlife = max(0.0, min(cv.adivflnc / dv.hldiv, cv.tlife))
+            cv.divlife = self.divertor_lifetime()
 
             # Centrepost lifetime (years) (ST machines only)
             if pv.itart == 1:
@@ -477,7 +477,7 @@ def calc_u_planned(self, output: bool) -> float:
             fwbsv.bktlife = min(cv.life_dpa / dpa_fpy, cv.tlife)  # DEMO
 
         # Divertor lifetime (years)
-        cv.divlife = min(cv.adivflnc / dv.hldiv, cv.tlife)
+        cv.divlife = self.divertor_lifetime()
 
         # Centrepost lifetime (years) (ST only)
         if pv.itart == 1:
@@ -1003,17 +1003,65 @@ def calc_u_unplanned_vacuum(self, output: bool) -> float:
         return u_unplanned_vacuum
 
     def avail_st(self, output: bool):
-        """Routine to calculate availability for plant with a Spherical Tokamak
-
-        :param output: indicate whether output should be written to the output file, or not
-        :type output: boolean
         """
-        # CP lifetime
+        Calculate the availability for a plant with a Spherical Tokamak.
+
+        This routine calculates the availability of a plant by considering various factors such as
+        the lifetime of different components, planned and unplanned unavailability, and maintenance cycles.
+
+        Parameters:
+        -----------
+        output : bool
+            Indicates whether the output should be written to the output file or not.
+
+        Detailed Description:
+        ---------------------
+        - The method calculates the Displacements Per Atom (DPA) per Full Power Year (FPY) based on the fusion power.
+        - It determines the lifetime of the first wall and blanket, divertor, and current drive.
+        - It calculates the time for a maintenance cycle and the number of maintenance cycles over the plant's lifetime.
+        - It computes the planned and unplanned unavailability of various components such as magnets, divertor, first wall and blanket, balance of plant, heating and current drive, and vacuum systems.
+        - The total availability of the plant is then calculated considering both planned and unplanned unavailability.
+        - The method also adjusts the lifetimes of components based on the calculated availability.
+        - Finally, it calculates the capacity factor and operational time of the plant.
+
+        If `output` is True, the method writes detailed availability information to the output file.
+
+        References:
+        -----------
+        - T. Franke 2020, "The EU DEMO equatorial outboard limiter — Design and port integration concept"
+          https://www.sciencedirect.com/science/article/pii/S0920379620301952#bib0075
+
+        Notes:
+        ------
+        - The method assumes certain constants and reference points for calculations.
+        - The method modifies the lifetimes of components to account for the calculated availability.
+        """
+
+        ref_powfmw = 2.0e3  # (MW) fusion power for EU-DEMO
+        f_scale = pv.fusion_power / ref_powfmw
+        ref_dpa_fpy = 10.0e0  # dpa per fpy from T. Franke 2020 states up to 10 dpa/FPY
+        dpa_fpy = f_scale * ref_dpa_fpy
+
+        if cv.ibkt_life == 0:
+            fwbsv.bktlife = min(cv.abktflnc / pv.wallmw, cv.tlife)
+        else:
+            fwbsv.bktlife = min(cv.life_dpa / dpa_fpy, cv.tlife)  # DEMO
+
+        # Divertor lifetime (years)
+        cv.divlife = self.divertor_lifetime()
+
+        # CP lifetime (years)
         cv.cplife = self.cp_lifetime()
 
+        # Current drive lifetime (assumed equal to first wall and blanket lifetime)
+        cv.cdrlife = fwbsv.bktlife
+
         # Time for a maintenance cycle (years)
-        # Lifetime of CP + time to replace
-        maint_cycle = cv.cplife + cv.tmain
+        # Shortest component lifetime + time to replace
+        shortest_lifetime = min(
+            fwbsv.bktlife, cv.divlife, cv.cplife, cv.cdrlife, cv.tlife
+        )
+        maint_cycle = shortest_lifetime + cv.tmain
 
         # Number of maintenance cycles over plant lifetime
         n_cycles_main = cv.tlife / maint_cycle
@@ -1071,10 +1119,49 @@ def avail_st(self, output: bool):
             1.0e0 - (u_planned + u_unplanned + u_planned * u_unplanned), 0.0e0
         )
 
+        # Modify lifetimes to take account of the availability
+        if ifev.ife != 1:
+            # First wall / blanket
+            if fwbsv.bktlife < cv.tlife:
+                fwbsv.bktlife = min(fwbsv.bktlife / cv.cfactr, cv.tlife)
+                cv.cdrlife = fwbsv.bktlife
+
+            # Divertor
+            if cv.divlife < cv.tlife:
+                cv.divlife = min(cv.divlife / cv.cfactr, cv.tlife)
+
+            # Centrepost
+            if pv.itart == 1 and cv.cplife < cv.tlife:
+                cv.cplife = min(cv.cplife / cv.cfactr, cv.tlife)
+
         # Capacity factor
         cv.cpfact = cv.cfactr * (tv.t_burn / tv.t_cycle)
 
         if output:
+            po.ocmmnt(self.outfile, "Plant Availability")
+            po.oblnkl(self.outfile)
+            po.ovarre(
+                self.outfile,
+                "Allowable blanket neutron fluence (MW-yr/m2)",
+                "(abktflnc)",
+                cv.abktflnc,
+            )
+            po.ovarre(
+                self.outfile,
+                "Allowable divertor heat fluence (MW-yr/m2)",
+                "(adivflnc)",
+                cv.adivflnc,
+            )
+            po.ovarre(
+                self.outfile,
+                "First wall / blanket lifetime (FPY)",
+                "(bktlife)",
+                fwbsv.bktlife,
+                "OP ",
+            )
+            po.ovarre(
+                self.outfile, "Divertor lifetime (FPY)", "(divlife)", cv.divlife, "OP "
+            )
             if tfv.i_tf_sup == 1:
                 po.ovarre(
                     self.outfile,
@@ -1113,6 +1200,13 @@ def avail_st(self, output: bool):
                 "OP ",
             )
             po.oblnkl(self.outfile)
+            po.ovarre(
+                self.outfile,
+                "Maintenance time for replacing CP (years)",
+                "(tmain)",
+                cv.tmain,
+                "OP ",
+            )
             po.ovarre(
                 self.outfile,
                 "Length of maintenance cycle (years)",
@@ -1174,7 +1268,8 @@ def avail_st(self, output: bool):
                 self.outfile, "Total plant lifetime (years)", "(tlife)", cv.tlife, "OP"
             )
 
-    def cp_lifetime(self):
+    @staticmethod
+    def cp_lifetime():
         """Calculate Centrepost Lifetime
 
         This routine calculates the lifetime of the centrepost,
@@ -1189,8 +1284,20 @@ def cp_lifetime(self):
             cplife = min(ctv.nflutfmax / (fwbsv.neut_flux_cp * YEAR_SECONDS), cv.tlife)
 
         # Aluminium/Copper magnets CP lifetime
-        # For now, we keep the original def, developped for GLIDCOP magnets ...
+        # For now, we keep the original def, developed for GLIDCOP magnets ...
         else:
             cplife = min(cv.cpstflnc / pv.wallmw, cv.tlife)
 
         return cplife
+
+    @staticmethod
+    def divertor_lifetime():
+        """Calculate Divertor Lifetime
+
+        This routine calculates the lifetime of the divertor based on the allowable divertor heat fluence.
+        :returns: Divertor lifetime
+        :rtype: float
+        """
+        # Divertor lifetime
+        # Either 0.0, calculated from allowable divertor fluence and heat load, or lifetime of the plant
+        return max(0.0, min(cv.adivflnc / dv.hldiv, cv.tlife))

tests/unit/test_availability.py

@@ -573,13 +573,20 @@ def test_avail_st(monkeypatch, availability):
     monkeypatch.setattr(cv, "tlife", 30.0)
     monkeypatch.setattr(cv, "u_unplanned_cp", 0.05)
     monkeypatch.setattr(tv, "t_burn", 5.0)
-    monkeypatch.setattr(tv, "t_cycle", 10.0)
+    monkeypatch.setattr(tv, "t_cycle", 9000.0)
+    monkeypatch.setattr(cv, "adivflnc", 10.0)
+    monkeypatch.setattr(dv, "hldiv", 10.0)
+    monkeypatch.setattr(cv, "ibkt_life", 0)
+    monkeypatch.setattr(cv, "abktflnc", 10.0)
+    monkeypatch.setattr(pv, "wallmw", 10.0)
+    monkeypatch.setattr(cv, "cplife", 5.0)
+    monkeypatch.setattr(cv, "cdrlife", 15.0)
 
     availability.avail_st(output=False)
 
-    assert pytest.approx(cv.t_operation) == 29.03225806
-    assert pytest.approx(cv.cfactr) == 0.82579737
-    assert pytest.approx(cv.cpfact) == 0.41289868
+    assert pytest.approx(cv.t_operation) == 15.0
+    assert pytest.approx(cv.cfactr) == 0.27008858
+    assert pytest.approx(cv.cpfact, abs=1.0e-8) == 0.00015005
 
     # Initialise fortran variables again to reset for other tests
     fortran.init_module.init_all_module_vars()
@@ -606,3 +613,23 @@ def test_cp_lifetime(monkeypatch, availability, i_tf_sup, exp):
     cplife = availability.cp_lifetime()
 
     assert pytest.approx(cplife) == exp
+
+
+def test_divertor_lifetime(monkeypatch, availability):
+    """Test divertor_lifetime routine
+
+    :param monkeypatch: Mock fixture
+    :type monkeypatch: object
+
+    :param availability: fixture containing an initialised `Availability` object
+    :type availability: tests.unit.test_availability.availability (functional fixture)
+    """
+
+    monkeypatch.setattr(cv, "adivflnc", 100.0)
+    monkeypatch.setattr(dv, "hldiv", 10.0)
+    monkeypatch.setattr(cv, "tlife", 30.0)
+
+    divlife_obs = availability.divertor_lifetime()
+    divlife_exp = 10.0
+
+    assert pytest.approx(divlife_obs) == divlife_exp