Numerical Investigation of Pool Boiling Under Ocean Conditions with Lattice Boltzmann Simulation. Part Ⅰ: Heaving Condition

Pool boiling is the heat-transfer mechanism of many heat exchangers inside ocean nuclear power plants working under the complex marine circumstances. Also, ocean conditions will create a new acceleration field other than gravity for the fluid, which induces some unique thermal–hydraulic characteristics. In this study, pool boiling under heaving conditions is numerically simulated using multiple relaxation time phase change lattice Boltzmann method. Firstly, the simulated results under static condition have been validated with recognized empirical equations, such as Rohsenow’s correlation at nucleate boiling, Zuber’s model, and Kandlikar’s model about critical heat flux (CHF). Then, pool boiling patterns, the boiling curve of time-averaged heat flux, transient heat flux, and heaving effects on different pool boiling regions are investigated. The results show that pool boiling curves of time-averaged heat flux between heaving conditions and static conditions with middle superheat degrees are similar. Heat transfer under heaving conditions at low superheat is somewhat enhanced, and it is weakened at high superheat, which leads to a slightly smaller critical heat flux with larger superheat compared with that under static conditions. Moreover, distinct fluctuation of the transient heat flux of pool boiling under heaving conditions is found for all boiling regimes. Furthermore, the heaving condition shows both positive and negative effects on pool boiling heat transfer at high-gravity and low-gravity regions, respectively. Besides, both the larger heaving height and shorter period time bring out more violent heaving motion and make a greater impact on pool boiling heat transfer.


