Numerical simulation of coupling heat transfer and thermal stress for spent fuel dry storage cask with different power distribution and tilt angles

Numerical simulation of coupling heat transfer and thermal stress for spent fuel dry storage cask with different power distribution and tilt angles 1 first-author 1 corresp * cjj-njut@njtech.edu.cn * cjj-njut@njtech.edu.cn cjj-njut@njtech.edu.cn https://orcid.org/0000-0002-7579-8024 2 https://orcid.org/0000-0002-8717-4112 1 College of Urban Construction, Nanjing Tech University, Nanjing 211800, China College of Urban Construction, Nanjing Tech University, Nanjing 211800, China School of Energy Science and Engineering, Nanjing Tech University, Nanjing 211800, China School of Energy Science and Engineering, Nanjing Tech University, Nanjing 211800, China Dry storage containers must be secure and reliable during long-term storage, and the effect of decay heat released from the internal spent fuel on the cask has become an important research topic. In this paper, a 3D computational fluid dynamics (CFD) model is presented, and the accuracy of the calculation is verified, with computational errors of less than 6.2%. The thermal stress of the dry storage cask was estimated by coupling it with a transient temperature field. The total power remained constant and adjusting the power ratio of the inner and outer zones had a small effect on the stress results, with a maximum equivalent stress of approximately 5.2 kPa, which occurred at the lower edge of the shell. In the case of tilt, the temperature gradient varied in a wavy distribution, and the wave crest moved from right to left. Altering the tilt angle affects the air distribution in the annular gap, leading to the shell temperature being transformed, with a maximum equivalent stress of 202MPa at the bottom of the shell. However, the equivalent stress in both cases was less than the yield stress (205MPa). Thermal stress CFD simulation Spent nuclear fuel Dry storage cask Introduction 1 With the development of nuclear power, the number of nuclear reactors has continued to increase. According to data from the IAEA, there were 437 operating nuclear power units worldwide by the end of 2021 [ 1 1 ] . The amount of spent fuel, highly radioactive fuel periodically removed from nuclear reactors, is increasing with the operation of reactors. Therefore, potential risk factors must be considered during storage, transportation, and preparation. The rapid development of nuclear power has led to a rapid increase in the cumulative amount of spent fuel, and spent fuel treatment technology is insufficient. Currently, long-term storage measures are available [, 2 2 , 3 3 ] . Long-term storage measures include dry storage and storage in the reactor pool (wet storage) [ 4 4 ] , as shown in Table 1 Table 1 . According to the current development trend, the reprocessing method does not prevent the accumulation of spent fuel. Moreover, the accumulation of a large volume of spent fuel results in a fulfilled capacity of the reactor, making refueling impossible. Therefore, off-reactor dry storage has become the primary treatment method for spent fuel [,, 3 3 , 5 5 , 6 6 ] . 1 1 Treatment method 1 1 Methods 1 1 Advantages 1 1 Limitations 1 1 Long-term storage 1 1 Wet storage 1 1 Storage in the reactor pool 1 1 High capacity, intensive storage, and easy operation 1 1 Storage capacity is limited and difficult to expand 1 1 Dry storage 1 1 Storage in a metal cask with thick concrete walls 1 1 Cheaper to build and maintain than wet storage, easy to transport 1 1 Incomplete recycling component technology 1 1 Reprocessing method 1 1 Chemical treatment 1 1 Improve the utilization capacity of nuclear fuel 1 1 No mature closed fuel cycle technology The design of spent fuel storage containers must fulfill the following requirements to prevent some of the effects of extreme accidents that occur during the transportation and storage of spent fuel:1) wrapping radioactive substances with a protective layer, 2) controlling the external radiation level, 3) avoiding dangerous accidents, and 4) preventing damage caused by heat. Its design and construction requirements require reliable storage of radioactive waste, radiation shielding, subcritical, and heat dissipation under any possible external load [,, 8 8 - 10 10 8-10 ] . During the design and manufacture of containers, possible problems in the accident state are considered, and the structural integrity of the container can be guaranteed under extreme conditions through simulations and performance testing [ 11 11 ] . Generally, the performance test comprises a drop test, impact resistance test, and earthquake resistance test. In addition to safety problems, it is also necessary to consider the impact of the aging of the container components caused by long-term storage environments on fuel recovery [ 12 12 ] , and the impact of the deformation of relevant components on spent fuel storage and recovery in accidents [,, 13 13 - 15 15 13-15 ] . With the rapid development of computers, CFD methods have been used to simulate various thermodynamic problems. Owing to the high radioactivity of spent fuel, performing experiments is difficult. Many studies have shown that the CFD method effectively and accurately predicts thermodynamic phenomena in the nuclear field [ 16 16 ] . Li et al. [ 17 17 ] established a three-dimensional model of a full-scale PWR and used CFD software to simulate the coolant mixing phenomenon. The results showed that reverse flow occurred in the core entrance section, and the mixing effect was enhanced. Wang et al. [ 18 18 ] demonstrated three internal structural schemes of the static pressure chamber and, through 3D CFD simulation, found that the mass flow rate of the central fuel assembly at the core inlet with vortex suppression plate or port skirt structure was greater than that of the peripheral fuel assemblies. Lee established a computational fluid dynamics (CFD) model of a fuel basket containing a fuel assembly; deduced the temperature difference relationship between the fuel basket and the internal components of the cask; and determined the temperature difference between the cylinder cavity and ambient air by using the concept of thermal resistance, which was verified by CFD analysis [ 19 19 ] . Based on the natural convection heat transfer model, Yoo H S selected the candidate and added a gas group that improved the heat transfer in dry cask storage systems through figure-of-merit analysis and evaluated the candidate gas through CFD modeling without the assumption of porous media [ 20 20 ] . Rodriguez et al. first simulated the inner surface of the container cavity to the surrounding air and then executed a full simulation of the cladding of the internal fuel assembly to the outer surface of the canister. The former served as the actual boundary condition of the latter and estimated the influence of the container on the nearby environment [ 21 21 ] . Poškas used the numerical modeling code ALGOR to model and analyze a newly loaded fuel cask (pre-stored in pools for 5 years) and a temporary fuel cask after 50 years of storage in open storage facilities under natural conditions. A local sensitivity analysis was also conducted for the parameters that primarily affect temperature distribution [ 22 22 ] . Lee et al. conducted a scale analysis to derive the scale ratio between the prototype fuel cask and the miniature model. The thermal analysis and test results of the prototype fuel cask and the semi-scale model were in satisfactory agreement. The thermal test results confirmed the thermal fluid flow similarity between the prototype cask and semi-scale model [ 23 23 ] . In summary, the literature primarily uses numerical simulations to investigate the temperature field. The spent fuel storage cask is generally filled with inert gas to improve heat transfer. The primary heat transfer methods are convection and radiation heat transfer, such that the outer wall temperature is lower than the internal component temperature [, 24 24 , 25 25 ] , resulting in uneven temperature distribution. Owing to the uneven temperature difference when heated, the inconsistent expansion and contraction of the container or fuel assembly induce thermal stress, which affects the long-term storage security of the cask. In addition, the spent fuel storage cask may tilt owing to the impact of storage location, which has potential safety hazards, and there are few relevant studies. Based on the research on temperature distribution, this study investigates the influence of the temperature field on the stress field and storage cask under different power distributions and inclination angles by using the unidirectional coupling numerical simulation of the temperature and stress fields of the dry storage cask. Model and method 1 Geometric Model 2 The dimensions of the cask model and concrete shell were based on Holtec's MPC-32 [ 26 26 ] . We reduced the inner structure of the storage cask to a single basket, and the spent fuel assembly was simplified to helium through a porous media model. According to the reference, the porosity of the solid component is 0.652 (solid volume/total volume=0.348) and approximately equal to the estimated value, whereas the correlation coefficient in the literature is obtained from experiments. The porosity in this study was 0.348, the viscosity coefficient was 687973 m -2 , and the inertia coefficient was 29.86 m -1 [ 26 26 ] . The overall model includes concrete cladding, an interior annular gap, a metal shell, and a basket ( Fig.1 Fig.1 ). The inclination angle was tilted θ (15°, 45°, 75°, and 90°) based on the central axis of symmetry ( Fig.1 Fig.1 ). Table 2 Table 2 presents the specific structural dimensions of the models. Because the overall structure is center-symmetric, a 1/8 model was used to spare computing resources. 1 1 Parameters 1 1 Value 1 1 Metal cask 1 1 External diameter (m) 1 1 1.806 1 1 Internal diameter (m) 1 1 1.780 1 1 Height (m) 1 1 4.850 1 1 Upper lid (m) 1 1 0.290 1 1 Thickness (m) 1 1 0.013 1 1 Basket thickness (m) 1 1 0.0096 1 1 Basket cell (m 3 ) 1 1 0.22 × 0.22 × 4.44 1 1 Concrete cladding 1 1 Height (m) 1 1 6.090 1 1 External diameter (m) 1 1 3.360 1 1 Internal diameter (m) 1 1 1.960 1 1 Thickness (m) 1 1 0.700 1 1 Upper lid 1 1 Height (m) 1 1 0.300 1 1 Diameter (m) 1 1 1.806 1 1 Inlet 1 1 Height (m) 1 1 0.250 1 1 Width (m) 1 1 0.380 1 1 Outlet 1 1 Height (m) 1 1 0.150 1 1 Width (m) 1 1 0.642 Physical model 2 Gravity was considered in the z-axis direction. The acceleration of gravity is 9.8 m/s 2 , the external air was 300 K and 101,325 Pa, and the internal helium was 294 K and 330,000 Pa. We made the following assumptions to simplify the calculation of the fluid working conditions: (1) there is transient heat transfer of the 3D model; (2) the fluid was steady, incompressible, and isotropic; (3) all solid materials are isotropic; and (4) the displacement at the liquid-solid coupling is the same. According to the four aforementioned assumptions, the following governing equation can be obtained: Continuous equation <math display="block" id="M1"><mrow><mfrac><mrow><mo>∂</mo><mi>u</mi></mrow><mrow><mo>∂</mo><mi>x</mi></mrow></mfrac><mo>+</mo><mfrac><mrow><mo>∂</mo><mi>v</mi></mrow><mrow><mo>∂</mo><mi>y</mi></mrow></mfrac><mo>+</mo><mfrac><mrow><mo>∂</mo><mi>w</mi></mrow><mrow><mo>∂</mo><mi>z</mi></mrow></mfrac><mo>=</mo><mn>0.</mn></mrow></math> Momentum conservation equation <math display="block" id="M2"><mrow><mi>u</mi><mfrac><mrow><mo>∂</mo><mi>u</mi></mrow><mrow><mo>∂</mo><mi>x</mi></mrow></mfrac><mo>+</mo><mi>v</mi><mfrac><mrow><mo>∂</mo><mi>v</mi></mrow><mrow><mo>∂</mo><mi>y</mi></mrow></mfrac><mo>+</mo><mi>w</mi><mfrac><mrow><mo>∂</mo><mi>w</mi></mrow><mrow><mo>∂</mo><mi>z</mi></mrow></mfrac><mo>=</mo><mo>−</mo><mfrac><mn>1</mn><mtext>ρ</mtext></mfrac><mfrac><mrow><mo>∂</mo><mi>p</mi></mrow><mrow><mo>∂</mo><mi>x</mi></mrow></mfrac><mo>+</mo><mi>v</mi><mrow><mo>(</mo><mrow><mfrac><mrow><msup><mo>∂</mo><mn>2</mn></msup><mi>u</mi></mrow><mrow><mo>∂</mo><msup><mi>x</mi><mn>2</mn></msup></mrow></mfrac><mo>+</mo><mfrac><mrow><msup><mo>∂</mo><mn>2</mn></msup><mi>v</mi></mrow><mrow><mo>∂</mo><msup><mi>y</mi><mn>2</mn></msup></mrow></mfrac><mo>+</mo><mfrac><mrow><msup><mo>∂</mo><mn>2</mn></msup><mi>w</mi></mrow><mrow><mo>∂</mo><msup><mi>z</mi><mn>2</mn></msup></mrow></mfrac></mrow><mo>)</mo></mrow><mo>,</mo></mrow></math> <math display="block" id="M3"><mrow><mi>u</mi><mfrac><mrow><mo>∂</mo><mi>u</mi></mrow><mrow><mo>∂</mo><mi>x</mi></mrow></mfrac><mo>+</mo><mi>v</mi><mfrac><mrow><mo>∂</mo><mi>v</mi></mrow><mrow><mo>∂</mo><mi>y</mi></mrow></mfrac><mo>+</mo><mi>w</mi><mfrac><mrow><mo>∂</mo><mi>w</mi></mrow><mrow><mo>∂</mo><mi>z</mi></mrow></mfrac><mo>=</mo><mo>−</mo><mfrac><mn>1</mn><mtext>ρ</mtext></mfrac><mfrac><mrow><mo>∂</mo><mi>p</mi></mrow><mrow><mo>∂</mo><mi>y</mi></mrow></mfrac><mo>+</mo><mi>v</mi><mrow><mo>(</mo><mrow><mfrac><mrow><msup><mo>∂</mo><mn>2</mn></msup><mi>u</mi></mrow><mrow><mo>∂</mo><msup><mi>x</mi><mn>2</mn></msup></mrow></mfrac><mo>+</mo><mfrac><mrow><msup><mo>∂</mo><mn>2</mn></msup><mi>v</mi></mrow><mrow><mo>∂</mo><msup><mi>y</mi><mn>2</mn></msup></mrow></mfrac><mo>+</mo><mfrac><mrow><msup><mo>∂</mo><mn>2</mn></msup><mi>w</mi></mrow><mrow><mo>∂</mo><msup><mi>z</mi><mn>2</mn></msup></mrow></mfrac></mrow><mo>)</mo></mrow><mo>,</mo></mrow></math> <math display="block" id="M4"><mtable columnalign="left"><mtr><mtd><mi>u</mi><mfrac><mrow><mo>∂</mo><mi>u</mi></mrow><mrow><mo>∂</mo><mi>x</mi></mrow></mfrac><mo>+</mo><mi>v</mi><mfrac><mrow><mo>∂</mo><mi>v</mi></mrow><mrow><mo>∂</mo><mi>y</mi></mrow></mfrac><mo>+</mo><mi>w</mi><mfrac><mrow><mo>∂</mo><mi>w</mi></mrow><mrow><mo>∂</mo><mi>z</mi></mrow></mfrac><mo>=</mo><mo>−</mo><mfrac><mn>1</mn><mtext>ρ</mtext></mfrac><mfrac><mrow><mo>∂</mo><mi>p</mi></mrow><mrow><mo>∂</mo><mi>z</mi></mrow></mfrac></mtd></mtr><mtr><mtd><mtext>    </mtext><mo>+</mo><mi>v</mi><mrow><mo>(</mo><mrow><mfrac><mrow><msup><mo>∂</mo><mn>2</mn></msup><mi>u</mi></mrow><mrow><mo>∂</mo><msup><mi>x</mi><mn>2</mn></msup></mrow></mfrac><mo>+</mo><mfrac><mrow><msup><mo>∂</mo><mn>2</mn></msup><mi>v</mi></mrow><mrow><mo>∂</mo><msup><mi>y</mi><mn>2</mn></msup></mrow></mfrac><mo>+</mo><mfrac><mrow><msup><mo>∂</mo><mn>2</mn></msup><mi>w</mi></mrow><mrow><mo>∂</mo><msup><mi>z</mi><mn>2</mn></msup></mrow></mfrac></mrow><mo>)</mo></mrow><mo>−</mo><mtext>gρ</mtext><mtext>.</mtext></mtd></mtr></mtable></math> Energy conservation equation <math display="block" id="M5"><mrow><msub><mrow><mi>ρ</mi><mi>c</mi></mrow><mtext>p</mtext></msub><mfrac><mrow><mo>∂</mo><mi>T</mi></mrow><mrow><mo>∂</mo><mi>τ</mi></mrow></mfrac><mo>=</mo><mi>λ</mi><mrow><mo>(</mo><mrow><mfrac><mrow><msup><mo>∂</mo><mn>2</mn></msup><mi>T</mi></mrow><mrow><mo>∂</mo><msup><mi>x</mi><mn>2</mn></msup></mrow></mfrac><mo>+</mo><mfrac><mrow><msup><mo>∂</mo><mn>2</mn></msup><mi>T</mi></mrow><mrow><mo>∂</mo><msup><mi>y</mi><mn>2</mn></msup></mrow></mfrac><mo>+</mo><mfrac><mrow><msup><mo>∂</mo><mn>2</mn></msup><mi>T</mi></mrow><mrow><mo>∂</mo><msup><mi>z</mi><mn>2</mn></msup></mrow></mfrac></mrow><mo>)</mo></mrow><mo>+</mo><msub><mi>q</mi><mi>v</mi></msub><mo>,</mo></mrow></math> where <math id="M6"><mrow><mi>u</mi><mo>,</mo><mi>v</mi></mrow></math> , and <math id="M7"><mi>w</mi></math> are the velocity components in the X, Y, and Z directions, respectively, <math id="M8"><mrow><mtext>m</mtext><mo>/</mo><mtext>s</mtext></mrow></math> ; <math id="M9"><mrow><mtext>g </mtext></mrow></math> is the acceleration of gravity, <math id="M10"><mrow><mtext>m</mtext><mo>/</mo><msup><mtext>s</mtext><mn>2</mn></msup></mrow></math> ; <math id="M11"><mrow><msub><mtext>c</mtext><mtext>p</mtext></msub><mtext> </mtext></mrow></math> is the specific volume at constant pressure, J/ (kg ·K); <math id="M12"><mtext>ρ</mtext></math> is the density, <math id="M13"><mrow><mtext>kg</mtext><mo>/</mo><msup><mtext>m</mtext><mn>3</mn></msup></mrow></math> ; <math id="M14"><mi>T</mi></math> is the temperature function, K; <math id="M15"><mi>τ</mi></math> is the time, <math id="M16"><mtext>s</mtext></math> ; and <math id="M17"><mrow><msub><mi>q</mi><mi>v</mi></msub></mrow></math> is the heat source, <math id="M18"><mtext>W</mtext></math> . Fluid-structure coupling equation: <math display="block" id="M19"><mrow><msub><mi>d</mi><mtext>f</mtext></msub><mo>=</mo><msub><mi>d</mi><mtext>s</mtext></msub><mo>,</mo></mrow></math> <math display="block" id="M20"><mrow><msub><mi>T</mi><mtext>f</mtext></msub><mo>=</mo><msub><mi>T</mi><mtext>s</mtext></msub><mo>,</mo></mrow></math> <math display="block" id="M21"><mrow><msub><mi>q</mi><mtext>f</mtext></msub><mo>=</mo><msub><mi>q</mi><mtext>s</mtext></msub><mo>,</mo></mrow></math> <math display="block" id="M22"><mrow><msub><mi>n</mi><mtext>f</mtext></msub><mo>·</mo><msub><mi>σ</mi><mtext>f</mtext></msub><mo>=</mo><msub><mi>n</mi><mtext>s</mtext></msub><mo>·</mo><msub><mi>σ</mi><mtext>s</mtext></msub><mo>,</mo></mrow></math> where f represents fluid; s represents solid; d is displacement; q is heat; and <math id="M23"><mrow><msub><mi>σ</mi><mi>f</mi></msub><mo> </mo><mtext>and</mtext><mo> </mo><msub><mi>σ</mi><mi>s</mi></msub></mrow></math> are the fluid and solid stresses, respectively. Numerical calculations require an appropriate turbulence model to reflect the actual situation. In this study, the air calculation domain was turbulent, and the helium gas inside was a natural convection. Hence, the k-omega turbulence model was selected from the turbulence model provided by Fluent. The SST-k-omega model considers the turbulent shear stress based on the general k-omega model; thus, it is more accurate than the general model [ 16 16 ] . Therefore, the SST-k-omega turbulence model was adopted in the simulation in this study, and the equation for this model is as follows: <math display="block" id="M24"><mrow><mfrac><mrow><mo>∂</mo><mrow><mo>(</mo><mrow><mtext>ρ</mtext><mi>k</mi></mrow><mo>)</mo></mrow></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mo>+</mo><mfrac><mrow><mo>∂</mo><mrow><mo>(</mo><mrow><mtext>ρ</mtext><mi>k</mi><msub><mi>u</mi><mi>i</mi></msub></mrow><mo>)</mo></mrow></mrow><mrow><mo>∂</mo><msub><mi>x</mi><mi>i</mi></msub></mrow></mfrac><mo>=</mo><mfrac><mrow><mo>∂</mo><mo> </mo></mrow><mrow><mo>∂</mo><msub><mi>x</mi><mi>j</mi></msub></mrow></mfrac><mrow><mo>(</mo><mrow><msub><mi>Γ</mi><mi>k</mi></msub><mfrac><mrow><mo>∂</mo><mi>k</mi></mrow><mrow><mo>∂</mo><msub><mi>x</mi><mi>j</mi></msub></mrow></mfrac></mrow><mo>)</mo></mrow><mo>+</mo><msub><mi>G</mi><mi>k</mi></msub><mo>−</mo><msub><mi>Y</mi><mi>k</mi></msub><mo>+</mo><msub><mi>S</mi><mi>k</mi></msub><mo>,</mo></mrow></math> <math display="block" id="M25"><mrow><mfrac><mrow><mo>∂</mo><mrow><mo>(</mo><mrow><mtext>ρ</mtext><mi>ω</mi></mrow><mo>)</mo></mrow></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mo>+</mo><mfrac><mrow><mo>∂</mo><mrow><mo>(</mo><mrow><mtext>ρ</mtext><mi>ω</mi><msub><mi>u</mi><mi>i</mi></msub></mrow><mo>)</mo></mrow></mrow><mrow><mo>∂</mo><msub><mi>x</mi><mi>i</mi></msub></mrow></mfrac><mo>=</mo><mfrac><mrow><mo>∂</mo><mo> </mo></mrow><mrow><mo>∂</mo><msub><mi>x</mi><mi>j</mi></msub></mrow></mfrac><mrow><mo>(</mo><mrow><msub><mi>Γ</mi><mi>ω</mi></msub><mfrac><mrow><mo>∂</mo><mi>ω</mi></mrow><mrow><mo>∂</mo><msub><mi>x</mi><mi>j</mi></msub></mrow></mfrac></mrow><mo>)</mo></mrow><mo>+</mo><msub><mi>G</mi><mi>ω</mi></msub><mo>−</mo><msub><mi>Y</mi><mi>ω</mi></msub><mo>+</mo><msub><mi>S</mi><mi>ω</mi></msub><mo>.</mo></mrow></math> In contrast with standard k-omega , turbulent viscosity <math id="M26"><mrow><msub><mi>μ</mi><mtext>t</mtext></msub></mrow></math> was calculated as follows: <math display="block" id="M27"><mrow><msub><mi>μ</mi><mi>t</mi></msub><mo>=</mo><mfrac><mrow><mtext>ρ</mtext><mi>k</mi></mrow><mi>ω</mi></mfrac><mfrac><mn>1</mn><mrow><mtext>max</mtext><mrow><mo>[</mo><mrow><mfrac><mn>1</mn><mrow><msup><mi>α</mi><mo>*</mo></msup></mrow></mfrac><mo>,</mo><mfrac><mrow><mi>S</mi><msub><mi>F</mi><mn>2</mn></msub></mrow><mrow><msub><mi>α</mi><mn>1</mn></msub><mi>ω</mi></mrow></mfrac></mrow><mo>]</mo></mrow></mrow></mfrac><mo>,</mo></mrow></math> <math display="block" id="M28"><mrow><msub><mi>F</mi><mn>2</mn></msub><mo>=</mo><mi>tanh</mi><msubsup><mo>∅</mo><mn>2</mn><mn>2</mn></msubsup><mo>,</mo></mrow></math> <math display="block" id="M29"><mrow><msub><mo>∅</mo><mn>2</mn></msub><mo>=</mo><mtext>max</mtext><mrow><mo>[</mo><mrow><mn>2</mn><mfrac><mrow><mo>√</mo><mi>k</mi></mrow><mrow><mn>0.09</mn><mi>ω</mi><mi>y</mi></mrow></mfrac><mo>,</mo><mfrac><mrow><mn>500</mn><mi>μ</mi></mrow><mrow><mtext>ρ</mtext><msup><mi>y</mi><mn>2</mn></msup><mi>ω</mi></mrow></mfrac></mrow><mo>]</mo></mrow><mo>,</mo></mrow></math> where <math id="M30"><mrow><msub><mi>G</mi><mi>k</mi></msub></mrow></math> represents the turbulent kinetic energy generated by the average velocity gradient, and <math id="M31"><mrow><msub><mi>Γ</mi><mi>k</mi></msub></mrow></math> and <math id="M32"><mrow><msub><mi>Γ</mi><mi>ω</mi></msub></mrow></math> represent the effective diffusion coefficients of k and ω, respectively. <math id="M33"><mrow><msub><mi>Y</mi><mi>k</mi></msub></mrow></math> and <math id="M34"><mrow><msub><mi>Y</mi><mi>ω</mi></msub></mrow></math> represent the dissipations of k and ω under the action of turbulence, respectively. <math id="M35"><mrow><msub><mi>S</mi><mi>k</mi></msub></mrow></math> and <math id="M36"><mrow><msub><mi>S</mi><mi>ω</mi></msub></mrow></math> are user-defined source items (zero in this model). <math id="M37"><mrow><msub><mi>σ</mi><mi>k</mi></msub></mrow></math> and <math id="M38"><mrow><msub><mi>σ</mi><mi>ω</mi></msub></mrow></math> are the turbulent Prandtl numbers of k and ω, respectively. <math id="M39"><mi>y</mi></math> is the distance to the adjacent surface, <math id="M40"><mrow><msub><mi>α</mi><mn>1</mn></msub></mrow></math> = 0.31. For spent fuel casks, the temperature distribution is principally affected by convection and radiant heat, with radiation being the main heat transfer mode for spent fuel assemblies. Commonly used radiation models include the Surface to Surface (S2S) model, Discrete Ordinates (DO) model, and ray tracing model. Among these radiation heat transfer models, the calculation accuracy of the S2S model is not as good as that of the DO or ray tracing models. To save computing resources, we did not adopt the ray tracing model. Additionally, the DO model is more appropriate than others when the geometric model used is comparatively complex. The transient calculation time step calculated by the formula (∆ t = <math id="M41"><mrow><mn>1</mn><mo>/</mo><mn>3</mn><mo>·</mo><mi>L</mi><mo>/</mo><mi>V</mi></mrow></math> , L : Characteristic size, V : Volume) was approximately 0.025 s, the convergence speed of the heat transfer calculation was slow, and the accuracy was high, but the calculation results were divergent when the same time step was used in the structure set. Finally, when the step size was 1 s, the heat transfer process converged faster than before. The resulting error was within an acceptable range, and the calculation convergence of the structural module was guaranteed. Boundary conditions and meshing 2 The SIMPLE algorithm was adopted based on the pressure-velocity coupling, and the energy equation was discretized using second-order upwind. The outlet was set as the pressure outlet, according to the calculation of the natural convection Reynolds number on the external metal shell, which was approximately 6000, and the inlet speed, which was approximately 0.46 m/s. The emissivity of the metal wall was 0.85, and the model calculation grid is shown in Fig.2 Fig.2 . Because of limited computing resources, 6.7 million grids were used to save computing time; the temperature is similar to 8.76 million grids and 10 million grids. The minimum orthogonal mass and skewness were used to test the element quality. The skewness of the grid element was ≤0.75, which is acceptable in the calculation. In this study, thermal stress was simulated by varying the power and tilt angle of the model. The total power ( <math id="M42"><mi>Q</mi></math> ) remained unchanged. The spent fuel grid model is divided into inner and outer zones ( Fig.3 Fig.3 ). The dark area is the inner zone, the remainder is the outer zone, the inner zone power is <math id="M43"><mrow><msub><mi>q</mi><mrow><mtext>in</mtext></mrow></msub></mrow></math> , and the outer zone power is <math id="M44"><mrow><msub><mi>q</mi><mrow><mtext>out</mtext></mrow></msub></mrow></math> . The <math id="M45"><mrow><mfrac><mrow><msub><mi>q</mi><mrow><mtext>in</mtext></mrow></msub></mrow><mrow><msub><mi>q</mi><mrow><mtext>out</mtext></mrow></msub></mrow></mfrac></mrow></math> ratio was changed ( Table 3 Table 3 ) to estimate thermal stress. 1 1 <math id="M46"><mrow><mfrac bevelled="true"><mrow><msub><mi>q</mi><mrow><mtext>in</mtext></mrow></msub></mrow><mrow><msub><mi>q</mi><mrow><mtext>out</mtext></mrow></msub></mrow></mfrac></mrow></math> 1 1 <math id="M47"><mrow><msub><mi>q</mi><mrow><mtext>in</mtext></mrow></msub><mrow><mo>(</mo><mrow><mtext>kW</mtext><mo>/</mo><mtext>cell</mtext></mrow><mo>)</mo></mrow></mrow></math> 1 1 <math id="M48"><mrow><msub><mi>q</mi><mrow><mtext>out</mtext></mrow></msub><mtext> </mtext><mrow><mo>(</mo><mrow><mtext>kW</mtext><mo>/</mo><mtext>cell</mtext></mrow><mo>)</mo></mrow></mrow></math> 1 1 <math id="M49"><mrow><mi>Q</mi><mo> </mo><mo stretchy="false">(</mo><mtext>kW</mtext></mrow></math> ) 1 1 0.5 1 1 0.58 1 1 1.15 1 1 30.0 1 1 1 1 1 0.94 1 1 0.94 1 1 30.0 1 1 2 1 1 1.36 1 1 0.68 1 1 30.0 In the mechanical calculation, we calculated the shell and basket, suppressed the other models, and set fixed supports on the upper and lower parts of the basket. In this study, the calculated thermal loads of the shell and basket were introduced into the model in a unidirectional coupling manner, and the calculation step was consistent with the heat transfer process. Model verification 2 The size of the model was analogous to that of Herranz et al. [ 26 26 ] , and the accuracy of the computation was verified by comparison with its temperature field. We compared the temperature distribution from the center to the outside of the model at a height of approximately 4.8 m with that in the literature and discovered that the temperature difference in the basket zone was 0.5%-6.2%, and the maximum temperature difference was 0.8%, at the height of 4.6-4.8 m. Although the materials used in the literature are undefined, the error in the metal container is within an acceptable range. In Fig.4 Fig.4 , in the annular gap outside the container, the temperature in the literature is lower than the simulated temperature because only the base of the concrete layer is set as the adiabatic condition in the simulation process (consistent with the literature); however, the computation results illustrate that the temperature of the outer wall of the container is close to the temperature of the concrete, which is approximately 550 K. This paper principally discusses the temperature and stress changes in the metal cask; therefore, the influence of the external concrete layer temperature on the subsequent calculation in the metal container can be disregarded. Figure 5 Figure 5 shows that the temperature in the metal container after steady-state computation is concentrated at the middle part, and the temperature in the inner zone of the component is higher than that in the outer zone, with the highest temperature of 621.81 K on assembly 14. The maximum temperature of the external concrete layer is approximately 550 K. We compared Fig.5 Fig.5 and Fig.6 Fig.6 and found that the base has a slight change in the external heat dissipation process of the container. All concrete layers should be set as adiabatic conditions to be consistent with the reference. In addition, Fig.6 Fig.6 reveals that in the case of partial convergence, the internal high temperature is mainly distributed in the middle of the component. With the natural convection of helium in the container and heat conduction of the assemblies, Fig.5 Fig.5 displays that the high-temperature part is transferred upward. The model is accurate based on the verification. Thermal stress 2 This study focuses on the stress changes caused by the temperature field, and the stress problem on the basket can be regarded as an axial plane stress problem, which is assumed to fulfill the four basic assumptions of elasticity [ 27 27 ] . This is: <math display="block" id="M50"><mrow><mrow><mo>{</mo><mrow><mtable columnalign="left"><mtr columnalign="left"><mtd columnalign="left"><mrow><msub><mi>σ</mi><mi>y</mi></msub><mo>=</mo><msub><mi>τ</mi><mrow><mi>x</mi><mi>y</mi></mrow></msub><mo>=</mo><msub><mi>τ</mi><mrow><mi>z</mi><mi>y</mi></mrow></msub><mo>=</mo><mn>0</mn></mrow></mtd></mtr><mtr columnalign="left"><mtd columnalign="left"><mrow><msub><mi>σ</mi><mi>x</mi></msub><mo>=</mo><msub><mi>σ</mi><mi>x</mi></msub><mrow><mo>(</mo><mrow><mi>x</mi><mo>,</mo><mi>z</mi></mrow><mo>)</mo></mrow></mrow></mtd></mtr><mtr columnalign="left"><mtd columnalign="left"><mrow><msub><mi>σ</mi><mi>z</mi></msub><mo>=</mo><msub><mi>σ</mi><mi>z</mi></msub><mrow><mo>(</mo><mrow><mi>x</mi><mo>,</mo><mi>z</mi></mrow><mo>)</mo></mrow></mrow></mtd></mtr><mtr columnalign="left"><mtd columnalign="left"><mrow><msub><mi>τ</mi><mrow><mi>x</mi><mi>z</mi></mrow></msub><mo>=</mo><msub><mi>τ</mi><mrow><mi>x</mi><mi>z</mi></mrow></msub><mrow><mo>(</mo><mrow><mi>x</mi><mo>,</mo><mi>z</mi></mrow><mo>)</mo></mrow></mrow></mtd></mtr></mtable></mrow></mrow></mrow></math> where <math id="M51"><mrow><msub><mi>σ</mi><mi>x</mi></msub><mo>,</mo><msub><mi>σ</mi><mi>y</mi></msub><mo>,</mo><msub><mi>σ</mi><mi>z</mi></msub><mo>,</mo><msub><mi>τ</mi><mrow><mi>x</mi><mi>y</mi></mrow></msub><mo>,</mo><msub><mi>τ</mi><mrow><mi>z</mi><mi>y</mi></mrow></msub><mo>,</mo><mtext>and </mtext><msub><mi>τ</mi><mrow><mi>x</mi><mi>z</mi></mrow></msub></mrow></math> are stress components. Consequently, the general physical equation of thermal stress can be simplified as <math display="block" id="M52"><mrow><mrow><mo>{</mo><mrow><mtable columnalign="left"><mtr columnalign="left"><mtd columnalign="left"><mrow><msub><mi>ε</mi><mi>x</mi></msub><mo>=</mo><mfrac><mn>1</mn><mtext>E</mtext></mfrac><mrow><mo>(</mo><mrow><msub><mi>σ</mi><mi>x</mi></msub><mo>−</mo><mtext>ν</mtext><msub><mi>σ</mi><mi>z</mi></msub></mrow><mo>)</mo></mrow><mo>+</mo><mi>α</mi><mi>T</mi></mrow></mtd></mtr><mtr columnalign="left"><mtd columnalign="left"><mrow><msub><mi>ε</mi><mi>z</mi></msub><mo>=</mo><mfrac><mn>1</mn><mtext>E</mtext></mfrac><mrow><mo>(</mo><mrow><msub><mi>σ</mi><mi>z</mi></msub><mo>−</mo><mtext>ν</mtext><msub><mi>σ</mi><mi>x</mi></msub></mrow><mo>)</mo></mrow><mo>+</mo><mi>α</mi><mi>T</mi></mrow></mtd></mtr><mtr columnalign="left"><mtd columnalign="left"><mrow><msub><mi>γ</mi><mrow><mi>x</mi><mi>z</mi></mrow></msub><mo>=</mo><mfrac><mrow><mn>2</mn><mrow><mo>(</mo><mrow><mn>1</mn><mo>+</mo><mtext>ν</mtext></mrow><mo>)</mo></mrow></mrow><mtext>E</mtext></mfrac><msub><mi>τ</mi><mrow><mi>x</mi><mi>z</mi></mrow></msub></mrow></mtd></mtr></mtable></mrow></mrow></mrow></math> where <math id="M53"><mrow><msub><mi>ε</mi><mi>x</mi></msub><mo>,</mo><msub><mi>ε</mi><mi>z</mi></msub><mo>,</mo><mtext>and </mtext><msub><mi>γ</mi><mrow><mi>x</mi><mi>z</mi></mrow></msub></mrow></math> are the strain components, <math id="M54"><mtext>ν</mtext></math> is Poisson's coefficient, <math id="M55"><mi>α</mi></math> is the thermal expansion coefficient, <math id="M56"><mtext>E</mtext></math> is the elastic modulus, and <math id="M57"><mi>T</mi></math> is the temperature variable. Substituting the plane equation into equation (16) yields <math display="block" id="M58"><mrow><mrow><mo>{</mo><mrow><mtable columnalign="left"><mtr columnalign="left"><mtd columnalign="left"><mrow><msub><mi>σ</mi><mi>x</mi></msub><mo>=</mo><mfrac><mtext>E</mtext><mrow><mn>1</mn><mo>−</mo><msup><mtext>ν</mtext><mn>2</mn></msup></mrow></mfrac><mrow><mo>(</mo><mrow><mfrac><mrow><mo>∂</mo><mi>i</mi></mrow><mrow><mo>∂</mo><mi>x</mi></mrow></mfrac><mo>+</mo><mtext>ν</mtext><mfrac><mrow><mo>∂</mo><mi>k</mi></mrow><mrow><mo>∂</mo><mi>z</mi></mrow></mfrac></mrow><mo>)</mo></mrow><mo>−</mo><mfrac><mrow><mtext>E</mtext><mi>α</mi></mrow><mrow><mn>1</mn><mo>−</mo><mtext>ν</mtext></mrow></mfrac><mi>T</mi></mrow></mtd></mtr><mtr columnalign="left"><mtd columnalign="left"><mrow><msub><mi>σ</mi><mi>z</mi></msub><mo>=</mo><mfrac><mtext>E</mtext><mrow><mn>1</mn><mo>−</mo><msup><mtext>ν</mtext><mn>2</mn></msup></mrow></mfrac><mrow><mo>(</mo><mrow><mfrac><mrow><mo>∂</mo><mi>k</mi></mrow><mrow><mo>∂</mo><mi>z</mi></mrow></mfrac><mo>+</mo><mtext>ν</mtext><mfrac><mrow><mo>∂</mo><mi>i</mi></mrow><mrow><mo>∂</mo><mi>x</mi></mrow></mfrac></mrow><mo>)</mo></mrow><mo>−</mo><mfrac><mrow><mtext>E</mtext><mi>α</mi></mrow><mrow><mn>1</mn><mo>−</mo><mtext>ν</mtext></mrow></mfrac><mi>T</mi></mrow></mtd></mtr><mtr columnalign="left"><mtd columnalign="left"><mrow><msub><mi>τ</mi><mrow><mi>x</mi><mi>z</mi></mrow></msub><mo>=</mo><mfrac><mi>E</mi><mrow><mn>2</mn><mrow><mo>(</mo><mrow><mn>1</mn><mo>+</mo><mtext>ν</mtext></mrow><mo>)</mo></mrow></mrow></mfrac><mrow><mo>(</mo><mrow><mfrac><mrow><mo>∂</mo><mi>i</mi></mrow><mrow><mo>∂</mo><mi>z</mi></mrow></mfrac><mo>+</mo><mfrac><mrow><mo>∂</mo><mi>k</mi></mrow><mrow><mo>∂</mo><mi>x</mi></mrow></mfrac></mrow><mo>)</mo></mrow></mrow></mtd></mtr></mtable></mrow></mrow></mrow></math> where <math id="M59"><mi>i</mi></math> and <math id="M60"><mi>k</mi></math> represent the displacements in the X -and Y directions. Substituting equation (17) into the equilibrium equation yields <math display="block" id="M61"><mrow><mrow><mo>{</mo><mrow><mtable><mtr><mtd><mrow><mfrac><mrow><msup><mo>∂</mo><mn>2</mn></msup><mi>i</mi></mrow><mrow><mo>∂</mo><msup><mi>x</mi><mn>2</mn></msup></mrow></mfrac><mo>+</mo><mfrac><mrow><mn>1</mn><mo>−</mo><mtext>ν</mtext></mrow><mn>2</mn></mfrac><mfrac><mrow><msup><mo>∂</mo><mn>2</mn></msup><mi>i</mi></mrow><mrow><mo>∂</mo><msup><mi>z</mi><mn>2</mn></msup></mrow></mfrac><mo>+</mo><mfrac><mrow><mn>1</mn><mo>+</mo><mtext>ν</mtext></mrow><mn>2</mn></mfrac><mfrac><mrow><msup><mo>∂</mo><mn>2</mn></msup><mi>k</mi></mrow><mrow><mo>∂</mo><mi>x</mi><mo>∂</mo><mi>z</mi></mrow></mfrac><mo>−</mo><mrow><mo>(</mo><mrow><mn>1</mn><mo>+</mo><mi>ν</mi></mrow><mo>)</mo></mrow><mi>α</mi><mfrac><mrow><mo>∂</mo><mi>T</mi></mrow><mrow><mo>∂</mo><mi>x</mi></mrow></mfrac><mo>=</mo><mn>0</mn></mrow></mtd></mtr><mtr><mtd><mrow><mfrac><mrow><msup><mo>∂</mo><mn>2</mn></msup><mi>k</mi></mrow><mrow><mo>∂</mo><msup><mi>z</mi><mn>2</mn></msup></mrow></mfrac><mo>+</mo><mfrac><mrow><mn>1</mn><mo>−</mo><mtext>ν</mtext></mrow><mn>2</mn></mfrac><mfrac><mrow><msup><mo>∂</mo><mn>2</mn></msup><mi>k</mi></mrow><mrow><mo>∂</mo><msup><mi>x</mi><mn>2</mn></msup></mrow></mfrac><mo>+</mo><mfrac><mrow><mn>1</mn><mo>+</mo><mtext>ν</mtext></mrow><mn>2</mn></mfrac><mfrac><mrow><msup><mo>∂</mo><mn>2</mn></msup><mi>i</mi></mrow><mrow><mo>∂</mo><mi>x</mi><mo>∂</mo><mi>z</mi></mrow></mfrac><mo>−</mo><mrow><mo>(</mo><mrow><mn>1</mn><mo>+</mo><mi>ν</mi></mrow><mo>)</mo></mrow><mi>α</mi><mfrac><mrow><mo>∂</mo><mi>T</mi></mrow><mrow><mo>∂</mo><mi>x</mi></mrow></mfrac><mo>=</mo><mn>0</mn></mrow></mtd></mtr></mtable></mrow></mrow></mrow></math> The geometric, physical, and equilibrium equations refer to the theory of elasticity [ 27 27 ] . The corresponding stress component can be calculated by substituting the temperature boundary conditions into Eqs. (17) and (18). In this study, SS304 steel was used as the wall surface of the metal cask in the simulation calculation, and zirconium was used as the lattice material. Table 4 Table 4 presents the specific parameters. 1 1 Name 1 1 Density <math id="M62"><mrow><mrow><mo>(</mo><mrow><mtext>kg</mtext><mo>/</mo><msup><mtext>m</mtext><mn>3</mn></msup></mrow><mo>)</mo></mrow></mrow></math> 1 1 Poisson's ratio 1 1 Young's modulus <math id="M63"><mrow><mo stretchy="false">(</mo><mo>×</mo><msup><mrow><mn>10</mn></mrow><mrow><mn>10</mn></mrow></msup><mtext>Pa</mtext></mrow></math> ) 1 1 Thermal conductivity( <math id="M64"><mrow><mtext>W</mtext><mo>/</mo><mrow><mo>(</mo><mrow><mtext>m</mtext><mo>·</mo><mtext>K</mtext></mrow><mo>)</mo></mrow><mo stretchy="false">)</mo></mrow></math> 1 1 SS304 1 1 7900 1 1 0.29 1 1 19.3 1 1 16.2 1 1 Zr 1 1 6505 1 1 0.34 1 1 9.33 1 1 22 Results and discussion 1 Influence of internal and external power 2 The transient calculation results are as follows. Assuming that the initial temperature of the metal container is 550 K after the fuel is stored for some time, we used the results of the approximate convergence (temperature fluctuation ± 2 K) for comparison, and the transient temperature field was used as the subsequent stress calculation load. Fig.7 Fig.7 illustrates the temperature distribution when the power ratio in the inner and outer zones was 0.5 and 2. At a power ratio of 0.5, the maximum temperature was 592.59 K, located on assembly 14. Nevertheless, the power of assembly 14 in the inner region was lower than that in the outer region. Moreover, the outward radiation power of component 3 was greater than that of the adjacent components (approximately two times), owing to the strong convective heat transfer on the front wall of component 3, located on the same side of the inlet and outlet. Fig.7(b) Fig.7(b) shows that when the power ratio is 2, the maximum temperature is 600.54 K, located in assembly 14. In contrast with the literature, the internal temperature of the component is higher than that of the external component before reaching the steady state, owing to the external component being more susceptible to the natural convection of helium gas and air in the external annular channel than the internal component. After changing the power, the temperature became lower than the power ratio of 1 (609 K). The temperature distribution in Fig.7(a) Fig.7(a) is closer to the container wall than that in Fig.7(b) Fig.7(b) , and the temperature in Fig.7(a) Fig.7(a) is significantly higher than Fig.7(b) Fig.7(b) near the top of the container in the inner region, which also reflects that when the internal power is smaller than the external power, the outer component not only dissipates heat but also conducts heat inward. Figure 8 Figure 8 illustrates the temperature distribution of the frame and the component with the highest temperature in the model under the analyzed power. Figure 9 Figure 9 shows the axial temperature distribution of the internal and external assemblies (14, 3) under different powers. Fig.8(a) Fig.8(a) transverse temperature change at the lower side near the bottom of the container is due to the natural convection influence of helium in the container. The exterior temperature is lower than the interior temperature, and helium flows inward from the outer wall. In Fig.8(b)(c) Fig.8(b)(c) , power distribution is more uneven than that in Fig.8(a) Fig.8(a) , which is affected by helium convection, and the high temperature shifts upward. Figure 9 Figure 9 displays that in the case of even power distribution, component gradient varieties in the inner and outer zones are primarily located at the upper and lower sides of the component. When the power distribution is different, the variation in the temperature gradient in the inner and outer zones of the upper side is significantly smaller than that in the case of uniform distribution. On the lower side, because of the influence of heat transfer in the outer zone under different power conditions, the internal temperature is always higher than that in the outer zone. The change in the temperature gradient in the outer zone is greater than that in the inner zone. The axial temperature gradient is caused by the filling gas in the cask and the cooling air in the concrete cladding. Figure 10 Figure 10 illustrates the partial temperature distribution outside the metal cask and through the shell at the height of approximately 4.8 m. The radial temperature gradient through the exterior concrete shell is minimal, which is different from the steady-state situation, owing to the minimal thermal conductivity of concrete and affected by air convection; thus, the temperature change is not obvious. The radial temperature changes along the cask wall and the annular gap are evident. The gradient between the basket and the cask wall is caused by the inactivity of helium gas and its heat conduction to the cask wall. The difference between the steady and the transient states is attributed to a distinct temperature increase at the annular gap due to the steady-state mainstream temperature being higher than the temperature of the outer boundary flow. In this study, the stress distribution inside the model was estimated by the equivalent force, the physical meaning of which is that when the equivalent force at a point of the stressed object exceeds the yield strength of the material, the material begins to undergo plastic deformation, which is defined by the following equation: <math display="block" id="M65"><mrow><msub><mi>σ</mi><mtext>e</mtext></msub><mo>=</mo><msqrt><mrow><mfrac><mn>1</mn><mn>2</mn></mfrac><mrow><mo>[</mo><mrow><msup><mrow><mrow><mo>(</mo><mrow><msub><mi>σ</mi><mn>1</mn></msub><mo>−</mo><msub><mi>σ</mi><mn>2</mn></msub></mrow><mo>)</mo></mrow></mrow><mn>2</mn></msup><mo>+</mo><msup><mrow><mrow><mo>(</mo><mrow><msub><mi>σ</mi><mn>2</mn></msub><mo>−</mo><msub><mi>σ</mi><mn>3</mn></msub></mrow><mo>)</mo></mrow></mrow><mn>2</mn></msup><mo>+</mo><msup><mrow><mrow><mo>(</mo><mrow><msub><mi>σ</mi><mn>3</mn></msub><mo>−</mo><msub><mi>σ</mi><mn>1</mn></msub></mrow><mo>)</mo></mrow></mrow><mn>2</mn></msup></mrow><mo>]</mo></mrow></mrow></msqrt></mrow></math> where <math id="M66"><mrow><msub><mi>σ</mi><mtext>e</mtext></msub><mo> </mo></mrow></math> is the equivalent force and <math id="M67"><mrow><msub><mi>σ</mi><mn>1</mn></msub><mo> </mo></mrow></math> , <math id="M68"><mrow><msub><mi>σ</mi><mn>2</mn></msub></mrow></math> , and <math id="M69"><mrow><mo> </mo><msub><mi>σ</mi><mn>3</mn></msub></mrow></math> are the first, second, and third principal stresses, respectively. For this model, the residual stress is disregarded, and only the thermal stress is considered. After calculations, we found that the thermal stresses on the shell surface and the basket changed little (varying by approximately 0.01 Pa) when the power ratios in the inner and outer zones differed ( Fig.11 Fig.11 and Table 5 Table 5 ). 1 1 Power ratio 1 1 Name 1 1 Total deformation (m) 1 1 Max equivalent stress (Pa) 1 1 Min equivalent stress (Pa) 1 1 1 1 1 Basket 1 1 2.59×10 -8 1 1 31.62 1 1 0.0073 1 1   1 1 Shell 1 1   1 1 5269.53 1 1 714.09 1 1 0.5 1 1 Basket 1 1 2.59×10 -8 1 1 31.63 1 1 0.0073 1 1   1 1 Shell 1 1   1 1 5269.53 1 1 714.09 1 1 2 1 1 Basket 1 1 2.59×10 -8 1 1 31.61 1 1 0.0073 1 1   1 1 Shell 1 1   1 1 5269.53 1 1 714.09 As shown in Fig.12 Fig.12 , where the maximum equivalent stress on the shell was 5269.5 Pa located at the intersection of the shell and concrete, the maximum circumferential stress was -6325.9 Pa (counterclockwise is positive) located on the lower side of the container wall with a total deformation of 2.59 × 10 -8 m. The maximum equivalent stress (31.62 Pa) of the basket was located at the bottom. Fig.12 Fig.12 shows that the stress was concentrated on the lower side of the shell. The temperature stratification on the lower side was primarily affected by the convection of the inlet air. Under the influence of the annular gap, the interior flow velocity of the air below is less than the inlet flow velocity, which cannot drive the inner air, and forming a vortex on the inner side is easy. However, a fluid with high speed at the inlet moves upward when heated. Regarding the thermal stress and physical parameters of the material (e.g., modulus of elasticity, coefficient of thermal expansion) and temperature change function, the choice of material physical parameters was relatively stable; thus, the total power remained unchanged, and the temperature gradient in the container directly affected thermal stress. The basket was affected by the temperature of the assembly, and the high temperature was concentrated above the basket, but the temperature difference between the upper and lower was minor (maximum of approximately 10 K), and the equivalent stress was comparatively small (approximately 6 Pa). As shown in Fig.12 Fig.12 , the stress change was apparent at the temperature gradient that altered the most. Therefore, we selected the specified path in Figs. 13 Figs. 13 , 14 14 , 15 15 to calculate the equivalent stress, normal stress (positive for tensile stress and negative for compressive stress), and circumferential stress. Figure 13(b) Figure 13(b) implies that when the power ratio of the inner and outer zones differs, the temperature variation trend is identical, and the temperature mutation is caused by the temperature accumulation on the upper side of the cask, which corresponds to Fig.13(a) Fig.13(a) Stress mutation is observed, but the trend is slight, and then all stresses increase evenly. Figure 13(a) Figure 13(a) the end of the path is also affected by the bottom stress, which makes the stress near the edge abrupt. As shown in Fig.14 Fig.14 , the path was at the bottom edge of the shell, and the stress was influenced by temperature fluctuations. The stress trend on the path first grew and then declined along the path direction, corresponding to its temperature distribution, and the temperature difference between the two sides was lower than the middle difference. Fig.15 Fig.15 illustrates that the temperature above the assembly was high, and the variety was minor. As the temperature decreased, stress increased. The temperature gradient changes in the inner assembly were small; thus, stress change was relatively stable. Therefore, we concluded that the temperature gradient is the primary factor of thermal stress on the shell, and the variation in the power ratio in the inner and outer zones with a constant total power only affects the temperature distribution of the assemblies, which has little effect on stress. The maximum equivalent stress estimated in the shell is less than the yield strength of the material (205 MPa); therefore, the container fulfills the safety requirements under this hypothetical condition. Effect of inclination angle 2 In the process of spent fuel shipping and storage, the inclination is bound to occur. Taking the internal and external power ratio of 1.0 as the basis case (θ=0°calculated above), we calculated the temperature field and the corresponding stress at inclination θ =15°, 45°, 75°, and 90° as follows: Figure 16 Figure 16 shows that shell temperature presents a wavy distribution affected by the tilt angle, with the peak and valley shifting from right to left and the peak appearing at the left edge. The bottom temperature is lower than the top temperature, and the top temperature slowly decreases when the tilt angle becomes larger. The leftmost velocity vector diagram displays that the inlet produces obvious vortices on the left side below θ=90°, near the left side of the velocity is higher than the right side. When this is close to the top of the container, the velocity of the air on the left side decreases under the influence of the annular gap and gradually increases toward the right outlet. Combined with Fig.12 Fig.12 , the container tilt causes variations in the internal air distribution, and air in the annular gap due to the thermal pressure influences the upward flow. Part of the hot air impinged on the wall, decreasing speed, which could not drive the inlet’s cold air to form a vortex below. With a gradual increase in the tilt angle, the flow velocity on the left side (the top in the vertical direction) gradually increased, and velocity stratification appeared, which enhanced the heat transfer on the shell, exhibiting the aforementioned varieties. Figure 17 Figure 17 demonstrates that in the tilting process, the stress was predominantly concentrated on the top and bottom edges of the shell. When the container was placed horizontally ( θ =90°), the deformation was uniformly distributed from the middle to the left, and the basket deformation with the increase in tilt angle gradually shifted from the bottom to upward. The basket maximum deformation emerged in the outermost surface. The outermost side of the basket was affected by the natural convection of the inner and outer sides, and the temperature was distributed on the left side (top of the vertical direction). Tables 6 Tables 6 and 7 7 show the maximum deformation, equivalent stresses, normal stresses (positive along the shell upward and negative downward), circumferential stresses (positive counterclockwise), and tensile stresses (negative compressive stresses) for the shell and basket. 1 1 θ 1 1 Total deformation (m) 1 1 Equivalent stress (Pa) 1 1 Normal stress(Pa) 1 1 Circumferential stress(Pa) 1 1 15° 1 1 6.68×10 -5 1 1 1.91×10 8 1 1 -1.66×10 8 1 1 -8.89×10 7 1 1 45° 1 1 6.68×10 -5 1 1 2.02×10 8 1 1 -1.70×10 8 1 1 -9.26×10 7 1 1 75° 1 1 6.64×10 -5 1 1 1.87×10 8 1 1 -1.86×10 8 1 1 -1.86×10 8 1 1 90° 1 1 3.12×10 -9 1 1 741.00 1 1 -458.84 1 1 948.30 1 1 θ 1 1 Total deformation (m) 1 1 Equivalent stress (Pa) 1 1 Normal stress (Pa) 1 1 Tensile stress (Pa) 1 1 15° 1 1 2.12×10 -11 1 1 30.40 1 1 -20.06 1 1 15.23 1 1 45° 1 1 2.40×10 -11 1 1 22.96 1 1 -19.65 1 1 11.83 1 1 75° 1 1 3.50×10 -11 1 1 7.77 1 1 5.43 1 1 8.09 1 1 90° 1 1 3.80×10 -11 1 1 5.21 1 1 6.60 1 1 6.65 We determined that the temperature gradient varieties directly involved stress. Figure 18 Figure 18 illustrates the circumferential and equivalent stresses along the path on the outer surface of the shell. Figure 18(a) Figure 18(a) shows the temperature variation along the path. When θ =45° and 75°, the temperature change on the path was distinct, and the temperature gradient on the lower side grew abruptly near the container. The equivalent stress initially decreased, then tended to be a steady value of approximately 30 kPa ( θ =90°, approximately 60 Pa), and quickly increased at the end. The trends of the different inclination angles were identical. The circumferential stress shifted slowly from top to bottom, and the circumferential stress increased with an increase in the temperature gradient at the end of the path. The change in wall temperature was mainly caused by the distinction of the air distribution, which leads to a distinction of the gradient below, increasing toroidal stress. The equivalent stress at the end of the path differs from the circumferential stress because the former is affected by the stress components in other directions. The end of the path is also influenced by the stress on other surfaces when it is located at the intersection of two surfaces. Figure 19(b) Figure 19(b) displays the temperature variations on the path of the assembly, and Fig.19(a) Fig.19(a) illustrates the path stress change on the corresponding basket. During the tilt process (except θ =0°, 90°), the temperature variation trend was identical, and the upper and lower temperature distinction was tiny (approximately 1-2 ℃). Due to the influence of the inclination angle, the external air distribution changes, and heat transfer is enhanced; thus, the basis case ( θ =0°) temperature was higher than in other conditions. When θ =90° (the container was placed horizontally), the fluid in the component could only flow to both sides, and the convective heat transfer ability was limited. Therefore, the overall temperature of the component higher than θ =15°, 45°, 75°, and stress was the least. With increasing tilt angle, the normal stress of inner assembly 14 (positive along the basket) decreased progressively and increased gradually along the radial direction of the path. The stress in the basket increased approximately with an increase in the temperature gradient. Figure 17 Figure 17 shows that the stress was concentrated at the top and bottom of the shell. Combined with the bottom deformation of θ =75° in Fig.20 Fig.20 , we found that the stress is concentrated near the edge. Therefore, the path in Fig.20 Fig.20 was selected as the top and bottom edges. In Fig.20(a) Fig.20(a) , the equivalent stress of path1 (bottom edge) exhibited a slight fluctuation and was smaller than that of path2. The equivalent stresses of Path1 and Path2 first declined and then grew with fluctuation variety. Fig.20(b) Fig.20(b) shows that the circumferential stress changed in a trough shape, and the circumferential stress of path2 was less than that of path1 at 0.3-0.4 m of the path. Path1 and path2 were at the intersection of surfaces, and the stress was not only affected by temperature variation but also by the stress of adjacent surfaces, leading to stress fluctuations. The aforementioned analysis of the tilt condition demonstrates that the horizontal state ( θ =90°) differed from the other states, and the stress magnitude was near to the basic case ( θ =0°). During the tilt process, the deformation steadily increased from below and was eventually uniformly distributed on the shell. When θ =45°, the maximum equivalent stress was 202 MPa (less than the yield stress of 205 MPa), and the container fulfills the safety requirements in the tilt process. The tilt predominantly alters the air distribution in the annular gap, which results in variations in the shell temperature distribution and various stresses. Therefore, it is concluded that altering the air distribution is a major cause of temperature distribution and stress variation. Notably, air distribution is also affected by the location, size of the inlet and outlet, fluid velocity, and shape of the air passage, which requires follow-up research. Conclusion 1 In this study, a 3D transient thermal coupling model was used to calculate the temperature and stress fields, and the following conclusions were drawn: 1. The temperature gradient is the primary cause of thermal stress on the wall, and the variation in the power ratio with constant total power only affects the temperature distribution of the component and has little effect on the stress (change is approximately 0.01 Pa). 2. Tilt alters the airflow distribution in the annular gap and enhances heat transfer, which is the principal factor of temperature distribution and stress change. 3. Maximum equivalent stress (202 MPa) is located at the shell and is less than the yield strength of the corresponding material (205 MPa), satisfying the safety requirements. In summary, we focused on the unidirectional coupling of heat transfer and thermal stress. We plan to extend the interaction of heat transfer and thermal stress to the transient stress calculation of the complete process, and the accident and anomaly effects will be considered. In addition, the similarity theory can be used to scale the model and reduce the calculations to realize the whole process analysis, so as to explore the stress change in a single spent fuel rod during the heating process. However, these problems require further study. References IAEA Annual Report for 2021 . ( International Atomic Energy Agency , 2021 ). https://www.iaea.org/publications/reports https://www.iaea.org/publications/reports Annual Report for 2021 M. Z. Bell , A. Macfarlane , Fixing the nuclear waste problem? The new political economy of spent fuel management in the United States . Energy Res. Soc. Sci. 91 , 102728 ( 2022 ). doi: 10.1016/j.erss.2022.102728 http://doi.org/10.1016/j.erss.2022.102728 Fixing the nuclear waste problem? The new political economy of spent fuel management in the United States Y. Can , B. B. Acara , Nuclear proliferation resistance assessment of fuel cycles closed with complete co-processing of spent fuel . Prog. Nucl. Energ. 150 , 104297 ( 2022 ). doi: 10.1016/j.pnucene.2022.104297 http://doi.org/10.1016/j.pnucene.2022.104297 Nuclear proliferation resistance assessment of fuel cycles closed with complete co-processing of spent fuel A. Wahid , T. Sungari , R. Ratio , Dynamic modeling and controlling of a spent nuclear fuel storage pool under periodic operation and station blackout conditions . Ann. Nucl. Energy. 166 , 108751 ( 2022 ). doi: 10.1016/j.anucene.2021.108751 http://doi.org/10.1016/j.anucene.2021.108751 Dynamic modeling and controlling of a spent nuclear fuel storage pool under periodic operation and station blackout conditions G. Rothwell , Spent nuclear fuel storage: What are the relationships between size and cost of the alternatives? Energy Policy. 150 , 112126 ( 2021 ). doi: 10.1016/j.enpol.2020.112126 http://doi.org/10.1016/j.enpol.2020.112126 Spent nuclear fuel storage: What are the relationships between size and cost of the alternatives? W. Wang , B. Wang , Y. Zhou , et al. , Study on the mass and heat transfer characteristics of the drying process of spent nuclear fuel cladding with nitrogen assisted . Nucl. Eng. Des. 391 , 111672 ( 2022 ). doi: 10.1016/j.nucengdes.2022.111672 http://doi.org/10.1016/j.nucengdes.2022.111672 Study on the mass and heat transfer characteristics of the drying process of spent nuclear fuel cladding with nitrogen assisted Y. S. Chen , Thermal analysis for the integrated spent fuel pool of the Chinshan plant in the decommissioning process . Ann. Nucl. Energy 119 , 163 - 174 ( 2018 ). doi: 10.1016/j.anucene.2018.05.005 http://doi.org/10.1016/j.anucene.2018.05.005 Thermal analysis for the integrated spent fuel pool of the Chinshan plant in the decommissioning process R. Ratiko , D.S. Wisnubroto , N. Nasruddin , et al. , Current and future strategies for spent nuclear fuel management in Indonesia . Energy Strateg. Rev. 32 , 100575 ( 2020 ). doi: 10.1016/j.esr.2020.100575 http://doi.org/10.1016/j.esr.2020.100575 Current and future strategies for spent nuclear fuel management in Indonesia D.-W. Lim , C.-W. Lee , J.-Y. Lim et al. , On the Particle Swarm Optimization of cask shielding design for a prototype Sodium-cooled Fast Reactor . Nucl. Eng. Technol. 51 ( 1 ), 284 - 292 ( 2019 ). doi: 10.1016/j.net.2018.09.007 http://doi.org/10.1016/j.net.2018.09.007 On the Particle Swarm Optimization of cask shielding design for a prototype Sodium-cooled Fast Reactor M. G. El-Samrah , A.F. Tawfic , Samir E. Chidiac , Spent nuclear fuel interim dry storage; Design requirements, most common methods, and evolution: A review . Ann. Nucl. Energy 160 , 108408 ( 2021 ). doi: 10.1016/j.anucene.2021.108408 http://doi.org/10.1016/j.anucene.2021.108408 Spent nuclear fuel interim dry storage; Design requirements, most common methods, and evolution: A review Y. Saito , J. Kishimoto , T. Matsuoka , et al. , Containment integrity evaluation of MSF-type cask for interim storage and transport of PWR spent fuel . Int. J. Pres. Ves. Pip. 117 , 33 - 41 ( 2014 ). doi: 10.1016/j.ijpvp.2013.10.007 http://doi.org/10.1016/j.ijpvp.2013.10.007 Containment integrity evaluation of MSF-type cask for interim storage and transport of PWR spent fuel C.-C. Hung , F. King , Y.-C. Yu , et al. , Corrosion assessment for spent nuclear fuel disposal in crystalline rock, using variant cases of hydrogeological modeling . Nucl. Sci. Tech. 31 , 116 ( 2020 ). doi: 10.1007/s41365-020-00822-6 http://doi.org/10.1007/s41365-020-00822-6 Corrosion assessment for spent nuclear fuel disposal in crystalline rock, using variant cases of hydrogeological modeling J. Oh , S. Kwag , J. Lee , A new design concept and seismic margin assessment for a spent fuel storage system . Nucl. Eng. Des. 326 , 150 - 161 ( 2018 ). doi: 10.1016/j.nucengdes.2017.11.011 http://doi.org/10.1016/j.nucengdes.2017.11.011 A new design concept and seismic margin assessment for a spent fuel storage system R. Poškas , P. Poškas , K. Račkaitis , et al. , A numerical study of thermal behavior of CASTOR RBMK-1500 cask under fire conditions . Nucl. Eng. Des. 376 , 111131 ( 2021 ). doi: 10.1016/j.nucengdes.2021.111131 http://doi.org/10.1016/j.nucengdes.2021.111131 A numerical study of thermal behavior of CASTOR RBMK-1500 cask under fire conditions H. Boucherit , A. Kaliatka , A. Lounis , Integrity control of an RBMK-1500 fuel rod locally oxidized under a bounding reactivity-initiated accident . Nucl. Sci. Tech. 30 ( 1 ), 1 - 11 ( 2019 ). doi: 10.1007/s41365-018-0538-2 http://doi.org/10.1007/s41365-018-0538-2 Integrity control of an RBMK-1500 fuel rod locally oxidized under a bounding reactivity-initiated accident M. Wang , Y. Wang , W. Tian , et al. , Recent progress of CFD applications in PWR thermal hydraulics study and future directions . Ann. Nucl. Energy 150 , 107836 ( 2021 ). doi: 10.1016/j.anucene.2020.107836 http://doi.org/10.1016/j.anucene.2020.107836 Recent progress of CFD applications in PWR thermal hydraulics study and future directions J. Li , M. Wang , D. Fang , et al. , CFD simulation on the transient process of coolant mixing phenomenon in reactor pressure vessel . Ann. Nucl. Energy 153 , 108045 ( 2021 ). doi: 10.1016/j.anucene.2020.108045 http://doi.org/10.1016/j.anucene.2020.108045 CFD simulation on the transient process of coolant mixing phenomenon in reactor pressure vessel M. Wang , L. Wang , X. Wang , et al. , CFD simulation on the flow characteristics in the PWR lower plenum with different internal structures . Nucl. Eng. Des. 364 , 110705 ( 2020 ). doi: 10.1016/j.nucengdes.2020.110705 http://doi.org/10.1016/j.nucengdes.2020.110705 CFD simulation on the flow characteristics in the PWR lower plenum with different internal structures D.-G. Lee , N.-H. Sung , J.-H. Park et al. , An assessment of temperature history on concrete silo dry storage system for CANDU spent fuel . Ann. Nucl. Energy 94 , 263 - 271 ( 2016 ). doi: 10.1016/j.anucene.2016.03.007 http://doi.org/10.1016/j.anucene.2016.03.007 An assessment of temperature history on concrete silo dry storage system for CANDU spent fuel H.S. Yoo , S.H. Yoo , E.S. Kim , Heat transfer enhancement in dry cask storage for nuclear spent fuel using additive high density inert gas . Ann. Nucl. Energy 132 , 108 - 118 ( 2019 ). doi: 10.1016/j.anucene.2019.04.018 http://doi.org/10.1016/j.anucene.2019.04.018 Heat transfer enhancement in dry cask storage for nuclear spent fuel using additive high density inert gas J. Benavides , G. Jimenez , M. Galbán , et al. , Methodology for thermal analysis of spent nuclear fuel dry cask using CFD codes . Ann. Nucl. Energy 133 , 257 - 274 ( 2019 ). doi: 10.1016/j.anucene.2019.05.026 http://doi.org/10.1016/j.anucene.2019.05.026 Methodology for thermal analysis of spent nuclear fuel dry cask using CFD codes R. Poškas , V. Šimonis , P. Poškas , et al. , Thermal analysis of CASTOR RBMK-1500 casks during long-term storage of spent nuclear fuel . Ann. Nucl. Energy 99 , 40 - 46 ( 2017 ). doi: 10.1016/j.anucene.2016.09.031 http://doi.org/10.1016/j.anucene.2016.09.031 Thermal analysis of CASTOR RBMK-1500 casks during long-term storage of spent nuclear fuel J.-C. Lee , K.-S. Bang , S.-H. Yu et al. , Similarity analysis of thermal-fluid flow for thermal testing using scaled-down model of spent fuel storage cask . Ann. Nucl. Energy 149 , 107791 ( 2020 ). doi: 10.1016/j.anucene.2020.107791 http://doi.org/10.1016/j.anucene.2020.107791 Similarity analysis of thermal-fluid flow for thermal testing using scaled-down model of spent fuel storage cask S. Alyokhina , A. Kostikov , Unsteady heat exchange at the dry spent nuclear fuel storage . Nucl. Eng. Technol. 49 ( 7 ), 1457 - 1462 ( 2017 ). doi: 10.1016/j.net.2017.07.029 http://doi.org/10.1016/j.net.2017.07.029 Unsteady heat exchange at the dry spent nuclear fuel storage Y. Wu , J. Klein , H. Zhou , et al. , Thermal and fluid analysis of dry cask storage containers over multiple years of service . Ann. Nucl. Energy 112 , 132 - 142 ( 2018 ). doi: 10.1016/j.anucene.2017.10.013 http://doi.org/10.1016/j.anucene.2017.10.013 Thermal and fluid analysis of dry cask storage containers over multiple years of service L. E. Herranz , J. Penalva , F. Feria , CFD analysis of a cask for spent fuel dry storage: Model fundamentals and sensitivity studies . Ann. Nucl. Energy 76 , 54 - 62 ( 2015 ). doi: 10.1016/j.anucene.2014.09.032 http://doi.org/10.1016/j.anucene.2014.09.032 CFD analysis of a cask for spent fuel dry storage: Model fundamentals and sensitivity studies A. I. Lurie , A. Belyaev , Theory of Elasticity . ( Springer Berlin , Heidelberg , 2005 .) Theory of Elasticity All authors contributed to the study conception and design. Material preparation, data collection and analysis were performed by Wei-Hao Ji, Jian-Jie Cheng, Han-Zhong Tao, Wei Li. The first draft of the manuscript was written by Wei-Hao Ji and all authors commented on previous versions of the manuscript. All authors read and approved the final manuscript. We are grateful to the High-Performance Computing Center of Nanjing Tech University for supporting the computational resources. 10.1007/s41365-023-01171-w 2022-09-16 2022-11-28 2022-12-02 2023-02-24 2023-03-20 © China Science Publishing & Media Ltd. (Science Press), Shanghai Institute of Applied Physics, the Chinese Academy of Sciences, ChineseNuclear Society and Springer Nature Singapore Pte Ltd. 2023 This is a post-peer-review, pre-copyedit version of an article published in Nuclear Science and Techniques. The final authenticated version is available online at: http://dx.doi.org/10.1007/s41365-023-01171-w http://dx.doi.org/10.1007/s41365-023-01171-w . Springer Nature or its licensor (e.g. a society or other partner) holds exclusive rights to this article under a publishing agreement with the author(s) or other rightsholder(s); author self-archiving of the accepted manuscript version of this article is solely governed by the terms of such publishing agreement and applicable law. 2023 Nuclear Science and Techniques 2 34