Cavitating flows of organic fluid with thermodynamic effect in a diaphragm pump for organic Rankine cycle systems

Diaphragm pumps often experience cavitation and subsequent fluid flow oscillation when delivering an organic fluid in small/micro scale organic Rankine cycle (ORC). The cavitation behaviour of diaphragm pumps has rarely been investigated for organic fluids so far. Three-dimensional, unsteady cavitating flows of organic fluid R245fa in a diaphragm pump were simulated with ANSYS 2019R2 CFX in suction stroke in terms of the k -u turbulence model, the ZGB cavitation model, rigid body motion model for onedimensional motion of valve and moving mesh technique for the first time. The thermodynamic effect in cavitation of R245fa was considered. The vapour volume fraction threshold for cavitation inception was determined, and the cavitation inception and cavitation developed states were identified, and vortex production and entropy generation rate during cavitation were clarified. Cavitation inception emerges at the edge of the valve seat, then on the valve surface. With cavitating development, the pressure and force on the valve, valve opening, and velocity oscillate violently due to vapour bubble collapse cycles. Expansion cavitation and flow induced cavitation happen in sequence at different crank rotational angles. The maximum temperature depression is 0.549 K in the cases studied. The volume-integrated entropy generation rate in the valve chamber correlates to cavitation states. © 2021 The Author(s). Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).


Introduction
The organic Rankine cycle (ORC) is a Rankine cycle by using an organic fluid with low boiling point as working medium to generate power from lower temperature sources such as biomass combustion, industrial waste heat and geothermal heat, etc [1e5]. Mechanical feed pumps are commonly employed to deliver an organic fluid to the evaporator in an ORC system.
The feed pump in ORC systems is featured with high head (20e4000) kPa and low mass flow rate (0.01e10) kg/s [6]. Two sorts of pump such as rotodynamic pumps and positive displacement pumps are suitable for that purpose. The rotodynamic pump has one or more rotating impellers to raise liquid pressure and velocity continuously, but the positive displacement pump delivers a liquid with a periodical increasing and decreasing volume intermittently. At the same head and flow rate, the positive displacement pump is subject to a better efficiency but higher maintenance costing and more complicated pipe systems than the rotodynamic pump.
Compared with rotodynamic pumps, positive displacement pumps are more suitable for ORC systems in terms of exergy destruction rate, heat absorption rate, thermal and exergy efficiency. However, positive displacement pumps have a higher NPSHr (net positive suction head required) (24 kPa) [16,17] or variable subcooling (4.4e20) C [15e18], and will be in high risk of cavitation in ORC systems. The cavitation in a mechanical feed pump for organic fluids can induce the instability of ORC system, especially impair the evaporator performance. Properly determining NPSHr or subcooling of the pump in handling an organic fluid is one important issue in the ORC system design and operation. Significant attention should be devoted to the cavitation in positive displacement pumps.
Organic fluids such as R113, R1114 and R245fa are thermosensitive in between water and nitrogen, as shown in Fig. 1 in terms of the thermodynamic parameter S, which is defined as [22].
where L is latent heat of a liquid at temperature T l , y l is specific volume of the liquid, c pl is specific heat capacity of the liquid, y v is specific volume of the vapour, D l is thermal diffusivity of the liquid, D l ¼ l l =ðr l c pl Þ, l l is thermal conductivity of the liquid. The thermodynamic parameter S can be considered as the proportional constant of the relationship between temperature depression and vapour bubble radius growth rate [22]. At a fixed temperature depression, the larger the thermodynamic parameter, the smaller the growth rate. At a fixed growth rate, the larger the thermodynamic parameter, the higher the temperature depression. The larger the thermodynamic parameter, the more thermo-sensitive the fluid. The thermodynamic effect in cavitation should be considered in determination of the cavitation performance of a pump when it delivers an organic fluid in an ORC system. In that context, the cavitation performance of a diaphragm pump of G20-E model when delivering organic fluid R245fa was studied analytically [23]. A NPSHr correction method for thermodynamic effect in cavitation was proposed and the onedimensional (1D) motion of the suction valve of the pump was simulated numerically by using two mechanical models to identify the cavitation at 100 kPa and 141 kPa pump inlet pressures along with preliminary experiments. The net positive suction head available (NPSHa) at those inlet pressures was calculated and the cavitation safety margin was addressed, the subcooling for the NPSHr, NPSHa and safety margin was then determined. However, in 1D valve motion simulations, the pressure drop across the gap between the valve and the seat is predicted poorly in terms of empirical flow coefficient. The flow details through the gap cannot be resolved with simple 1D flow models.
Recently, the cavitation with thermodynamic effect has been increasingly studied with computational fluid dynamics (CFD) method for thermo-sensitive liquids flowing in nozzles, hydrofoils, and centrifugal pumps. Additionally, a few CFD simulations of cavitation in reciprocating pumps were conducted when handling cold water without thermodynamic effect [24e29]. However, there has not been any CFD investigation into cavitation performance and flow details in reciprocating pumps used in ORC systems when delivering organic fluids with thermodynamic effect in cavitation so far.
In the present article, three-dimensional (3D), unsteady cavitating flows of organic fluid R245fa with thermodynamic effect in a diaphragm pump for an ORC system are tackled by using CFD simulations in ANSYS 2019R2 CFX software. The work includes cavitation model constants calibration with organic fluid R114 in venturi flow systems in the literature, thermodynamic effect involvement in CFX, rigid body 1D motion modelling for the valve, moving mesh handling in CFX, cavitation inception identification and developed cavitation characterization, vortex production and entropy generation discussion. The work is original and can be meaningful to study on cavitating flows of organic fluid, valve design of diaphragm pump, and pump operation in ORC systems.