INTRODUCTION
Marine applications, such as ocean nuclear power plants (ONPPs), power propulsion of ships, offshore small modular reactors, floating liquefied natural gas (FLNG) platforms, and seawater desalination apparatuses, have been developed by a great leap forward in these years with the largescale exploitation and utilization of ocean resources (Yan, 2017;Wang et al., 2019;Tian et al., 2020). Unlike the land-based condition, ocean applications would be affected by marine motions, which are caused by wind forces, tides, and ocean currents. The marine motion would induce additional inertial forces on the system and thus change the thermal-hydraulic characteristics of the energy and power plants. Also, the additional inertia forces introduced by rotational movement could be divided into normal inertial force, tangential inertial force, and additional Coriolis force (Yan, 2017) as follows: F in − m(a n + a t + a co ), where F in is the additional inertia force, a n is the normal acceleration, a t is the tangential acceleration, and a co is the Coriolis acceleration; also, translational motions will bring out additional forces along the corresponding direction. For research convenience, marine motions are generally decomposed to three translational motions (heaving, swaying, and surging) and three rotational motions (rolling, pitching, and yawing) (Wang et al., 2019), which are illustrated in Figure 1. Among these, heaving and rolling conditions are most commonly encountered, which will induce one kind of inertial force parallel to gravity and three kinds of inertial forces, respectively, and their effects on the work fluid inside the marine applications would produce unique thermal-hydraulic characteristics, which will be studied in this study series with Part Ⅰ focused on heaving conditions and Part Ⅱ (Zou et al., 2021) concentrated on rolling conditions. In recent decades, studies about thermal-hydraulic characteristics of marine applications under heaving conditions have been focused on flow circulation or flow boiling , which were mainly applied to the marine reactors (O'Rourke, 2010). These investigations have included critical heat flux (CHF), heat-transfer coefficient (HTC), bubble behaviors, flow patterns, flow instability, and pressure drop. Among these studies, CHF, HTC, and bubble behaviors could be referenced for the research of pool boiling mechanisms under heaving conditions. Most studies have shown that heaving motion could decrease CHF of flow boiling. In early 1966, Isshiki (1966) found that CHF decreased linearly with the increase of heaving acceleration, and Otsuji and Kurosawa (1982) got a similar result that heaving motion reduced CHF because of long bubble slugs causing the dry-out of the liquid film. Then, Ishida et al. (1990) pointed out that the CHF ratio was proportional to the 1/4th power of the gravitational acceleration ratio in their model and verification experiment of flow boiling under heaving conditions. Hwang et al. (2011) simulated CHF under heaving conditions using commercial codes, and the results showed that the CHF was much lower than that under static conditions because of gravity acceleration oscillation and the decrease of mass flow rate. Gui et al. (2020) theoretically investigated CHF characteristics in a narrow rectangular channel by considering the additional forces due to heaving and rolling motions, and it was shown that the CHF decreased because of inlet flow oscillation. Liu et al. (2018) developed a combining model investigating the CHF and the flow instability in ocean motion and found that heaving motion had little effect on CHF.
As for HTC and bubble behaviors of boiling under ocean conditions, almost all works were focused on the rolling condition, and thus, investigations about effects of heaving conditions were scarce. Xi et al. (2015) theoretically computed that the time-averaged HTC was reduced under heaving conditions compared to that under a static state. Recently, Hong et al. (Ishida et al., 1990) acquired the bubble departure size in forced convective flow under static and heaving conditions, and they found that the bubble departure diameter increased by increasing the heaving acceleration. The additional inertial force caused by heaving motion has been always parallel to gravity force, which would lead to the bubbles withstanding the high gravity and low gravity states. The study of Macro et al. (Di Marco and Grassi, 2002) about pool boiling under different gravitational accelerations showed that bubbles became smaller at high gravity, and they grew bigger at low gravity. Similar results could also be found in the simulation by Ma et al. (2017) and experiments by Raj et al. (2010). Besides, the gravity influences on CHF have also been studied in previous research studies (Raj et al., 2010;Ma et al., 2017;Fang et al., 2019), and the results showed that higher gravity within a certain range could increase CHF.
Relatively few studies have been published on pool boiling under ocean conditions, especially the investigations under heaving motion. Actually, many marine heat exchangers adopted pool boiling (Yu, 2020), and their unique thermal-hydraulic characteristics under complex ocean conditions would need to be uncovered. The lattice Boltzmann method is a mesoscopic numerical method potential on multiphase flow and phase transition (Huang et al., 2019;Huang et al., 2021). Recently, simulations using the lattice Boltzmann method (LBM) about pool boiling have been a hot topic (Cheng et al., 2014;Li et al., 2016), and it is very effective to incorporate the inertial forces caused by ocean conditions. Among multiphase models of LBM, pseudopotential LBM and phase field LBM have been applied to phase change heat transfer (boiling, evaporation, and condensation). Zhang and Chen (2003) took the first LBM simulation of nucleate boiling using the pseudo-potential model with a new scheme of the force term. Later, Hazi and Markus (2009) developed a pseudo-potential LBM phase change model based on a local balance law for entropy by introducing a source term to the temperature equation, which applied the equation of state (EOS) to drive spontaneous phase change without artificial hypothesis, and this treatment has been intensively November 2021 | Volume 9 | Article 771758 2 adopted. Subsequently, Gong and Cheng (2012) proposed an improved phase change source term and obtained the whole pool boiling curve (Gong and Cheng, 2015). Besides, Li et al. (2015) developed a hybrid thermal pseudo-potential LBM model by taking a finite-difference scheme to compute the temperature field and investigate the effects of surface wettability on boiling heat transfer. On the other hand, Dong et al. (2009), Safari et al. (2013, and Tanaka et al. (2011) devised their phase change model from the category of phase field LBM, respectively. With the help of these theoretical models' development, LBM simulations of phase change heat transfer have ushered in the explosive growth. During these studies, a pseudo-potential model has been most widely used because of its high efficiency, and the whole boiling curve including natural convection, nucleate boiling, CHF, and film boiling was investigated, together with the bubble behaviors of nucleation, growth, coalescence, and departure (Gong and Cheng, 2015;Li et al., 2015). This method could also be used to simulate boiling heat transfer under marine conditions by taking into account the inertial forces.
In this study, pseudo-potential LBM is adopted to study pool boiling under heaving conditions. In the Model Description Section, we discuss briefly the model used in this study and the strategy of including the inertia force. Then, the computation setup and boundary conditions are given, and the model is validated with recognized empirical equations of pool boiling under static conditions in the Computation Setup, Boundary Conditions, and Model Validation Section. Moreover, results and discussion are implemented in the Results and Discussion Section, where the mechanism of pool boiling under heaving conditions is shown, with conclusions drawn in the Concluding Remarks Section. Besides, it is noteworthy that the methodology and results in this study serve a reference for the investigation of pool boiling under rolling conditions in Part Ⅱ of this study series (Zou et al., 2021), which provides a comprehensive study of pool boing under complex marine conditions together.

MRT Pseudo-potential Multiphase LBM
The evolution equation of LBM with the MRT collision operator is given by the following (D'Humières et al., 2002): where f i (x, t) is the discrete density distribution function at position x and time t with particle velocity e i , Λ M −1 SM is the collision matrix, S is a non-negative diagonal matrix relevant to the relaxation time given by S diag , and M is an orthogonal transformation matrix (D'Humières et al., 2002). In this study, S diag (1.0, 1.1, 1.2, 1.0, 1.0, 1.0, 1.0, ). Γ M −1 (I-S/2)M is the collision matrix for forces. The relaxation time τ v is related to fluid kinematic viscosity with τ υ υ/(c s 2 δ t ) + 0.5, and the density distribution function f is projected into the moment space with M, m = Mf, and m eq Mf eq , in which m eq is defined as follows: The detailed collision and streaming processes in Eq. 2 are executed as The last term Fin Eq. 2 is the discrete force implemented into LBM, and the method proposed by Guo et al. (2002) is adopted in MRT as where ω i is the weighting coefficient with ω 0 4/9, ω 1-4 1/9, and ω 5-8 1/36; c s is lattice sound speed given by c s 2 c 2 /3 in the D2Q9 model with c δ x /δ t being the lattice speed; and v is the real fluid velocity. The discrete lattice velocity vectors e i are given by the following: Moreover, F is the total force, including inter-particle force F int , fluid-solid interaction force F s , gravitational force F g , and inertial force F h caused by heaving motion.
The inter-particle interaction force F int is given as where β 1.16 is a weighting factor depending on the equation of state used; it is chosen with simulation tests on Maxwell reconstruction in the equation of state . G (x + e i δt) denotes the interaction strength of the fluid, which is defined as φ(x) is the effective density computed as where c 0 6 in the D2Q9 model and p is the pressure calculated with the equation of state. Here, the Peng-Robinson equation of state (PR EOS) is adopted, which is given as where a 0.45724R 2 T 2 c /p c , b 0.0778RT c /p c , T c is the critical temperature, p c is the critical pressure, ρ c is the critical density, R is the gas constant, and w is the acentric factor. In this simulation, a 2/49, b 2/21, R 1, and w 0.344 are chosen.
The fluid-solid interaction force F s is calculated as where G s is the parameter denoting the fluid-solid interaction strength, and s(x) is a binary function (with s(x) 1 for solid and s(x) 0 for fluid).
The gravitational force is defined as where g is the gravitational acceleration. The macroscopic density ρ and velocity u of the fluid are given as The real fluid velocity v is computed as

Energy Equation Model
The evolution equation of the temperature distribution function is given by the following (Gong and Cheng, 2012): where τ T is the relaxation time related to thermal diffusivity α with τ T α/(c s 2 δ t ) + 0.5 and g i eq (x, t) is the equilibrium distribution function for temperature, which is given as Here, we can take v 0 to simulate heat conduction inside the solid.
The source term ϕ resulting from phase change is defined as follows (Zhang and Cheng, 2017): The temperature is calculated as

Inertial Force Under Heaving Motion
Heaving motions are shown in Figure 2, and extra inertial acceleration can be figured out with the change of heaving displacement h, which has sine relationship with time, where h m is the maximum heaving height, t p is the heaving period, and a h is the heaving acceleration. h m 1900δ x and t p 100000δ t are adopted in this study without otherwise specified, and the maximum value of a h is 0.5g. It should be noted that the simulation is based on the inertial system under heaving motion, and the direction of additional inertial force is opposite to the heaving acceleration, which is determined as

COMPUTATION SETUP, BOUNDARY CONDITIONS, AND MODEL VALIDATION
In this section, we will first introduce the computation setup and boundary conditions. Then, the model is validated with the Laplace law and heat-transfer curve with recognized empirical equations of pool boiling under static conditions.

Computation Setup and Boundary Conditions
The 2D computation domain shown in Figure 3 is an 800 × 1,000 lattice region (corresponding to L 1 × L 2 16.5l 0 × 20.6l 0 , where l 0 is the capillary length), which was found to be able to ensure accuracy through grid dependency checking. The liquid (below L h 400) is heated in an adiabatic solid container with a heater (L x × L y 600 × 100) placed downside, and the top of the container is open. Some researchers have studied the influence of liquid height on pool boiling and found that the liquid level at 5l 0 is enough for studying HTC and CHF (Pioro et al., 2004), and the liquid level of 6.2l 0 is adopted in this study.
The computation is initialized with the saturated liquid and vapor at T 0 0.8T c with property parameters as given in Table 1. The kinematic viscosity ] and thermal diffusivity α will directly influence relation time, so the values are obtained with some simulation tests based on a reasonable Pr ]/α. Then, the dynamic viscosity is calculated with μ ραPr. On the other hand, when the thermal conductivity ratio k l /k v is determined, c p,l /c p,v can be calculated by α k/ρc p . Latent heat of vaporization h fg 0.4392 is calculated with the method proposed by Liu and Cheng (2013). The physical parameters ζ (including υ, α, c p , and c v ) of the fluid in the interface regime are computed by linear interpolation as The convective boundary (Lou et al., 2013) is applied on the upper boundary, and a constant wall temperature T sw T 0 is adopted on the left and right sides. A heater is placed in the center of the bottom, with physical parameters given as (ρc p ) s 23.6 and November 2021 | Volume 9 | Article 771758 5 α s 4, corresponding to k s /k l 80 and k s /k v 1,360. The heater temperature downside is fixed at T b with the other position being adiabatic. Besides, the conjugate thermal boundary is adopted on the heater surface as follows (Li et al., 2014):   where i is the opposite direction to i. The thermal mass ratio of the solid and fluid is defined as ξ (ρc p ) s /(ρc p ) f .

