Material derivative in fluid dynamics

The material or substantial derivative links the Eulerian and Lagrangian descriptions by describing the rate of change of a fluid property (e.g., velocity and density) experienced by a moving particle: (7)where is the material derivative of the property ϕ, u is the velocity field, and is the gradient ϕ. This equation captures the temporal change in ϕ at a fixed location and the change due to the motion of the fluid itself, aligning with the Eulerian viewpoint while describing the Lagrangian particle dynamics.

Continuity and momentum equations

Reynolds Transport Theorem (RTT) plays a crucial role in converting between system (Lagrangian) and control volume (Eulerian) approaches. This allows the application of conservation laws (mass, momentum, and energy) over a control volume that encompasses the moving system, effectively bridging the gap between the descriptions: (8)where V(t) and S(t) are the time-dependent volume and surface enclosing the fluid or structure, u is the velocity of the fluid relative to surface S; and n is the outward unit normal to S.

Coupling mechanism via virtual work

The coupling of fluid and structure in FSIs through the virtual work principle involves ensuring that the sum of all virtual work contributions equals zero, which reflects the equilibrium state. This principle can be expressed through the virtual work equation, which, in its generalized form for FSIs, is as follows: (9)where: V_f and V_s denote fluid and solid volumes, respectively. ρ_f and ρ_s are the fluid and solid densities, respectively. a_f and a_s are the acceleration fields of the fluid and solid, respectively. and δu_s represent the virtual displacements of the fluid and solid, respectively. denotes the stress tensor of the solid. denotes the virtual strain tensor of the solid. S_f is the surface on which the external tractions t_f act on the fluid. b_f represents the body force per unit mass (e.g., gravity) acting on a fluid. S_fs is the interface surface between the fluid and solid. f_c is the contact force of the fluid-solid interface per unit area. δu_fs represents the virtual displacement at the fluid-solid interface.

The penalty function method is applied at the fluid–structure interface (S_fs) to enforce the compatibility and no-penetration constraints as (10)where k_c denotes the penalty stiffness coefficient that enforces the constraint, Δu_fs denotes the relative (penetration) displacement between the fluid and solid at the interface, which should be minimized.

This equation ensures that any penetration between the fluid and solid is resisted by a force proportional to the penetration depth, effectively coupling the fluid and solid domains through the virtual work performed by the contact forces. This mechanism allows for consistent and physically accurate interactions between the fluid and solid phases in the simulations, enabling the practical analysis of complex FSI within the HEL method framework.

Configuration of optimization

The response-spectrum method is a foundational analytical tool in the seismic analysis of liquid-containing structures for nuclear reactor systems. In this study, we employ a parametric approach to investigate the sensitivity of the seismic response to key structural and liquid parameters using response-spectrum analysis, a method widely applied in seismic optimization studies [62]. This study is focused on vertical storage tanks, which are critical components of nuclear facilities, and evaluates three primary parameters: liquid level, tank height, and wall thickness. This parametric study aims to quantify the influence of structural parameters on the seismic response, while simultaneously assessing the impact of liquid parameter variations on the seismic behavior of vertical tanks.

Figure 2a illustrates the study’s structural analysis and optimization flowchart. The structural system under investigation comprises a tank body, skirt supports, and upper and lower butterfly heads. The focus of this research is on the wall thickness, which varies from 5 mm to 60 mm, explicitly targeting the sections of the tank in direct contact with the contained liquid. The tank height is adjusted from 7 m to 12 m; however, these modifications do not change the support structures. Specifically, the height of the skirt support structure is 0.6 m, and its thickness is 50 mm throughout the analysis. The design strategically sets the liquid level within the tank to accommodate a range between 0 m and 6.4 m, ensuring that it covers all potential tank volumes. Utilizing the central composite design (CCD), a comprehensive set of parametric design models is generated. These models are evaluated using the kriging response surface methodology, which is a robust technique for modeling and analyzing complex systems [63, 64]. The evaluation criteria encompass a variety of indicators, such as modal shape and frequency, maximum equivalent stress sustained by the structure, and total displacement under seismic conditions.

Fig. 2Full-screen View

(Color online) Response of liquid-containing structures during seismic event using the response-spectrum method. (a) Structural analysis and optimization flow chart; (b) Boundary conditions of the response-spectrum analysis; (c) Horizontal and vertical response spectrum of 5% damping; (d) Dominant structural modes with significant modal mass ratios in principal directions