Mean pump performance
A positive displacement diaphragm pump of G20-E model was selected to feed organic fluid R245fa to the evaporator. The pump liquid end structure is sketched in Fig. 2a. The diaphragm is driven by the piston via the hydraulic oil in the left chamber of the diaphragm. When the piston moves to the left, the suction valve is opened and the discharge valve is closed, then the liquid is sucked into the chamber, and the pump is in suction stroke. As the piston moves to the right, the suction valve is closed but the discharge valve is opened, the pump is in discharge stroke.
The pump volume flow rate, inlet and outlet pressure and NPSHr vary with time. Based on the pump data charts, the mean pump flow rate, Q , and NPSHr are formulated as when pumping water at room temperature [23].

Instantaneous pump performance
The diaphragm is flexible rubber material, fluid-structure interaction analysis or visualization experiment is required to capture its motion and deformation characteristics [30]. Since our aim is to study the cavitation in the suction valve, the diaphragm actual motion is out of our scope. Here an idealized diaphragm motion and deformation is proposed as shown in Fig. 2b  During the motion, the diaphragm keeps in flat cone shape with two radii R 1 ¼ 5.75 mm, R 2 ¼ 15.75 mm and a variable cone height/ displacement x, x2½0; s. As x ¼ s=2, the diaphragm becomes a circular plane. The top surface with R 1 ¼ 5.75 mm is with the same speed as the piston, v. The cone surface stroke is reduced to zero at R 2 ¼ 15.75 mm from s at R 1 ¼ 5.75 mm linearly. Thus, the stroke profile is expressed by 8 > < > : At the beginning of suction stroke, i.e. t ¼ 0 or 4 ¼ Ut ¼ 0, and neglecting the 2nd order terms [31,32], the diaphragm displacement equation is written as The moving velocity of the diaphragm reads as The instantaneous flow rate in the suction pipe during 4 ¼ ½0; p is written as and the equation can be simplified into the following form [23].
The mean flow rate in the suction pipe in the period of 4 ¼ ½0; 2p is given by [23].
If R 1 and R 2 , s and U are known, the instantaneous and mean flow rates can be estimated with Eqs. (7) and (8), respectively. The geometrical parameters of the diaphragm are tabulated in Table 1. 3. CFD computational models

Cavitation model adopted
Numerical simulation of cavitating flows has appeared since 1960s. Cavitation models are essential for numerical simulations of such flows. Cavitation models can be classified into four categories: (1) interface tracking models, (2) homogeneous two-phase models, (3) multiscale models, and (4) stochastic models. The homogeneous two-phase cavitation model compromises between cavitation physics and computational costing and has found mostly extensive applications. In this category, the equation-of-state model, arbitrary mass transfer rate model, the Rayleigh-Plesset equationbased mass transfer rate model, and the nucleation cavitation models are four common models.
In the homogeneous two-phase cavitation model, the fluid is considered as a homogeneous, two-phase mixture with variable density. The continuity equation, momentum equations, energy equation, transport equation for number of bubbles or vapour/ liquid volume fraction are solved simultaneously in the model.
In [34], a Rayleigh-Plesset equation-based mass transfer rate model, called ZGB cavitation model, was proposed. In the model, the continuity equation of the liquid and vapour is written as where a is vapour volume fraction, which is defined as the ratio of the vapour volume in a given total volume, r v and r l are density of vapour and liquid, respectively; u i are time-averaged velocities in x i coordinates in a Cartesian coordinate system, x i are Cartesian coordinates, i ¼ 1,2,3; _ m is a source term for vapour/liquid phase change rate, and expressed by where a 0 is nucleation site volume fraction in a liquid, r is density of the mixture of vapour and liquid, r ¼ ar v þ ð1 À aÞr l , R b is vapour bubble radius, p is static pressure of the mixture of vapour and liquid. Since r l zr and let R b ¼ constant, the phase change rate is in the following form where F vap and F con are empirical model constants, R b0 is the radius of initial vapour bubbles or the radius of nucleation sites. In ANSYS 2019R2 CFX, which is finite volume method based CFD package for general propose of fluid flow and heat transfer, the cavitation model expressed by Eqs. (9) and (11) is implemented, with the default model constants are F vap ¼ 50, F con ¼ 0.01, R b0 ¼ 10 À6 m, a 0 ¼ 5 Â 10 À4 for isothermal cavitation of water [34], here a 0 ¼ , N b is number of bubbles per unit volume, and 14 should yield based on the known R b0 and a 0 . In the article, this model was adopted but with new model constants applicable to the organic fluid R245fa.

Involvement of thermodynamic effect
The ZGB cavitation model has been extended to include thermodynamics effect and there are two methods for involving the effect. The first method is vapour bubble radius growth rate correction, and the second method is temperature-dependent thermal property constant involvement. In the first method, the unsteady heat transfer equation [35e38] or steady heat transfer equation [39] is established in the interface between a spherical vapour bubble and the liquid to calculate the bubble radius growth rate correction due to liquid temperature depression. Also, the thermophysical and transport property constants are temperaturedependent, and the energy equation is activated. The model constants F vap and F con should be calibrated against experimental data of cavitating flows with thermodynamic effect. In the second method, the thermal property constants are temperaturedependent, and the energy equation is activated for temperature field along with calibrated model constants F vap and F con [40e46].
Cavitation regimes such as inertia-control, heat-diffusion control and intermediate states described in [47] are not concerned in the bubble radius growth rate correction methods [35e39], and they are subject to confused physical concepts. Thus, the second method is adopted here.
To include the thermodynamic effect, firstly, the thermal property constants of organic fluid R245fa should be variable with temperature. Based on the scattered thermophysical and transport property constants generated by REFPROP, which is a program, developed by the National Institute of Standards and Technology (NIST) to calculates the thermophysical and transport properties of industrially important fluids and their mixtures [48]. The R245fa liquid and vapour densities r l , r v , specific heat capacities c pl , c pv , dynamic viscosities m l , m v , thermal conductivities l l , l v , and vapour pressure p v are best fitted by the following expressions in terms of liquid bulk temperature T lb or local temperature T l ranged in 180 K and 420 K where ðr l =r v Þ max is the maximum density ratio of liquid R245fa to its vapour, even though ðr l =r v Þ max is not property constant, it is required by ANSYS CFX when the ZGB cavitation model is launched. Our experience witnessed that if all the property constants in Eq. (12) are expressed by local liquid temperature T l , the numerical simulation of a cavitating flow is unstable, or even cashed. Therefore, expect the vapour pressure, the rest property constants remain as function of bulk or far field liquid temperature T lb only.
To facilitate the interpolation of specific enthalpy and entropy in ANSYS CFX in terms of temperature, the reference states both liquid R245fa and its vapour are listed in Table 2. In ANSYS CFX, the specific enthalpy at others than the reference states are calculated by using specific heat capacities, c pl , c pv and local liquid temperature T l .