Validation of the Model With the Laplace Law
Consider first the problem of the formation of a stationary bubble with a radius r in a saturated liquid. The pressure jump Δp across the bubble interface depends on the Laplace law as where p in and p out represent the pressure inside and outside the bubble, respectively, and σ is the surface tension. Figure 4A shows the coexistence of a stationary bubble (blue color) with radius r 40 surrounded by its saturated liquid (red color) at T s 0.8T c in a 200 × 200 domain with period boundaries on all borders. As presented in Figure 4B, by changing the initial radius of the bubble, the simulated pressure jump Δp is proportional to 1/r, which is in good agreement with the Laplace law, and the surface tension is obtained as σ 0.2454.

Simulation of Pool Boiling Under Static Conditions and Validated With Empirical Equations
To analyze bubble scale process, characteristic parameters are chosen as follows (Son et al., 1999;Gong and Cheng, 2017): where l 0 is the capillary length and u 0 and t 0 are, respectively, the characteristic velocity and time.
The corresponding dimensionless parameters are then defined as l* l/l 0 , t* t/t 0 , and the dimensionless wall superheat is defined as the Jakob number, Ja c p,l (T w − T sat )/h fg .
Moreover, the dimensionless time-averaged heat flux q′ is defined as follows (Ma and Cheng, 2019):  where L x is the heater length, Δt is the evolution time of pool boiling, and q bo is the reference boiling heat flux with q bo μ l h fg /l 0 . The boiling curves with time-averaged heat flux q′ calculated during the heaving period t p (40,000-140,000δ t ) is shown in Figure 5A. For the nucleate boiling region, the Rohsenow's correlation (Rohsenow, 1951) is used for comparison, where c sf is an experimental coefficient. C l l 0 real /l 0 LB is the length scale factor between LB simulation and reality, and C l 3.89 × 10 −5 is calculated in this study. The red dashed line in Figure 5A is the curve of Rohsenow's model with c sf 0.0059, and c sf is between 0.00225 and 0.027 for fluid-heating surface combinations in the literature (Holman, 2010). The simulation result fits well with Rohsenow's correlation.
The simulated critical heat flux (CHF) is compared with the CHF model proposed by Zuber (Zuber, 1959), which is given by the following: where K 0.12 to 0.157 and K 0.13 is adopted for comparing (shown as the horizontal dashed line on the top in Figure 5A). The simulated CHF is 8.0% higher than Zuber's model. Moreover, it should be mentioned that in Kandlikar's CHF model (Kandlikar, 2001), K is calculated from the contact angle as In this study, the contact angle of the surface is 45°, and the simulated CHF agrees very well with Kandlikar's CHF model. Pool boiling patterns with typical bubbles' evolution during different boiling regimes under static conditions are presented in Figure 5B. It is shown that there is no bubble nucleation on the heater in the natural convection regime where the heat flux is low. As the wall superheat (Ja) increases, bubbles start to nucleate on the wall and grow or merge to a size big enough to departure from the surface. The heat flux in the nucleate boiling regime increases with a higher slope as shown in Figure 5A. Then, vigorous bubbles' coalescence causes the occurrence of local dry-out, and the heat flux reaches the maximum value at Ja 0.186. After CHF, it enters into the transition boiling regime where a larger vapor film patch covers the heater, and the boiling heat flux decreases by increasing wall superheat. Furthermore, as the wall superheat goes to Ja 0.3, the boiling pattern changes to film boiling where a stable vapor film with two vapor column forms on the surface and restrains the liquid supplementing to the hot surface.