The loads impacting the structural system include gravity and hydrostatic pressure, with the structural base secured by fixed supports, as shown in Fig. 2b. Seismic response analysis is conducted based on the response spectra of the equipment close to the ground, as shown in Fig. 2c [65-67]. While conventional structural design often utilizes standard code-based spectra, such as those from ASCE to 7-16 [57], in this study, we employ a specific spectrum obtained from ITER documentation. This approach is justified by the unique requirements of nuclear facilities, which necessitate more stringent seismic safety considerations that align with international standards for nuclear safety, such as those published by the International Atomic Energy Agency (IAEA) [68].

Figure 2c illustrates the seismic response spectrum derived from ITER documentation and used in this analysis. The spectrum incorporates two seismic levels: Seismic Level 1 (SL-1) and Seismic Level 2 (SL-2), which are defined by the IAEA for nuclear power plant design. SL-2, equivalent to the Safety Shutdown Earthquake, represents the maximum potential earthquake for which the plant must safely shut down. SL-1, corresponding to the Operating Basis Earthquake, denotes the seismic level at which the plant should remain operational. In the context of ITER, SL-2 is used to design safety-important structures, systems, and components, whereas SL-1 is applied to conventional structures and systems. Due to the rotational axis symmetry inherent in the structural design, the seismic responses in the X and Y directions are assumed to remain consistent under the same seismic intensity. This assumption is crucial as the focus is on the methodologies, not on the impact of varying response-spectrum amplitude or parameters. Selecting a universal response spectrum ensures that the findings are widely applicable and relevant to the method’s generalizability while adhering to the specialized requirements of nuclear facility design.

To further characterize the structural behavior under seismic conditions, a modal analysis was performed. Figure 2d shows the modal analysis results, highlighting the dominant structural modes in the principal directions. The analysis reveals the modes with significant modal mass ratios, providing insight into the structure’s dynamic behavior under prestressed conditions.

The material properties of both the solid and fluid components, which are critical for the seismic analysis, are presented in Table 3. Figure 3 shows the minimum and average true stress–strain curves for the 304 stainless steel, providing critical insights into the material’s behavior under dynamic conditions.

Fig. 3Full-screen View

(Color online) Minimum and average true stress–strain of 316L

The response-spectrum analysis method for assessing the seismic response of storage tanks in nuclear reactors has provided insights into their dynamic behavior. We used complete quadratic combination and square root of sum of squares to aggregate periodic vibration modes, while the missing mass method summed rigid body modes. These methods enable analysis by accounting for multiple vibrational modes. The analysis, considering frequencies from 0 Hz to 240 Hz, includes periodic and rigid modes. The selected frequency range ensures that the structural natural frequency remains distant from the seismic frequency of 33 Hz, effectively reducing the impact of seismic events.

Results and discussion

The sensitivity of the structural response to various parameters under different design conditions is depicted in Fig. 4. This sensitivity distribution diagram illustrates the positive and negative correlations of parameters such as tank height and liquid level with the seismic response indicators. The tank height significantly influences the evaluation and optimization parameters, whereas the liquid level has minimal impact on the structural modes. This observation aligns with the results obtained from the prestressed modal analysis, indicating that the static load response predominantly determines the structural response under seismic conditions.

Fig. 4Full-screen View

(Color online) Sensitivity distribution diagram

The analysis results are clearly depicted on the response surface, as illustrated in Fig. 5. This visualization delineates the influence of each study parameter on the seismic performance targets. These insights provide a quantitative basis for the ongoing optimization of tank structural designs, with the aim of enhancing their safety and efficacy under seismic conditions.

Fig. 5Full-screen View

(Color online) Structural response from response-spectrum analysis with error bars representing variability due to liquid-level height. (a) First mode frequency; (b) Total displacement; (c) Equivalent stress

The tank wall thickness and height response surfaces against seismic performance indicators, such as the first-order mode, structural displacement, and stress, are depicted in Fig. 5a-Fig. 5c. These visualizations show the influence of parameters on seismic performance targets. Figure 5a shows that tank body height dominates the first-order natural frequency under prestress. When the height is high enough, the first-order natural frequency will be close to or lower than the rigid body frequency of 33 Hz, implying that earthquakes have negligible effect on low-height tanks. Combining Fig. 5b and Fig. 5c, we observe that the normalized seismic response under stress prestressing mode can be represented by maximum relative displacement and maximum equivalent stress. This result does not combine structural self-weight condition responses, only seismic ones. These insights provide a basis for optimizing tank structural designs to enhance safety under seismic conditions. By analyzing these parameters, this investigation characterizes the tank’s seismic response, contributing to understanding dynamic behavior of liquid-containing structures in nuclear facilities.