Calibration of the cavitation model constants
The thermodynamic parameter S of water, nitrogen, organic fluids R245fa, R113, and R114, is shown in Fig. 1 in terms of liquid temperature T l . The S value of a liquid specifies the extent of the thermodynamic effect in cavitation. The extent of the thermodynamic effect of the liquid R245fa is in between nitrogen and water and comparable to R114 or R113, especially R114. This means that the model constants F vap and F con in the ZGB model Eq. (11) should be smaller than the default F vap ¼ 50 and F con ¼ 0.01 for 20 C water.
For water at around 90 C flowing over a hydrofoil, these constants are F vap ¼ 10 and F con ¼ 0.002 [35,40]. For liquid hydrogen flowing over a quarter calibre hydrofoil, F vap ¼ 112 and F con ¼ 0.006 [37]. For liquid nitrogen, the constants F vap ¼ 5 and F con ¼ 0.001 can result in good agreement with experimental observations in temperature and pressure profiles of the cavitating flows around a 2D quarter calibre hydrofoil [38,39,44,45]; for liquid air, the constants are F vap ¼ 0.2 and F con ¼ 0.012 [46].
Since the thermodynamic effect in cavitation for organic fluid R245fa is different from that for hot water, liquid hydrogen and nitrogen, as shown in Fig. 1, the model constants applicable to R245fa should be calibrated based on existing experimental data. Unfortunately, the cavitation experimental data on R245fa could not be found in the literature, but the experiment on R114 cavitating flows in a venturi was available. Because the S curve of R114 is very similar to R245fa, as shown in Fig. 1, the experimental data of well-developed cavitating flows of R114 in a venturi [49] was employed to calibrate the model constants F vap and F con . The detail of the calibration is elucidated in Appendix. It is shown that the two model constants depend on bulk liquid temperature and Reynolds number.
The diaphragm pump operates with liquid R245fa at 8.4 C bulk temperature only, but the Reynolds numbers at the pump inlet and valve gap demonstrate significant variations with crank rotational angle f as shown in Fig. 3 based on 1D mechanical model described in [22]. The curves in the figure suggest that the influence of Reynolds number on model constants F vap and F con should be taken into account. It is assumed that the correlations of two model constants F vap and F con with Reynolds number shown in Fig. A5b are applicable to cavitating flows of organic fluid R245fa in the suction Table 2 The thermophysical and transport property parameters of R245fa at reference condition in CFD simulations.

Item
Parameter Numerical figure m v (Pa s) chamber of the diaphragm pump. In terms of the Reynolds number at the pump inlet Re h , the correlations are expressed by where the Reynolds number at the pump inlet is defined as Eq. (13) has been coded with CFX expression language (CEL) in CFX-Pre to implement two model constants into the cavitation model.

Computational models, boundary and initial conditions
The liquid end of the diaphragm pump is composed of suction chamber, suction valve, spring, retainer, pumping chamber, discharge valve and discharge chamber. There is a liquid flow in the valve chamber and pumping chamber only when the pump is in the suction stroke. Since cavitation commonly emerges in the suction chamber during the suction stroke, the valve chamber, suction valve and pumping chamber are included in the computational model and the rest components are neglected, as shown in Fig. 4. Because the valve chamber, suction valve and pumping chamber are symmetrical about the plane through the centre lines of the suction and discharge cylindrical chambers, only half of the geometries are taken into account in CFD simulations, as illustrated in Fig. 5.
In the figure, a 50 mm long suction pipe is connected with the valve chamber inlet to help the chamber with a smooth liquid flow and easy boundary condition implementation. There are three fluid domains, namely suction pipe fluid domain, valve chamber fluid domain and pumping chamber fluid domain. One fluid-fluid interface is built between the suction pipe and the valve chamber, and the other between the valve chamber and the pumping chamber.
In these fluid domains, the liquid is treated as a homogeneous two-phase mixture flow of liquid R245fa and its vapour. The two phases are incompressible and share the same velocity and pressure. The mixture undergoes a turbulent flow in the fluid domains during the suction stroke and yields the unsteady Navier-Stokes equations and energy equation, k À u turbulence model, and vapour volume fraction transport equation Eq. (9) with the source term Eq. (11). These equations and their corresponding models can be found in [52] and are omitted here.
The suction valve is regarded as rigid body with mass m v and the valve 1D motion along the suction chamber axial direction is determined by the Newton's second law of motion. The forces acting on the valve include the spring force, the fluid drag force and the force due to the pressure difference across the valve. The first force is calculated by using the Hooke's law in terms of its stiffness k and total displacement h which equals to the valve displacement or opening or lift because the reference position of the valve is defined at the position where the valve just contacts with the seat, see Fig. 4. Since the suction chamber is horizontal when the diaphragm pump is in operation, the force due to the gravity acceleration is not considered. The motion equation of the valve reads [52].
where f p is the force due to the pressure difference across the valve, f R is the drag fore applied by the liquid on the valve, the last term in the right-hand side of Eq. (14) represents the spring force, h 0 is the pre-compressed displacement from the free state of spring, see  Table 1.
At the suction pipe inlet, the known constant static pressure, temperature and turbulence intensity are given, and the gradient of  fluid flow variables is zero. The total pressure is fluctuation during a suction stroke. It is difficult to define a total pressure during cavitation because the flow rate can vary with cavitation state. Thus, a constant mean inlet pressure is imposed at the inlet, and the flow rate is decided by the diaphragm motion prescribed and the cavitation state automatically. Zero gradients for flow variables mean that there a fully developed flow from the suction pipe to the pump inlet.
In the plane of symmetry of the model, a symmetrical boundary is held. The diaphragm is adiabatic and subject to a moving boundary condition. In that case, the axial displacement of the diaphragm has to be specified. Eq. (3) is included into the model by using CEL in CFX-Pre. The rest boundaries are adiabatic and no-slip smooth wall, and the scalable wall function implementation is imposed.
As the valve and the diaphragm are moving, a dynamic mesh is needed in both the suction chamber and the pumping chamber fluid domains. The displacement diffusion algorithm is employed to increase the mesh cell volumes near smaller volumes than the average ones in the fluid domains with model exponent of 0.75.
Initial conditions are required in simulations of the unsteady cavitating flow in the domains. In the three fluid domains, the axial velocity, liquid temperature, relative pressure and turbulence intensity are provided with 0.001 m/s, 8.4 C and 141 kPa or 89 kPa depending on operational conditions studied, and 5 %, respectively.