RESULTS AND DISCUSSION
In this section, the pool boiling heat-transfer curve of timeaveraged heat flux between heaving conditions and static conditions is first computed. Then, heat-transfer characteristics of pool boiling under heaving conditions are given, including pool boiling patterns at different regimes, the fluctuation of transient heat flux, and the effects of heaving effects on pool boiling at different boiling regions.

Pool Boiling Patterns Under Heaving Conditions
After the above verification of the model with empirical equations under static conditions, pool boiling under heaving conditions is simulated in this section with a heaving period of 100000δ t with a maximum additional inertial acceleration of 0.5g 0 . Representative pool boiling patterns during one stable heaving period (t 3 × 10 5 δ t ∼ 4 × 10 5 δ t ) are presented in Figure 6. At the nucleate boiling regime as shown in Figure 6A, bubbles are observed to have a normal size at the beginning; then, they become smaller at the half supergravity period from 0t p to 0.5t p . Then, it enters the half low-gravity period from 0.5t p to 1.0t p , and bubbles return to the normal size and keep growing bigger for low gravity. According to the analysis given by Friz (1935), the departure diameter of bubbles is proportional to g −0.5 , which is also validated in other simulations (Gong and Cheng, 2012;Fang et al., 2017;Li et al., 2017). Therefore, high gravity leads to small bubbles, and low gravity brings out large bubbles, which has been experimentally observed (Di Marco and Grassi, 2002;Raj et al., 2010), and it is consistent with this simulation result.
The results for pool boiling at CHF under heaving conditions as presented in Figure 6B are similar to that at the nucleate boiling regime. More bubbles exist because of higher superheat, and the bubbles also undergo the cycles of smaller sizes at the half supergravity period and bigger sizes at the other half low-gravity period. As for film boiling regime in Figure 6C, the change of bubble size is also alike, and it needs to be noticed that the shape of the vapor film is much influenced by the heaving motion.