We identified the optimal structural parameters for tank design: a height of 7 m corresponding to highest natural frequency, wall thickness of 30 mm for structural integrity, and liquid level of 4.8 m (75% of internal height). These parameters, selected based on preliminary tests, provide a foundation for future tank design and safety assessments. The modal mass analysis shown in Fig. 6 reveals the concentration of the modal mass within the first two modes on a logarithmic scale. This concentration indicates the initial modes’ critical role in the seismic response, with the dominant mode at a characteristic frequency of 145 Hz classified as rigid owing to its high frequency relative to the periodic modes.

Fig. 6Full-screen View

(Color online) The relative modal mass as a function of eigenfrequency

The tank’s maximum displacement and equivalent stress levels are analyzed, and the specific locations of these maxima are identified and documented in Table 4. These findings provide a basis for evaluating structural safety and reliability according to the criteria of the RCC-MR [69] code or ASME VIII Div. 1 [70] and Div. 2 [71]. This evaluation is crucial to ensure the integrity and safety of nuclear reactor components under seismic conditions.

Implementation of the HEL method model

This study utilizes the HEL method to evaluate the dynamic responses of a liquid-containing structure under seismic excitation. This modeling approach divides the system into two computational domains: the Euler domain, which represents the fluid dynamics within the tank, and the Lagrangian domain, which simulates the structural dynamics of the tank, as shown in Fig. 7. The integration of these domains enables a detailed representation of FSI under seismic loading, capturing both the complex motion of the fluid and structural deformation.

Fig. 7Full-screen View

(Color online) Seismic time-history analysis model for structural dynamics. (a) Boundary conditions and markers for the three observation points; (b) Lagrangian domain; (c) Eulerian domain

The Euler domain simulates fluid behavior under seismic excitation through a finely meshed Eulerian grid, capturing wave propagation, sloshing, and pressure distribution. The Lagrangian domain models the tank’s structural response using a refined mesh for precise simulation of deformation, stress distribution, and stability during seismic events.

Three observation points monitor the seismic response: the top observation point, liquid-level side point, and dynamic peak point. The top observation point, located at the tank’s highest edge, tracks global displacement and deformation under seismic excitation, providing insights into structural performance and stability.

The liquid-level side point corresponds to the maximum liquid level under dynamic conditions, moving with changes in liquid height caused by seismic forces. The point position tracks within ±25% of the static initial liquid height to represent the liquid’s fluctuating behavior throughout the seismic response.

The dynamic peak point represents the location within the Lagrangian domain where the maximum structural response occurs under seismic excitation. This point varies with time depending on the dynamic interaction between the fluid and structure, capturing the peak displacement or stress experienced by the tank walls during the seismic event. Identifying this point is crucial for assessing the structural safety and integrity of the tank as it provides insights into the areas experiencing the highest strain and potential failure risks.

By carefully defining these three observation points, the HEL method can effectively capture the global and localized seismic responses of the system, thereby enabling a comprehensive evaluation of the FSI. This approach provides a robust and secure framework to accurately simulate the complex dynamics of liquid-containing structures subjected to seismic loading.

The primary function of the Euler domain is to simulate fluid behavior under seismic excitation, emphasizing the need for a finely meshed Eulerian grid to accurately capture intricate fluid phenomena. Conversely, the structural response of the tank encapsulated within the Lagrangian domain necessitates a similarly detailed mesh to ensure the precision of the structural dynamics simulations.

The synergy between these domains through Euler-Lagrange coupling offers a nuanced understanding of the interplay between fluid dynamics and structural integrity under seismic loading. This methodology enhances the accuracy of seismic response predictions and provides a comprehensive framework for analyzing the complex interactions inherent in fluid–structure systems.

Grid independence analysis

Grid independence analysis is critical in numerical simulations to ensure mesh size does not unduly influence results. The objective is to identify an optimal mesh resolution balancing computational efficiency and solution accuracy. We used the response surface method to assess mesh size impacts in Eulerian and Lagrangian domains on FSI dynamics within the HEL method. The simulation involves a tank filled to 75% with liquid under gravitational loading to analyze fluid–structure interaction. The HEL method captures this dynamic interaction effectively. Figure 8 shows the grid independence analysis results, illustrating mesh size effects on calculation cost and structural response, focusing on maximum displacement and stress.

Fig. 8Full-screen View

(Color online) Response of liquid-containing structures during the seismic event using the HEL method. (a) Calculation cost; (b) Total displacement; (c) Equivalent stress