Mesh size and time-step independence
Two meshes, i.e., fair and fine meshes, were created and three time-steps were tried to determine the independence of simulated results on mesh size and time-step. The mesh size, mesh information, mesh metrics, time-steps adopted, and six cases are listed in Table 3. The fair mesh pattern and the meshes in the suction chamber and pumping chamber fluid domains are illustrated in Fig. 5.
The mesh pattern is unstructured and includes the elements such as wedges (56.6 %), tetrahedrons (39.0 %), hexahedrons (1.9 %) and pyramids (2.5 %) for the fair mesh. In the valve chamber and pumping chamber fluid domains the meshes are attached with 10layer and 8-layer boundary layer mesh near the walls. Since the gap between the seat and the valve is quite small and the pumping chamber is thinner in the axial direction than the radial direction,  Table 3.
The valve had to be subject to 0.1 mm initial/dead lift to the seat when the meshes were created or a CFD simulation was impossible. A few lifts smaller than 0.1 mm were attempted, however, all the corresponding CFD simulations crashed. Therefore, 0.1 mm is the minimum initial/dead lift for CFD simulations of cavitating flow in this diaphragm pump at 480 rpm.
Six CFD cases were launched to examine effects of both mesh size and time-step on the motion of the valve at p 1 ¼ 141 kPa inlet pressure and the corresponding results are shown in Fig. 6. That the Courant number is less than 2 is a necessary condition for convergence while solving time-dependent partial differential equations. In the cases (fair, 50) and (fine,50), the volume averaged Courant numbers are as large as 1.08 and 1.17, and the relationships of the valve lift, velocity, minimum pressure on the seat, minimum pressure on the valve, and force acting on the valve by the liquid against the crank rotating angle f show a significant oscillation compared with the rest four cases where the volume averaged Courant number is smaller than 0.60. In the rest four cases, the valve lift, velocity, and force acting on the valve exhibit consistent trends of variation, suggesting a mesh size and time-step independency achieved.
Two minimum pressures in the cases (fine, 100) and (fine, 200) are slightly lower than the cases (fair, 100) and (fair, 200), the case (fine, 100) should serve as the time-step and mesh size independent case. Hence, the mesh and time-step in the case (fine, 100) were applied to the CFD simulations here.

Cavitation inception
Six transient CFD simulations were performed at inlet pressures p 1 ¼ 141, 89, 88, 87, 86, and 85.2 kPa by employing the computational models stated in Section 3 and fine mesh in terms of 100 of number of time-steps to identify cavitation behaviour of the valve. The minimum liquid static pressure p min , minimum liquid temperature T lmin , and maximum temperature depression DT lmax (¼T lb -T lmin ) and the maximum vapour volume fraction a max on the seat and valve were observed to clarify cavitation inception and temperature depression in the valve chamber. The results at the crank rotational angle 4 ¼ 106.2 are summarised in Table 4. Since the numerical solution process is unstable and convergent solutions no longer exist at p 1 <85.2 kPa, thus, the inlet pressure has to be limited to 85.2 kPa. The minimum pressures at the seat and valve p min , the resultant fluid dynamic force on the valve f , valve instant velocity V, valve lift h, and instant flow rate through the gap between the seat and the valve q are plotted in Fig. 7 as a function of 4 at five inlet pressures.
In Table 4, the cavitation is triggered by reducing pump mean inlet pressure p 1 . On the seat, as soon as p 1 is reduced to 89 kPa from 141 kPa, the minimum pressure on the seat p min is below the local vapour pressure p v , indicating cavitation occurrence. When p 1 is reduced further to 85.2 kPa, the cavitation develops on the seat and ends up with the maximum temperature depression DT lmax ¼ 0.274 K and the maximum vapour volume fraction a max ¼ 0.358 on the seat. On the valve, however, p 1 is not lower than p v until p 1 ¼ 86 kPa. At p 1 ¼ 85.2 kPa, the cavitation on the valve is so bad that DT lmax ¼ 0.400 K and a max ¼ 0.489. This fact suggests that the cavitation occurs at first on the seat and then develops on the valve quickly and significantly.
In Fig. 7, the variables p min , f , V, h, and q are in pulsation when the rational angle is larger than 60 at p 1 86 kPa. This phenomenon is closely related to the vapour bubble collapse effects in cavitation. The oscillating instantaneous flow rate profile during cavitation has been observed by using PIV in a reciprocating pump [53].
Usually, the cavitation inception is defined with occurrence of a visible vapour bubble in the lowest pressure zone in a liquid. For example, when a vapour bubble with 1 mm streamwise length occurs in a liquid in either experimental observations or CFD simulations, the cavitation inception is considered to occur [54]. The other proposal suggests 1 ml vapour volume is used to define cavitation inception [55,56]. Also, 0.01 % vapour volume fraction threshold, which is 0.0001 as volume fraction value, is for cavitation inception [57]. This value seems too small to be a threshold. 1 ml vapour volume is too large to our case and 1 mm streamwise is determined difficultly from CFD simulation results. Hence, we keep Table 3 The mesh information and cases for examining the independence of CFD results on mesh size and time-step.