Comparison of the Boiling Curve Between Heaving Conditions and Static Conditions
In this section, we take heaving period t p 10 × 10 4 δ t and maximum heaving height h m 1900δ x . The boiling curve with the time-averaged heat flux during the heaving period t p (40,000-140,000δ t ) is presented in Figure 7, which includes the whole boiling regimes of natural convection, nucleate boiling, CHF, and film boiling. It is shown that the heaving effects have a very slight effect on the time-averaged heat flux. At convection and film boiling regions, the time-averaged heat flux curves between heaving and static conditions may almost coincide with each other because of the perfect periodicity of heaving motion. As for the nucleate boiling region and CHF under heaving motion, heat transfer at low superheat is somewhat enhanced, and it is weakened at high superheat, which leads to slightly smaller CHF with larger superheat compared with that under static conditions. This result is attributed to the combined effects of heaving motion on bubble detachment and coalescence. Detachment of isolated bubbles at low superheat is improved by heaving motion, which increases the time-averaged heat flux. On the other hand, there are times when low gravity appears under heaving conditions, the bubbles are bigger with the lower buoyancy, and more bubbles stay and merge on the surface. Also, coalescence of bubbles at high superheat is drastic under heaving conditions, and rewetting of the liquid to the heater surface is hindered, which leads to a somewhat lower CHF. The effects of heaving motion on the whole boiling curve under different boiling regimes will be analyzed as follows.

