{'fraction_non_alphanumeric': 0.018333561362703515, 'fraction_numerical': 0.005655128380535413, 'mean_word_length': 4.738811829364041, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Computation of Flow in Screw Compressors

A CFD model enabling accurate and computationally affordable simulation of unsteady flow in screw compressors has been developed. This paper focuses on computational aspects, including real-gas CFD using hybrid structured/unstructured moving grids, and specifics of grid generation for moving rotors and their communication with the discharge plenum.


Introduction
A compressor is an essential part of the thermodynamic cycle in chiller systems. Rotary screw compressors are typically used in applications requiring more flow than is produced by reciprocating compressors but less than is produced by centrifugal compressors. They are produced in sizes that range from 10 to several 1000 cubic feet per minute. For chillers they are used for systems with the size ranging from 10 to about 500 tons of cooling capacity.
Screw compressors are positive displacement machines. The gas compression process is a continuous sweeping motion. Two or more helical screw rotors are intertwined in such a way that they create a closed volume between rotor lobes, casing and endwall. The rotation of the rotors decreases this volume continuously compressing the fluid. It is discharged through openings in the casing and endwall that are connected to a radial and axial port of the plenum, respectively. The ports open when a lobe passes one of these openings and releases flow into the discharge plenum. The flow discharge is very abrupt creating unsteady flow with large pressure pulsations. These create significant vibrations and noise issues.   To whom any correspondence should be addressed. Intensive research on gas pulsations in screw compressors have been carried out in both academia and industry. Fujiwara and Sakurai [1] measured gas pulsation, vibration, and noise in a screw compressor. Koai and Soedel [3] developed an acoustic model in which they analyzed flow pulsations in a twin screw compressor and related the pulsations to the compressor performance. Sangfors [4] and Huagen et al [2] developed mathematical models for the prediction of gas pulsations in screw compressor suction and discharge chambers. More recently, Kovačević and his coworkers [5], [6] developed a RANS-based three dimensional computational fluid dynamic model for the analysis and prediction of screw compressor performance.
In this paper, we present our three-dimensional CFD model for screw compressors used in chiller applications. We discuss the critical elements, which include rotor and plenum grid generation, nonmatching interfaces, the real gas implementation and non-reflective boundary conditions for specifying average outflow pressure. The structured rotor grid generation method builds upon earlier methods such as the one presented in Kovacevic [10], while providing an alternative approach for controlling grid quality in the region between rotors.
The CFD model is validated for a typical screw compressor. Comparison of the computations with experimental data shows that the presented method accurately predicts the pressure spectra. The computational model provides insight into the three dimensional flow field and can be used to evaluate compressor design modifications for pulsation reduction.

Computational framework
The in-house CFD solver used in these simulations is based on a second-order-accurate (both in space and time), explicit-in-time numerical method for compressible flow equations [12], with the komega RANS turbulence model. For refrigerants, the solver has been extended to model real gas properties via look-up tables. Moving grids are handled using a standard methodology [13]. Both, structured and unstructured grids can be used. The structured solver provides higher speed and accuracy, while the unstructured solver provides higher flexibility in grid generation. Communication between the two solvers is realized with a non-matching interface as the grids are generally nonmatching.
For real gas modeling, our in-house density-based compressible flow solver has been extended with NIST real gas property tables based on pressure and temperature. This approach avoids the construction of inverse tables for refrigerants. In order to analyze the screw compressor for a range of operating conditions, i.e. different discharge plenum pressures, a non-reflective outflow boundary condition that allows the specification of the mean pressure is required. The investigated computational model consists of the two screw rotors, radial and axial port and plenum with check valve and exit pipe, shown in Figure 2. In this particular case, the radial and axial ports are gridded with a structured mesh and the plenum and check valve are gridded with an unstructured mesh. Inflow conditions of total pressure and temperature are imposed at the axial inlet of the screw rotors. A time averaged static pressure is imposed at the outlet of the pipe.

Rotor grids
The surface of the screw compressor rotors is of complex geometrical shape. However, the basic surface of each rotor is defined by very few parameters: the profile of one lobe, the axis location, the number of lobes, the casing radius and the skewing angle along rotation axis. Figure 3 shows the surfaces and an axial cut through the rotors of a typical screw compressor.
The difficulty in the rotor grid generation arises from the fact that the two rotors are located very close to each other. In fact, they are touching each other at a few points. In the present approach, a very small clearance is always assumed between the two rotors. The rotors are rotating and an unsteady computation requires grids for each time step. The difference of the coordinate location of each node between successive time steps determines the grid velocity. If the nodes move relative to each other the grid deforms. To reduce moving grid errors which are proportional to grid velocity and cell deformation, the grid sequence has to be smooth in time. A structured meshing approach seems to be most appropriate as it gives the best control upon the grid node distribution.
Although the rotor surface is three dimensional, the cross section cut through both rotors normal to the rotation axis in Figure 3 reveals a two dimensionality of the geometry. Moving this cut along the axis produces an apparent rotation of the two surface lines. The rotation repeats itself -the geometry contains a periodicity in the axial direction. Therefore, the general grid generation strategy is to create a sequence of two-dimensional structured O-grids around each rotor and assemble them along the rotation axis into a three-dimensional grid. In the computation, the rotation of the rotor grids is achieved by manipulating grid indices and moving grid planes along the rotation axis.
There are two ways to connect the two O-grids around the individual rotors to each other (see Figure 4). A one-to-one connection is desirable as it preserves the fluxes at the interface. Such a connection introduces a relationship between the grid points on the surface of the two rotors constraining the grid generation algorithm. The circumferential length of the two rotors differs and the grid nodes on the rotors need to move to maintain the one-to-one connection between them. This leads to very deformed grids in the gaps, potentially eliminating the advantage of flux conservation at the interface that the one-to-one approach provides. A higher grid count circumferentially and normal to the rotor surface smoothes the grid movement but increases the computational cost.
On the other hand, the grid generation simplifies if the two rotor grids are connected with a nonmatching interface and the fluxes are interpolated. For this approach, the O-grids slide relative to each other at the interface and the deformation of the individual rotor grids is significantly reduced. In addition, a lower grid resolution is possible with this approach, although the high lobe curvature may lead to a local grid overlapping.
The grid generation process for the rotors is then as follows: a sequence of 2D Euler grids for the standalone rotors is created. Then an interface line is defined in the region where two rotors interact; this line is designed so that it changes smoothly in time. The standalone grids are then corrected in the interacting region by using interpolation along the wall normal gridlines. Finally, these 2D grids are assembled into 3D grids.  Figure 4. Grid between the rotors in the region highlighted in Figure 3 -for matching (left) and non-matching (right) interface.

2D grids for individual rotors
The 2D grid is generated from a grid for a single lobe defined by the area between the casing and the profile as shown in Figure 5, left. This area is gridded with 73x17 grid points where the points on the casing are spaced equidistantly and the points on the profile are spaced using a two sided exponential function. The first and last few points are aligned radially. Note that 73 points of circumferential resolution for one lobe is a good compromise between computational cost and geometry resolution, as can be seen in Figure 4. This one-lobe grid is then copied circumferentially. With 5 lobes for the male rotor there are 360 cells circumferentially. The grid rotates each time step 1 degree. A higher time resolution can be achieved by creating intermediate grids. In our computations we went up to 8 substeps for one degree of rotation.

2D grids for interacting rotors
Next step is to define a line that will separate the two interacting grids. It starts and ends at the stationary points where the rotor casings overlap ( Figure 6). This interface line should not cross the rotor surfaces and should change smoothly when the rotors rotate.
The interface line is found through the construction of a high resolution one-to-one connected auxiliary domain. An example is shown in Figure 6. Part of the domain is bounded by the rotor surfaces on which von Neuman type boundary conditions are applied. An elliptic solver is used that moves grid points on the rotor surface creating a smooth Euler type of grid.  Figure 6. Auxiliary domain for overlap region.
An example of an resulting grid is shown on the left in Figure 7. Note, that for some positions of the rotors a large mismatch of 5 to 1 cells is observed. This is the result of the elliptic solver placing the interface line in the middle of the domain. An improved interface line can be obtained by constraining it with the rotor casing radius. This is shown on the right of Figure 7 where now only 2 cells interface 1.
The final interacting rotor grids are interpolated from the standalone rotors and the interface line. Points are interpolated only along the lines that originate from the rotor surface. Practically, the individual rotor grids are 'squeezed' normal to the rotor surface ( Figure 8). First, the intersection of the grid line with the interface line is found. Then the length of the grid line from the intersection point toward the rotor surface is computed. Once this length is known points are distributed equidistanly or, if a viscous grid is need, using a two-sided exponential distribution for which a first cell height is specified on the rotor surface. A final 2D viscous grid with 73 x 33 gridpoints per lobe is shown in Figure 9.

3D rotor grids
In the final step, the 2D grids for interacting rotors are stacked in the axial direction to form a 3D grid. There are two ways of stacking the grids. When adding a 2D grid the circumferencial index can be kept constant or it can be rotated with the grids. The 'axial' gridlines are either parallel to the axis or follow the geometry features. Imagine moving the 2D grid plane axialy through the 3D grid: for the first method only the rotor surface rotates and the grid lines normal to the surface are 'squeezed' by the geometry features. For the second method, the grid rotates with the rotor and most of the grid only undergoes a solid body rotation. This is the reason why this stacking method has been chosen, as one can demonstrate that this reduces the numerical errors associated with moving grids. In addition, this method also allows easier enforcement of boundary conditions in the structured solver. The number of cells used to resolve one lobe determines the periodicity in the axial direction. When running the code in parallel, we split the grid in axial slices and provide each processor with one periodic slice. Each processor gets all the information necessary to reconstruct the full rotation of the rotor. If a smaller time step is needed intermediate grid can be generated or interpolated from the current grid. Splitting the grid in slices is not only efficient for the reconstruction of the grid but also for code efficiency in general. For a structured code there is usually one 'inner' index for which the data is most compact alligned in the cache of the computer. Optimal speed can be achieved if this dimension is the largest which for the rotor grid is the circumferential direction.

Plenum grid
Downstream of the rotors we model the radial and axial ports, plenum, check valve and exit pipe. Some parts of the plenum are gridded separately and MMI interpolation is used to connect the parts. Generally, the port grid resolution should be comparable to the rotor grid resolution. However, the rotor grids slide above the ports and refining the grid in the ports in the direction of rotation increases the temporal resolution of the port opening and as such the accuracy of the interpolation in a region where large fluid jets and flow gradients appear. The plenum is gridded with an unstructured Euler grid with a general cell size such that at least the first two frequencies of the pressure spectra are resolved. A finer grid is used in the narrow channels of the check valve.

Dynamic Mis-Matched Interface (MMI)
The flow solver exchanges data on a Mis-Matched Interface (MMI) by interpolating the flow variables. The additional complexity is that the boundaries of the two surfaces with mismatched grids that are interfacing with each other do not coincide and change in time. Thus the boundary of the interpolating region is difficult to specify. In our approach the interpolating region is defined by a flow-through criteria. A solid wall boundary condition is applied where the criteria is not fulfilled. This procedure is dynamic since at each time step the metric of the rotor grids is recalculated and the boundary of the interpolating region is re-determined.
Examples of MMI interfaces are shown in Figure 10. The axial port is shown on the left. The interface is flat but the grid on both patches is of very different density. The interface between the two rotors is shown on the right. Here the grid of both patches is of similar density but the interface is three dimensional with large gaps and overlaps. The challenge is to define an unique criteria that finds the correct flow through area and an accurate interpolation stencil for any of these interfaces. For moving interfaces, the physical location of the donor points as well as the interpolation stencil has to change in time smoothly. The detection criteria should not be binary; it should allow to distinguish between imperfect choices.
The user defines two surfaces which contain a flow through area. For every node from one surface, the MMI routine projects it onto a plane defined by the cell face of the other surface. After discarding the points for which the distance to the plane is above a certain threshold a detector is used to decides which cell face is closest. This is shown on the left of Figure 11 on the example of a triangular face. The projection P' of the point P creates with the nodes of the triangular face the areas S a , S b and S c . The sum S a +S b +S c is equal to the triangle area S only if the projection P' is inside the triangle. For points outside, the ratio (S a +S b +S c )/S departs from unity quickly (Figure 11, right). The projection point for which the ratio is closest to unity is the best choice. The area ratio S a /S also determines the weight with which the flow variable in point A is weighted for a linear interpolation.

Real gas formulation
Since the screw compressor is used for air conditioning applications the real gas properties of the refrigerant need to be considered. For this our in-house density based compressible flow solver has been extended with real gas property tables. The look-up tables are based on NIST data and all thermodynamic quantities are defined as functions of pressure and temperature. The essence of the implementation in our solver is that pressure and temperature are updated directly from the conservative variable changes and no inverse tables are used.
The governing flow equations are solved with conservative variables to provide a solution change . This solution change is then used to update the pressure and temperature by Which can be written as The derivatives and fluid properties here are obtained from the looked up table using the pressure and temperature from the current iteration. Once the pressure and temperature is updated, the fluxes in equation (1) are computed with the thermodynamic properties obtained with the updated pressure and temperature.

5
Outflow boundary condition Acoustic propagation is sensitive to boundary conditions. Non-reflective inflow and outflow conditions are essential for an accurate noise prediction. The non-reflective boundary condition treatment follows Hixon et al ( [7]). The correction to the flow variables is decomposed into a mean flow correction to maintain the desired mean flow quantities and a non-reflective correction for the transient waves (Giles [8]): The correction to the mean quantities at the boundary is constructed using the 1-D characteristics of Thompson ([9]) which is based on the mean flow difference from the target values: Validation The ability to predict noise in a screw compressor has been validated for a simplified model. The model consists of the main components: two screw rotors and the plenum with radial and axial ports, check valve and exhaust pipe. The working fluid is refrigerant R134a.
In the real machine the flow enters the screw rotors axially as well as radially through a suction inlet. This has been omitted and inflow conditions have been specified at the axial inlet to the screw rotors. We found that a crucial parameter for noise generation in a screw compressor is the start of discharge pressure, i.e. the pressure in the rotor lobes at the location of the ports just when the ports open. The pressure buildup in the rotors is significantly affected by the total conditions at the inlet. In the computation we adjusted these conditions such to match the start of discharge pressure. At the outflow, a time averaged pressure has been specified over a non-reflective boundary condition.
The rotor grids are viscous wall integration grids and adiabatic viscous wall boundary conditions are specified at the rotor casing. As described the rotor grids are constructed such that axial endwall of the screws rotors is not resolved. Therefore an Euler wall boundary conditions is specified on the screw rotor endwall. The ports also have wall integration viscous grids and the plenum is gridded with an Euler grid and slip wall conditions are applied.
A typical run for 10 screw rotor revolutions with computational grids with about 9.3 and 2 million structured grid points for the screw rotors and ports, respectively, and 0.120 million unstructured grid points for the plenum, requires about 2 days using 50 processors on a 2013 Linux cluster. This relatively low turnaround time enables CFD analyses of this type feasible within a product design cycle.
An instantaneous snapshot of the pressure field is shown in Figure 2. One can see the pressure buildup along the rotor lobes which are under an angle to the rotors axis. In the computation these pressure stripes move axially towards the plenum with an increase of pressure until they reach the radial port. Pressure waves bounce within the lobes as flow discharges.
Unsteady pressure history on the plenum wall has been stored for each time step at the locations labeled 1,2,3,4 in Figure 2. In a post-processing step these pressure histories have been decomposed in power spectra for the fundamental lobe passing frequency and its harmonics. The power spectra for the four locations are averaged and compared to corresponding experimental data in Figures 12-14. The dashed lines indicate a +/-5dB range around the test data for reference. Repeatability of the test data is +/-1dB. Figure 12 compares CFD results to test data for a nominal operating condition. It is seen that the fundamental and two subsequent harmonics are captured within measurement accuracy. Figures 13 and 14 show sensitivity to change in the back-pressure and rotor speed, respectively. It is seen that the CFD accurately captures sensitivity to these parameters. In particular, Figure 13 shows that the fundamental increases markedly with increased back pressure, and Figure 14 shows that higher harmonic content increases with speed.

Conclusions
A three-dimensional computational model of a screw compressor has been developed and validated for gas pulsation prediction.. It has been shown that the methodology accurately predicts the first three harmonics with computational cost suited for application in the product design cycle.
Several novel numerical items have been addressed and presented. These include the rotor and plenum grid generation, real gas property look up tables, non-reflective boundary conditions and mismatching interface interpolation. The hybrid grid approach with structured grids for the screw rotors and unstructured for the plenum has proven to be an efficient and accurate way to compute screw compressor flow field in an industrial environment.