Mesh information
Time using of vapour volume fraction as the threshold as done in [57] but put forward 0.10 vapour volume fraction as the cavitation inception threshold. From Table 4, the inlet pressure for the cavitation inception at 480 rpm should be in 89e88 kPa. A linear interpolation was performed based on the vapour volume fractions at the two inlet pressures, the inlet pressure for 0.10 vapour volume fraction was 88.6 kPa. Then, a transient simulation of cavitating flow at p 1 ¼ 88.6 kPa was launched, and the corresponding results are listed in Table 4 as well. At this inlet pressure, a 0.1048 vapour volume fraction is achieved, declaring the cavitation inception occurrence.

Developed cavitation
Developed cavitation depends on crank rotational angle and inlet pressure. At a fixed rotational angle, in Table 4, with decreasing inlet pressure, the minimum pressure p min , vapour

Cavitating flow pattern
The contours of vapour volume fraction, iso-surfaces of vapour volume fraction of 0.1, mixture static pressure, velocity and temperature at p 1 ¼ 85.2 kPa are illustrated in Fig. 9 at 4 ¼ 52.2, 106.2, and 142.2 . Based on the contours of vapour volume fraction in the 1st and 2nd rows, the cavity firstly occurs in the seat which is a 45 chamfer, then grows downstream and another cavity appears on the valve surface where there is an intersection between the spherical surface and the cylindrical surface.
In the 3rd row, there is a sharp drop in the static pressure across the gap between the set and the valve initially, then the lowpressure zone expands, finally the zone size decreases owing to increasing pressure in the pumping chamber.
In the 4th row, the mixture highest velocity rises, and the higher velocity region expands with increasing 4. The ratio of the highest velocity to the mean velocity is around 10, the flow through the gap is basically an annular jet.
In the 5th row, the mixture temperature in the cavity varies slightly with 4, resulting in 0.549 K temperature depression in maximum. This temperature depression is much smaller than 1.5e5.1 K depression in R114 cavitating flow in a venturi [49]. The liquid R245fa is subject to a weaker thermodynamic effect in cavitation under operating conditions considered in comparison with R114 in the venturi [49].
There has not been visual experimental data on the cavity in the valve of diaphragm pump. Therefore, the biggest cavity predicted at 4 ¼ 106.2 is qualitatively compared with the visualized cavity (white area) in a poppet valve in hydraulic oil systems [58] in Fig. 10. An isolated shedding cavity is found near the valve downstream the gap in the experiment. In the prediction, however, there is a sheet cavity over the valve surface but also an isolated shedding cavity off the valve and downstream the gap. Both the prediction and the experiment share a similar pattern in the cavity.

Discussion
Large scale ORC power plants have been relatively mature, and the current research is focused on small/micro scale ORC power plants (i.e., a few kWs) for waste heat recovery (e.g., energy recovery from IC engines) or distributed power generation, etc. For large scale ORC power plants, centrifuge pumps are employed to reduce risks of cavitation. Additionally, the cavitation can be supressed by a high inlet pressure. As a result, cavitation in centrifugal pumps is rarely an issue for large scale ORC power plants. However, for such small/micro scale ORC power plants, the feed pump is subject to a low mass flow rate. Diaphragm pumps are commonly used for small/micro scale ORC power plants due to their oil-free feature and good availability. However, it should be noted that e.g., diaphragm pumps are not optimised for refrigerant application, which is why this work aims to investigate potential areas of improvement, ultimately providing guidelines for developing suitable liquid pumps for small scale ORC power plants.
This work aims to have further understanding of the cavitation effect primarily at the pump inlet. The pump inlet pressure of an ORC system is a function of the condensing temperature and the system mass (liquid height etc). The expander load and thus high side pressure has negligible impact on the low-pressure side of the system, therefore, the significantly added complexity of an overall system approach would provide little further insight into the Table 4 The minimum liquid pressure and temperature, maximum temperature depression and vapour volume fraction on the seat and valve at various inlet pressures predicted by CFD simulations.