Transient Heat Flux at the Convection Regime Under Heaving Conditions
It is shown in Figure 7 that heaving motion has a slight effect on the time-average heat flux of pool boiling because of its periodic nature, and this phenomenon will be analyzed based on the variance of transient heat flux q t ′, which is defined as The transient heat flux at the convection regime is given in Figure 8. It can be found that the transient heat flux fluctuates with time periodically under heaving conditions. The gray stippled line is plotted with heaving height h 0, and the gray dashed line marks the high-gravity region. By comparing the curve of heaving height h and transient heat flux under heaving conditions, we find that after it comes to the high-gravity stage, the heat flux rises, and when it enters the low-gravity region, the transient heat flux starts to decrease. The peaks of transient heat flux always appear after it comes to the low-gravity region, and  November 2021 | Volume 9 | Article 771758 9 the valley will be found after it becomes high gravity. In other words, the heat flux will be minimum shortly after the containers fall to the average height (h 0) and will be maximum after the containers rise to the average height. The phase offset between transient heat flux at convection and heaving motion is caused by the inertia of liquid flow. On the other hand, transient heat flux under static conditions will be no obvious pattern, and its variance range is much lower.
In order to explain the effects of heaving motion on transient heat flux, graphs of velocity vector and pressure field are given in Figure 9, which shows the moments at h 0. From t 0t p to 0.5t p , it undergoes a high-gravity process and the convection velocity increases at increasing pressure field, and thus, the transient heat flux rises. On the other hand, during the half period at low gravity from t 0.5t p to 1.0t p , the convection velocity decreases at decreasing pressure field, and thus, the transient heat flux goes down. Since the convection cycle is caused by gravity and density difference, the transient heat flux increases at the half high-gravity period because of convection promotion and vice versa at the half low-gravity period.

Transient Heat Flux at the Nucleate Boiling Regime Under Heaving Conditions
The transient heat flux at the nucleate boiling regime is plotted in Figure 10. The gray dashed line marked the high-gravity region. At this regime, the heat flux is mainly influenced by bubble nucleation and liquid convection. At the high-gravity region, the liquid pressure is high, which leads to difficult bubble nucleation. On the other hand, the buoyancy is large and the liquid convection is strong. It is interesting to observe that the variation of transient heat flux is a little more hysteretic than that at the convection region. The peaks of transient heat flux are always in the low-gravity region, shortly after the system rises to  From Figure 11, it can be concluded that large gravity leads to high pressure in the liquid, and the nucleation is impeded, so the heat flux is low at the high-gravity region. On the contrary, more bubbles and bigger bubbles are generated at the low-gravity region, corresponding to a high heat flux. Because the nucleation and growth of bubbles need time, the peak will be later than the moment with minimum gravity and finally will be found to be just before the containers rise to the average height (h 0 for here).