Three key metrics evaluated mesh adequacy: calculation cost, maximum displacement, and maximum stress. A static analysis under hydrostatic pressure and gravity was performed alongside the dynamic HEL simulation, serving as a baseline for comparison, represented by the blue curve in the X-Z plane in Fig. 8. The results show mesh resolution’s importance in the Lagrangian domain for accurate structural response. In both methods, maximum displacement and stress occurred at the tank bottom’s geometric center, though magnitudes varied between analyses due to different loading conditions. The calculation cost of the HEL method FSI simulation, shown in Fig. 8a, depends heavily on mesh sizes in both domains. When the Eulerian domain has small mesh while the Lagrangian domain remains coarse, computational time increases significantly with potential convergence issues. This discrepancy is observed because the Eulerian domain needs finer resolution for accurate fluid dynamics modeling. As shown in Fig. 8b and Fig. 8c, while simulation shows approximate convergence at larger grid sizes around 800 mm, the convergence is weak and rapidly diverges with further refinement. This phenomenon suggests that coarse meshes may not be sufficient to resolve the intricate FSI accurately, particularly in transient analyses, as larger mesh sizes smooth out critical interaction details (leading to inaccurate predictions of structural responses).

As the mesh is refined, the grid sizes of Eulerian and Lagrangian domains must be matched and sufficiently refined for convergence and accuracy in HEL simulations. A mesh size of approximately 200 mm in both the domains balances the calculation cost and accuracy. Optimal results are achieved when the Eulerian domain uses a mesh size of 100 mm × 100 mm × 100 mm and the Lagrangian domain uses 80 mm × 80 mm × 80 mm. This correspondence between grid sizes is crucial for accurate fluid–structure coupling and reliable outcomes.

Figure 9 presents tank liquid levels from empty to high, designated as L1 through L4. This classification enables assessment of how liquid levels affect dynamic response and structural stresses under seismic conditions. The analysis is crucial, as liquid presence can alter stress distribution and deformation patterns. This classification helps formulate safety measures and enhance liquid storage facility design in nuclear reactors.

Fig. 9Full-screen View

(Color online) Classification of the structure at various liquid levels. (a) Level L1 (empty), 0.0 m; (b) Level L2 (low), 1.6 m; (c) Level L3 (half), 3.2 m; (d) Level L4 (high), 4.8 m

The modeling area consists of two domains, as shown in Fig. 7: the “Tank” represents the Lagrangian domain, while the “Whole fluid domain” represents the Euler domain encapsulating it. As detailed in Table 5, the grid sizes in both domains are matched to maintain solution accuracy. The Lagrangian domain contains two layers of grids in the thickness direction to capture the tank’s structural dynamics during seismic response. Liquid phases L1 through L4 in Table 5 serve only as initialization parameters and reflect material volume fractions within the mesh. The Euler domain’s grid size remains constant while fluids flow between grids, interacting with structure and fluid, whereas the Lagrangian domain follows material motion. The Courant—Friedrichs–Lewy (CFL) condition stays below one to ensure numerical stability.

Ground motion acceleration synthesis

Ground motion selection is crucial in time-history dynamic analysis, with multiple recorded motions used for average structural responses. We used synthetic records from the ITER-specific response spectrum, ensuring consistency with the response-spectrum analysis and confirming direct comparability between the HEL and response-spectrum methods. The ITER spectrum reflects seismic hazard considerations for nuclear facilities that recorded motions may not represent. Synthetic records enable control over frequency content and duration, ensuring alignment with the design spectrum. For this study’s specific site and hazard level, recorded ground motions may inadequately match the target spectrum without excessive scaling. The study aims to compare the HEL method with the response-spectrum method, using synthetic records to maintain consistency while accurately reflecting the facility’s design basis. This approach enables direct evaluation of the HEL method in fluid–structure interaction analysis.

Following the ground motion acceleration synthesis methodology (described in “Method”), in this study, we develop the intensity function using the Saragoni & Hart envelope [72], which characterizes seismic wave maximum amplitude over time. This generates representative seismic input within the I(t) envelope, as shown in Fig. 10. The time-acceleration curves for SL-1 scenario show horizontal peak ground acceleration of 0.060 g and vertical peak acceleration of 0.045 g. The generated motions maintain an average error below 10%, demonstrating high fidelity in synthetic accelerogram representation.

Fig. 10Full-screen View

(Color online) Synthetic ground motion and corresponding response spectra. (a) Acceleration time history in the horizontal direction; (b) Corresponding response-spectrum envelope; (c) Acceleration time history in the vertical direction; (d) Corresponding response-spectrum envelope