Inlet
Seat Valve  cavitation phenomenon. For this reason, there is no benefit to add the experimental/simulation results of the overall system in this paper. Currently all cavitation mitigation techniques (subcooling, liquid column height etc) come with ORC system design compromises, with even a small decrease in NPSHr leading to an overall system improvement.
In the article, the valve motion was involved in CFD simulations as a rigid body, described with Eq. (14). Fig. 11 presents a comparison between the CFD simulations and the 1D model predictions [23] at p 1 ¼ 141 kPa inlet pressure, including instant flow rate through the valve gap, valve lift and velocity and force acting on the valve by the fluid. In the two cases, the same diaphragm motion, such as Eq. (5), is applied, then the identical instant flow rate q is resulted as shown in Fig. 11a. But the rest parameters in Fig. 11bed calculated by the CFD simulations are far below the predictions by the 1st-and 2nd-order 1D models in [23]. The empirical flow coefficient used in the models [23] or CFX modelling method for flows around a valve might be responsible for the discrepancies. In CFX modelling of flows around a valve, a residual or dead gap must remain to make all the fluid domains connected. This gap reduces the pressure difference across the valve initially. In [27], Fluent was employed to simulate cavitating flows in a reciprocating pump suction chamber. A dead gap existed, but was closed with a wall, with increasing crank angle, the wall is switched to an internal boundary to allow the liquid to go through the gap. How to realize this switch function in the gap entrance needs furth investigations in ANSYS CFX.
The model constants proposed in the article vary with Reynolds  number at the pump inlet as expressed as Eq. (13). Their variation is plotted as a function of f in a CFD simulation in Fig. 12. Basically, F vap varies in the range of 4.653e4.711, and F con in the range of 0.00579e0.00588. In a reciprocating pump, there are two types of cavitation, i.e., the cavitation due to expansion of the pumping chamber at the beginning of suction stroke, and the cavitation at the maximum flow rate, which were visualized by using high-speed camera [59,60]. The former is called expansion cavitation and the latter is named as flow induced cavitation [59,60]. These two sorts of cavitation have been identified in the simulations, as shown in Fig. 8. The proposed modelling methods, and the adopted flow models have reflected the actual flow condition and the intrinsic feature of cavitation in the valve chamber.
For reciprocating pumps, their mean flow rate Q -NPSHa curves can be measured, and the NPSHa at 3 % drop in the mean flow rate in comparison with non-cavitation mean flow rate is considered as their NPSHr [12,61]. The mean flow rate Q -NPSHa curve based on CFD simulations is shown and compared with an expected ideal Q -NPSHa in Fig. 13. Here, NPSHa is based on the inlet pressure p 1 , the vapour pressure and liquid density at bulk temperature, p v ðT lb Þ, r l ðT lb Þ, and expressed by The mean flow rate rises rather than decreases after the inlet pressure p 1 85.5 kPa or NPSHa 0.70 m. After p 1 <85.2 kPa, the simulation in CFX is diverged with the ZGB cavitation model. As a result, the expected Q -NPSHa curve is unachievable. Hence, an updated cavitation model is desirable for the cavitating flows through the variable gap of the suction valve in a diaphragm pump.
The sudden change in flow direction and wall curvature can induce vortex in a flow field. The contours of fluid absolute helicity H are present in Fig. 14 where u ! is fluid velocity vector, and scalar H reflects the vortex generation or transportation in the flow field. The region with H>0 means that there is either vortex generation or transportation there. The figure illustrates that vortices are generated at the seat and propagated downstream along the valve surface. This phenomenon resembles to the observed vortex propagating downstream along the surface of a stationary poppet valve with sharp corner seat at a fixed valve opening [62]. Based on the colour in the contours, H value at p 1 ¼ 85.2 kPa is the maximum, then followed by the p 1 ¼ 89, 87, and 86 kPa cases, without any correlations. The volume-averaged absolute helicity H m and the volume-integrated H were calculated, but H m is nearly constant against NPSHa and subject to a very small coefficient of determination R 2 as shown in Fig. 15, thus H m does not show any correlation with NPSHa.
The specific entropy generation rate (EGR) in flow and temperature fields stands for the loss of exergy in the fields and indicates the devaluation of the energy in an irreversible process from a thermodynamics point of view. The specific EGR is defined by the following expressions [63e65]. (17) where is the total specific entropy generation rate, and are the specific entropy generation rates in the flow field by viscous dissipation and turbulence, respectively, is the specific entropy generation rate by temperature gradient in the temperature field, r is the mixture density, k is the turbulence kinetic energy, u is the specific rate of dissipation of turbulence kinetic energy, i and j are the coordinate direction index, i; j ¼ 1,2,3 for three-dimensional flow fields, x i and x j represent the Cartesian coordinates in the i and j directions, u i and u j are the velocity components in the i and j directions, and e ij is the time-averaged velocity strain rates. In Fig. 14, the largest specific EGR is found on the seat, where vortices are generated. The maximum specific EGR values rise from   The volume-integrated entropy generation rate in the valve chamber at crank angle 4 ¼ 106.2 was calculated at various inlet pressures. The volume-integrated entropy generation rate in the valve chamber, , is expressed by (18) where V is the total fluid volume in the valve chamber, note that , in which is the dissipation rate of energy defined in fluid mechanics. A -NPSHa plot is presented in Fig. 15.
varies little with NPSHa at NPSHa!0.81 m (p 1 !87 kPa). With cavitation development, i.e., NPSHa<0.81 m (p 1 <87 kPa), starts to rise quickly. This effect is similar to the -NPSHa of centrifugal pump predicted [45,66]. These scattered pointes can be correlated to NPSHa with the following formula: Since the coefficient of determination in Eq. (19) is significantly large, is correlated with NPSHa reasonably well. Fluid absolute helicity reflects the characteristic of transport of vortex lines in a fluid domain. It was assumed that this characteristic may be associated with cavitation intensity in the suction valve. The present results demonstrated that this assumption is improper, and the absolute helicity represents transport property of vortex during cavitation as discovered in [61].
The volume-integrated entropy generation rate is related to velocity shear strain rate , which is associated with viscous shear stress. Based on the new cavitation criterion that cavitation occurrence is decided by the principal stress in a liquid rather than the threshold of liquid static pressure alone [67,68]. Since the shear stress contributes to the principal stress in the liquid, undoubtedly, ffiffiffiffiffiffiffiffiffiffiffiffi 2e ij e ij q should have a link to cavitation situation in the suction valve.
In Section 4.1, we propose that a vapour volume fraction of 0.1 is employed as the criterion of cavitation inception in the suction valve. A comparison of predicted cavitation in the suction valve by using different criteria is made in Table 5. The cavitation inception in the suction valve is underpredict based on the criteria of 1 mm vapour volume streamwise length [54], and 1 ml vapour volume [55,56], but overpredicted by using a vapour volume fraction of 0.0001 [57]. Nonetheless, these criteria and their thresholds for the suction valve need to be validated by visual experiments in future.
In Fig.8(a)e(c), the a max on the seat is larger than that on the valve at p 1 ¼ 87, 86 kPa, but lower than that on the valve at p 1 ¼ 85.2 kPa. This effect is relevant to pressure profiles on the seat and valve at those pressures. The pressure contours on the seat and valve are illustrated in Fig. 16 at p 1 ¼ 87, 86, 85.2 kPa. The lowest pressure zone occurs only on the seat at p 1 ¼ 87, 86 kPa. However, an additional lowest pressure zone is seen on the valve at p 1 ¼ 85.2 kPa. This lowest pressure zone should be responsible to the larger a max appears on the valve at p 1 ¼ 85.2 kPa.
The article is subject to limitations. Firstly, even though the cavitation model constants F vap and F con were calibrated by using the cavitation flows of organic fluid R114 in a venturi, the cavitation behaviour of R114 might not be exactly identical to R245fa, additionally the cavitating flow in a stationary venturi can differ from the unsteady cavitating flow through the gap between the moving valve and the stationary seat. The use of the two model constants should be with caution. The results predicted here should be compared with the corresponding observations and measurements in future.
Secondly, organic fluids suffer from high concentrations of noncondensable gas, i.e., air [49], which can prompt or intensify cavitation nucleation. In the cavitation model constants calibration, the cavitation nucleus population density in R114 was regarded to the same as that in water. The actual non-condensable gas   (19) concentration in organic fluids needs to be measured and involved in cavitation models in future. Thirdly, the ZGB cavitation model is essentially in inertia-control cavitation regime. New cavitation models in the other cavitation regimes should be developed in future.
Fourthly, the idealized motion was postulated for the diaphragm and then used as moving boundary condition in CFD simulations. The influence of the diaphragm motion in CFD simulations need to be addressed in future.
At last, the valve opening, instant velocity and the force acting on the valve by the liquid predicted by CFD simulations differs considerably from those provided by 1D mechanical model [23]. To confirm these differences, experiments on valve motion are desirable in future.