Transient Heat Flux at CHF Under Heaving Conditions
The transient heat flux at CHF fluctuates violently under static conditions, as shown in Figure 12A, and it is noisy. And as seen in Figure 12B, the transient heat flux at CHF under heaving condition is different from that in nucleate boiling regime, because the peaks coincide with the high-gravity region. Analyzing the bubbles' behavior in Figure 13, it is found that the superheat at CHF is always very high and bubble nucleation is strong all the time at this moment. On the other hand, since the bubbles are small at the high-gravity region, the water rewetting is easier and less area is occupied by vapor, so it corresponds to the peak of heat flux at the high-gravity region. Besides, the pressure gradient is larger with higher gravity, and the bubble's buoyancy is improved, so the heat flux is higher with easier bubble departure. On the contrary, when the gravity is low, more area is occupied by vapor and the bubble's departure is also influenced, and the heat flux is low.

Transient Heat Flux at the Film Boiling Regime Under Heaving Conditions
At the film boiling regime, the fluctuation of transient heat flux is more regular, as shown in Figure 14, since the bubble behavior is periodic at this moment. The peaks represent the bubble's departure from the vapor film, and since bubbles grow very regularly at film boiling, there exist many small peaks in the transient heat flux. Moreover, the fluctuation of transient heat flux under heaving conditions is larger at the low-gravity region and smaller at the high-gravity region. The reason is different-size bubbles' departure will leave vapor films of different thicknesses.

Effects of Heaving Parameters
The influence of heaving parameters on nucleate boiling is studied, and the transient heat flux curves are given in Figure 15, where   Figures 15A-D), and these on the same column have the same heaving period ( Figures 15A-D). It can be found that both larger heaving height and shorter period time bring out more violent heaving motion, and the fluctuation amplitude will be larger. Besides, the fluctuation period of heat flux is the same with the heaving height, which can also be reflected by comparing Figures 15A,C with Figures 15B,D. It needs to be noticed that the valley may be very deep for short heaving period time as seen in Figure 15B, because more and larger bubbles burst to nucleate and grow, and then vapor films form on the heater surface temporarily.

CONCLUDING REMARKS
In this study, lattice Boltzmann simulations of pool boiling under heaving conditions are performed. First, the simulated boiling curve under static conditions has been validated with Rohsenow's correlation for nucleate boiling, Zuber's model, and Kandlikar's model for CHF. Then, heat-transfer characteristics and bubble behaviors of pool boiling under heaving conditions are investigated, and the differences with that under static conditions have been analyzed. Moreover, heaving effects on bubble behaviors and heat transfer at different pool boiling regions are studied. The following conclusions can be drawn from this study: 1) Pool boiling under heaving condition will go through periods of high gravity and low gravity, which leads to unique bubble behaviors of nucleation, growth, coalescence, and departure. Moreover, pool boiling heat transfer under heaving condition is somewhat enhanced at low superheat, and it is weakened at high superheat, which leads to slightly smaller CHF with larger superheat compared with that under static condition. 2) The variance of transient heat flux under heaving conditions is directly related with bubble behaviors, and its fluctuation is very distinct. The transient heat flux and heaving motion keep good synchronization for all pool boiling regimes (convection, nucleate boiling, CHF, and film boiling), and their phase offset is different according to pool boiling intensity and bubble behaviors under large-gravity or lowgravity stages. On the contrary, the transient heat flux of pool boiling under static conditions is relatively chaotic. 3) Large gravity under heaving motion leads to high pressure in the liquid, and bubble nucleation is impeded. On the contrary, more bubbles and bigger bubbles are generated at the low-gravity region, corresponding to a high heat flux. 4) Both a larger heaving height and a shorter period time generate more violent heaving motion, and the fluctuation amplitude of transient heat flux for pool boiling will be larger. Besides, the valley of transient heat flux may be very deep for a short heaving period time. Frontiers in Energy Research | www.frontiersin.org November 2021 | Volume 9 | Article 771758