The structural form in this section aligns with the structure described in Table 4. The primary analysis condition includes 75% of the tank’s internal height specified in Table 5, corresponding to the L4 liquid level. This condition examines the structural performance under significant load. The structure is also analyzed under three other liquid-level conditions to ensure comprehensive understanding across operational scenarios, ensuring robust conclusions.

Results and discussion

The analysis of fluid sloshing and structural response under seismic excitation provides insights into the dynamic interaction between liquid and tank structure. Figure 11 presents curves showing the temporal variation in maximum sloshing height of liquid within the water tank, analyzed using the HEL method. This graph captures the liquid’s dynamic response under seismic excitation. The results indicate that at 4 s, the L4 liquid sloshing height reaches its maximum value, corresponding to the intensity of the synthetic ground motion in Fig. 10. In contrast, for the L2 and L3 liquid levels, the maximum values occur around 8 s. This delay may be due to liquid mass inertia and damping effects within the tank structure. The tank of the L1 liquid level is empty, and therefore, no sloshing height exists. These findings demonstrate seismic loading’s impact on fluid behavior. The variations in sloshing heights across different liquid levels emphasize the need for assessing fluid–structure interactions during seismic events for optimizing tank design and structural integrity.

Fig. 11Full-screen View

(Color online) Maximum sloshing wave heights over time under seismic loading based on the HEL method (Smoothed results)

A series of snapshots of fluid sloshing within the tank under seismic excitation at two liquid levels were captured at intervals, as shown in Fig. 12. The minimum value for explicit stress in the cloud chart was set at 1 MPa to enhance visibility of fluid sloshing effects on structural stresses. These images show stress distribution changes on the tank structure over time. The circumferential stress is uniform below the liquid surface due to hydrostatic pressure. Above the liquid surface, liquid sloshing causes notable changes in structural stress distribution through fluid-wall interaction.

Fig. 12Full-screen View

(Color online) Temporal evolution of structural stress distribution in seismic time-history analysis. (a) Level L3; (b) Level L4

The stress and displacement responses are extracted at 0.01 s intervals, as shown in Fig. 13, using an upper envelope method to standardize structural response curves. Figure 13b and Fig. 13c show higher stress levels at elevated liquid volumes compared to reduced volumes. Using the HEL method, stress values at intermediate to high liquid levels (L2-L4) are approximately double those from response-spectrum analysis. In scenarios without liquid, both methods yield comparable stress magnitudes. This difference occurs due to two factors: increased liquid sloshing amplitude at higher levels enhances fluid inertial forces (shown in Fig. 13d, Fig. 13e, and Fig. 13f), and hydrostatic pressure intensifies at higher fluid levels. While stress levels near the liquid surface do not reach peak levels seen elsewhere, they remain high relative to maximum stress points.

Fig. 13Full-screen View

(Color online) Temporal variation of the maximum stress and displacement of the structure at different liquid levels. (a) Stress response envelope; (b) Stress at the liquid-level side point; (c) Stress at the dynamic peak point; (d) Displacement at the top observation point; (e) Displacement at the liquid-level side point; (f) Displacement at the dynamic peak point. Dotted lines indicate results from the response-spectrum method

As the seismic event progresses, the relative displacement of the structure peaks between 8 s to 10 s, despite reduced acceleration. This indicates a nonlinear and delayed displacement response, showing cumulative effects under seismic influences. The structural response reflects a complex interplay of factors accumulating over time. This is supported by maximum values observed for L2 and L3 liquid-level sloshing, as shown in Fig. 11. The alignment of these maxima with structural displacement peaks demonstrates fluid dynamics’ influence on structural behavior during seismic events. The HEL analysis shows that the average relative displacement at three observation points equals the response-spectrum analysis result, while the peak-to-peak value is approximately twice that of the response spectrum analysis result. This finding in Fig. 13d to Fig. 13f aligns with stress results in Fig. 13b and Fig. 13c for the liquid-containing tank, highlighting the importance of considering both instantaneous and cumulative responses in structural analysis.

Further analysis shows that the structural acceleration response peaks at 2 s (refer to acceleration time-history curve and I(t) envelope in Fig. 10) represent the maximum instantaneous dynamic load from seismic activity. Intriguingly, as shown in Fig. 13d, Fig. 13e, and Fig. 13f, the relative displacement is minimized at this point, indicating that the maximum acceleration occurs at lower displacement values, likely due to phase lag between acceleration and displacement. This minimal displacement suggests a dynamic counter-adjustment within the structure, explored through the HEL method across varying liquid levels.