Conclusion
In the paper, 3D, unsteady cavitating flows of organic fluid R245fa through the valve chamber of a diaphragm pump for ORC systems were investigated numerically in suction stroke with CFD simulations, which adopted the k-u turbulence model, the ZGB cavitation model with calibrated two model constants, the rigid body motion model for 1D motion of valve and moving mesh technique. The thermodynamic effect in cavitation of R245fa was considered approximately by using the two model constants calibrated with R114 cavitating flows in a venturi, liquid temperaturedependent vapour pressure, and solving energy equation. The cavitating flows of R245fa were simulated at a series of inlet pressures of 141, 89, 88, 88.6, 87, 86 and 85.2 kPa, and the cavitation inception and developed states were captured. The vortex production and entropy generation in cavitation were then discussed. The two calibrated cavitation model constants depend on Reynolds number and liquid bulk temperature. Cavitation inception happens on the chamfer of the valve seat based on 0.10 vapour volume fraction threshold, then on the valve surface along with significant oscillation in pressure and force on the valve, and valve opening, and velocity. Cavitation states are also related to crank rotational angle. Particularly, the expansion cavitation and flow induced cavitation occur at a low flow rate and a high flow rate in sequence at a certain inlet pressure. The temperature depression depends on inlet pressure and crank rotational angle with the maximum depression of 0.549 K. Vortices are generated at the seat and propagated downstream on the valve surface, but the volumeaveraged or volume-integrated fluid absolute helicity or itself does not correlate with cavitation states. The volume-integrated entropy generation rate correlates with cavitation states to some extent. The cavitation models in the other cavitation regimes, effects of non-condensable gas concentration on cavitation of organic fluid, and experiment work on the cavitation performance and valve motion of the diaphragm pump are on demand in future.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
R114 liquid and vapour densities r l , r v , specific heat capacities c pl , c pv , dynamic viscosities m l , m v , thermal conductivities l l , l v , and vapour pressure p v are approximated by the following expressions in terms of bulk liquid temperature T lb or local liquid temperature T l ranged in 240 K and 350 K Snice the numerical simulations of cavitating flows of R114 in the venturi is steady and stable in most cases, the specific heat capacities c pl , c pv and vapour pressure p v are all in terms of local liquid temperature T l .
The reference state of R114 is chosen to be at À6.7 C temperature and 66.82 kPa pressure, the vapour reference specific enthalpy and Molar mass are 143859.98 J kg À1 and 170.9 g mol À1 . The liquid reference specific enthalpy and entropy and the vapour reference specific entropy are set to be zero.

A2 CFD simulations
Because the experimental venturi is with circular cross-section, its geometry and the flow are circumferentially symmetrical, a sector with Df ¼ 5 central angle was used in CFD simulations. To guarantee a uniform boundary condition, the fluid domain is extended to 5D upstream and 10D downstream, respectively, as shown Fig. A2.
The flow models are composed of the Reynolds time-averaged Navier-Stokes equations, k À u two-equation turbulence model and the ZGB cavitation model as well as the energy equation. The boundary conditions include inlet boundary condition, no-slip smooth wall condition, symmetrical condition, and outlet boundary condition, see Fig. A2. The inlet boundary condition is subject to a given static pressure, bulk liquid temperature and zero gradients for the fluid velocities and turbulence variables. At the outlet, a known mass flow rate is given.
In [49], experimental conventional cavitation number K, approach or far field liquid velocity u 0 and bulk liquid temperature T lb are provided. The conventional cavitation number K for dynamic similarity of cavitating flows is based on the vapour pressure at bulk liquid temperature, and defined as [49].
Since T lb has been given, the liquid density r l ðT lb Þ and vapour pressure p v ðT lb Þ at T lb can be calculated with Eq. (A1) in a straight way, then the inlet liquid static pressure p 0 and mass flow rate m f are determined by the following expressions 8 > > < > > : where Df is the central angle of the fluid domain, D4 ¼ 5 , and p 0 is applied to the inlet but m f to the outlet. The experimental inlet condition in [49] and the counterpart in CFD simulations are listed in Table A1. In the table, there are 12 experimental cases that are simulated successfully in comparison with 14 cases in total [49]. The remaining two experimental cases are with 25.56 C bulk liquid temperature, a converged numerical solution cannot be obtained from them owing to significant thermodynamic effect at that high bulk temperature.
The mesh is wedge cell-dominated with a small number of tetrahedral cells. The inflation mesh is created near the walls to have a better resolution to the boundary layer on the walls. Three meshes were generated, i.e. mesh A (104985 cells, 8 layer inflation near wall), B (212184 cells, 22 layer inflation) and C (291889 cells, 32 layer inflation), which are with an averaged mesh quality of 0.378.
where p w is fluid pressure on the venturi wall. The predicted pressure coefficient agrees well with the measurement in the throat and diffuser inlet region, but slightly poor in the middle part of the diffuser wall. Since there is a jump in the curvature between the throat and the diffuser, a sharp drop in the coefficient is observed in the inlet to the diffuser. The pressure coefficient profiles predicted with mesh B and C are overlapped along the wall, and the mesh size independence is achieved in these meshes. Thus, mesh C is adopted in CFD simulations of cavitating flows of R114 in the venturi because these simulations need less computational effort.

A3 Correlations for the model constants
In [49], the observed cavity lengths L cav in cases 1e9 were 2 3 4 inch (69.9 mm), but 1 2 inch (12.7 mm), 1 1 4 inch (31.8 mm) and 4inch (101.6 mm), respectively. To secure a pair of suitable model constants F vap and F con for each experimental case in Table A1, the follow steps are adopted: 1) A simulation of isothermal two-phase flow of R114 and its vapour (zero volume fraction) with inlet boundary conditions and the mass flow rate applied to the outlet in a case is launched until a converged numerical solution is reached. 2) A pair of temporary F vap and F con are assigned, and the ZGB cavitation model and energy equation are turned on, then a simulation of cavitating flow is activated by using the converged solution in 1) as an initial condition. 3) If a converged solution is obtained, then the cavity length is measured in CFX-Post with the dimensional scaler. If not, a new pair of suitable F vap and F con are assigned, and go to 2). 4) If the predicted cavity length coincides with observed one within ±3 mm error, the suitable F vap and F con are redeemed, and the case is finished, or a new pair of F vap and F con are given, the procedures 2) and 3) are repeated, until the suitable F vap and F con are realised.
The determined suitable F vap and F con in 12 cases are tabulated in Table A2. The typical appearance of cavitated R114 predicted in cases 6 and 10 is compared with the observation is illustrated in Fig. A4. Since the visualization pictures of cavity fail to present cavity thickness information, just the cavity length was employed to judge the reasonability of the determined F vap and F con .
The model constants F vap and F con vary with bulk liquid temperature T lb and Reynolds number Re (¼u 0 Dr l =m l ), as shown in Fig. A5. Both F vap and F con decrease with increasing T lb and Re. Thus, the thermodynamic effect becomes more significant with increasing bulk liquid temperature and velocity.

A4 Temperature depression and minimum cavitation number
The predicted temperature depression DT l and minimum cavitation number K min are plotted as a function of bulk liquid temperature T lb and Reynolds number Re in Fig. A6. Even though the predicted temperature depression DT l rises with ascending bulk liquid temperature like the experimental DT l , they are in considerably different slopes. The model underestimates temperature depression at a lower bulk liquid temperature but overestimates it at a higher temperature.
The predicted minimum cavitation number K min is defined by the pressure in the cavity at local liquid temperature, and expressed as [49].
where Dp v ¼ p v ðT lb Þ À p v ðT l Þ is vapour pressure depression, K min is a parameter indicating cavitation similarity, i.e., if cavitation is similar in two flow systems at a Reynolds number for the same liquid, then K min is constant, or vice versa. As shown in Fig. A6, the experimental K min is nearly unchanged across tested bulk liquid temperatures and Reynolds numbers. The predicted K min , however, rises significantly with increasing bulk liquid temperature T lb and Reynolds number Re. The ZGB cavitation model in the inertiacontrol cavitation regime cannot predict a correct K min across a range of bulk liquid temperature or Reynolds number even though the cavity length is estimated properly.
In the calibration, the initial vapour bubble radius R 0 ¼ 10 À6 m and the density of cavitation nuclei number a 0 ¼ 5 Â 10 À4 for water are held, respectively. For R114, however, the initial vapour bubble size and density of cavitation nuclei number may be different from water. The influence of initial vapour bubble size and density of cavitation nuclei number on two model constants needs to be clarified in future.
The calibration of cavitation model constants for cavitating flows over a hydrofoil can be conducted based on an optimization algorithm [50,51]. However, the simulation of cavitating flow in a venturi often cashed. This property makes the automatic optimization of two model constants very difficult. Consequently, the new methods proposed in [50,51] were not employed. Obviously, this issue is worth being investigated in future.    A4. The comparison of cavitation appearance in experiment and CFD simulation at identical cavitation number and approach velocity, the left column is for case 6, the right column for case 10, the white and black pictures are cavitation visualization in the venturi [49], the colourful pictures are CFD simulation in the present article.  far field fluid pressure in a Venturi, Pa p 1 mean fluid pressure at the pump inlet, kPa p v saturated vapour pressure of liquid at temperature, kPa p w fluid pressure on the venturi wall, Pa q instantaneous flow rate of the diaphragm pump, l/min Q mean pump flow rate, l/min Re, the experimental data was credited to [49].
W. Li and Z. Yu Energy 237 (2021) 121495 R instantaneous radius of flat cone of diaphragm, mm R 1 R 2 two radii of diaphragm in flat cone shape, mm R b radius of vapour bubbles, m R b0 initial radius vapour bubbles or nucleation sites, m Re Reynolds number of the Venturi, Re ¼ u 0 Dr l = m l )

Re gap
Reynolds number in the gap between the seat and the valve body Re h Reynolds number at the pump inlet, Re h ¼ u h d h r l = m l specific entropy of liquid or vapour, J/(kg K) total specific entropy generation rate, W/(m 3 K) specific entropy generation rates in flow field by viscous dissipation, W/(m 3 K) specific entropy generation rates in flow field by turbulence, W/(m 3 K) specific entropy generation rate by temperature gradient in temperature field, W/(m 3