Abstract
Rotating stall axial compressor is a difficult research field full of controversy. Over the recent decades, the unsteady tip leakage flows had been discovered and confirmed by several research groups independently. This paper summarizes the research experience on unsteady tip leakage flows and stability enhancement in axial flow compressors. The goal is to provide theoretical bases to design casing treatments and tip air injection for stall margin extension of axial compressor. The research efforts cover the tip flow structure at near stall that can explain why the tip leakage flows go unsteady and the computational and experimental evidences that demonstrate the axial momentum playing an important role in unsteady tip leakage flow. It was found that one of the necessary conditions for tip leakage flow to become unsteady is that a portion of the leakage flow impinges onto the pressure side of the neighboring blade near the leading edge. The impediment of the tip leakage flow against the main incoming flow can be measured by the axial momentum balance within the tip range. With the help of the theoretical progress, the applications are extended to various casing treatments and tip air recirculation.
1. Background and Motivation
Modern high-performance axial compressors in gas turbines, especially in aeroengines, are mostly unshrouded due to the high rotating speed of the shafts, which means that the tip clearance is necessary. The tip clearance brings additional complication into the corner of rotor blade and the casing. In 1996, Lakshminarayana [1] published a comprehensive illustration on the flow structures within rotor blade passages. The shock, the tip leakage flow (TLF), and the tip leakage vortex (TLV) are the three flow structures that are unique to the blade tip region. Two boundary layers, the end-wall boundary layer at casing and the blade surface boundary layer, provide a battle ground for all these flows interacting with each other. Please note that TLF is different from TLV in that in some cases only a portion of the TLF would be revolved into TLV and the rest of the TLF would either mix with the incoming main flow or leak again through the tip gap of the neighboring blade.
Lakshminarayana’s illustration is meant to present the situation at compressor’s design condition. It would become even more complicated when the compressor is throttled to close-to and near stall conditions, which will be the main theme of this paper. The tip clearance not only changes the flows structures at the design condition but also alters the Stall Inception and stalling behaviors. We noticed the debate on whether or not the stalling mechanism requires a nonzero tip clearance. Vo et al. [2] hypothesized that spikes (one of the two basic types of stall precursors) emerged when the interface between the TLF and the main flow spilled out of the leading edge of rotor blade. This hypothesis thus can only be true when the tip clearance is not zero. Pullan et al. [3] argued that spikes can be recognized as a vortex tube spanning from blade’s suction surface near the tip to the casing, a result of blade leading edge separation due to high incident angle of the incoming main flow, which can happen even with zero tip clearance.
As a matter of fact, most modern axial compressors are tip critical; that is, their rotating stall is initiated at the tip region. In other words, as long as the tip clearance is not zero, the tip leakage flow takes part in Stall Inception process. Thus, the hypothesis from either side of the debate is applicable to the compressors we are dealing with daily. Pullan’s argument, the flow separation at rotor blade’s suction surface near leading edge as the cause of the rotating stall, is a continuation of the traditional Emmons’ model [4] that is widely accepted in engineering practices.
More than ten years ago, we started to study the tip flow structures at near stall and estimate the location of the interface between the incoming main flow (MF) and TLF. At first, we were able to repeat the finding of others [4–7] that TLF would become oscillatory at operating points close to stall limit, while the whole compressor is completely stable. This brings an unexpected benefit, making the tip air injection become closed-loop without sensing spikes. We then realized that this had been done several years before us without knowing the existence of the oscillatory tip leakage flow [8–10]. In order to estimate the MF/TLF interface location, axial momentum equation within the tip region was then integrated over a sequence of control volumes, which ended up with a curve of cumulative axial momentum distribution along blade’s axial chord. This curve resembles a bell shape, whose peak point marks the “boundary” between main flow dominated region and the TLF dominated region. This makes it possible to compare the ability to extend the stall margin for various casing treatments, which is very helpful in screening out of bad designs from a large pool of casing treatment candidates in early design stage.
Note that the purpose of this paper is not to review the entire research field of compressor instability. Day [11] presented an excellent review recently as a Scholarly Lecture in the annual conference of International Gas Turbine Institute in 2017. In this paper, the tip leakage flow structure at close-to and near stall is reviewed in Section 2, followed by the features of oscillatory TLF in both low- and high-speed compressors (Section 3). The role of axial momentum and the “bell-shaped” curve are explained in Section 4. The applications of both the “bell-shaped” curve and the oscillatory TLF are reviewed in Section 5. Conclusions are given in the last section for convenience, together with discussions on future work.
2. The Tip Leakage Flow Structure in Close-to and Near Stall
It is necessary to clarify a few terms before reading further. As a compressor is throttled from large flow rate to stall, it would experience several steps. Figure 1 illustrates them along a typical compressor characteristic. Several points are identified in the figure. The point of peak efficiency (PE) is where the efficiency reaches the maximum, while the design point (DP) is the point that the compressor is designated to work. While these two often referred to the same point, there are occasions where these two are different. For instance, when a compressor has to work on a designated flow rate to match the other components in the same system (e.g., a gas turbine), its DP would be different from PE. For the cases studied in this paper, DP and PE are the same because all the cases are single rotors. This ensures that the flow angles around the blade are well organized, which forms a fair ground for comparison against the rest of the off-design points.
As a compressor is throttled to the stall limit along its characteristic, it goes through a segment on the curve with totally steady flow, followed by a segment with only the TLF becoming unsteady (close-to stall) before it reaches the point of near stall (NS). The NS point is a point right before stall happens. In experiments, a compressor at NS may ramp down into stall by itself without further closing the throttle. In CFD computation, NS is tricky to identify because it depends on how the state of NS is approached. When the compressor is throttled by gradually increasing back pressure step by step, NS is identified as the last point at which the time-accurate unsteady simulation can “converge” into a periodic solution, because no more solutions could be found when increasing the back pressure any further. This is a common practice in the compressor CFD community. The segment with only the tip leakage flow becoming unsteady is called “close-to stall” in this paper. In 1996, Graf [5] first observed periodicity and motion of the tip leakage vortex in his study based on CFD simulations of the flow through a single-stage compressor utilizing a 3D time-accurate Reynolds-averaged Navier-Stokes solver. In 2001, Mailach et al. [4] published test results of periodic fluctuation of tip leakage vortex in a large-scale low-speed compressor, which was named “Rotating Instability” thereafter. Bae et al. [6] studied the phenomenon in a cascade wind tunnel and explained the vortex oscillation with Crow Instability, a model that was used to explain the unstable vortex pairs behind the aircraft wings. Deng et al. [7] found that this vortex oscillation strongly depended on the size of the tip clearance, and one of its necessary conditions was that the trajectory of the vortex core impinged the pressure side of the neighboring blade.
The process through which a compressor ramped down into stall is called “Stall Inception.” There are two basic stall precursors in Stall Inception, modal waves and spikes. The distinguishing among “close-to stall,” “near stall,” and “Stall Inception” is very important, as stall precursors only happen in Stall Inception. Near stall is a point that belongs to the segment of close-to Stall, in which no spike would be formed yet. Therefore, according to the opinion of this paper’s authors, the oscillation of tip leakage vortex or the Rotating Instability is NOT part of Stall Inception, although it may trigger spikes or modal waves that later eventually initiate rotating stall.
The flow structure in the tip region is the key to understand the mechanism behind the flow phenomena. The early efforts included Adamczyk et al. (1993) [20] and Suder and Celestina (1996) [21]. At off design conditions, a large “blockage,” that is, a low-energy region within the blade passage, was found due to vortex breakdown caused by the shock. Their excellent work was so much ahead of their time that no much more knowledge in this regard had been gained since then until unsteady oscillatory TLF was found and confirmed.
The modern powerful tools of unsteady Navier-Stokes solvers empower us to revisit this flow structure at a level much closer to stall and from a viewpoint of unsteady flow. One of the rotors was Darmstadt Rotor 1, which is a transonic rotor of a single-stage transonic compressor rig at Technische Universität Darmstadt [22, 23]. The comparison of experimental and computational characteristics is depicted in Figure 1 [12]. Note that the steady simulation fits well with the test curve. While overpredicting the total pressure rise up to 2.5%, the unsteady simulation was able to make the computed stall limit very close to that of the test result (less than 1.5% of the test flow rate at stall). Figure 2 depicts the regions of influences by tip leakage flows. At PE, the loading is distributed over about 80% of the blade chord. The shock is attached to the leading edge and the tip leakage vortex flows out of the blade passage. At NS, the loading fluctuates within the first half of the blade chord. The shock is detached from the blade leading edge and interacts with the tip leakage vortex, causing a large “blockage” as seen in Adamczyk’s paper [20]. Note that the pictures at NS are only one instant of the oscillatory movement.
(a) Point PE
(b) Point NS
Particles were released within the tip gap to trace the trajectories of TLF in numerical results. Figure 3 compares them at PE and at one of the instants at NS. They are the 3D views of those in Figure 2. For the one at NS, a schematic is plotted in Figure 4 to summarize the flows at NS for clarity. It can be seen that the TLF are divided into two parts along the blade chord based on the location of shock on the suction surface. The first portion of TLF forms the vortex core and flows through the blade passage, while the second portion goes across the blade passage and hits the pressure side of the neighboring blade. It then splits there, a part of which leaks over the blade tip one more time.
(a) Point PE
(b) Point NS
3. The Unsteady Features of Oscillatory Tip Leakage Flow
The fact that the tip leakage vortex impinges the pressure side of the neighboring blade is believed to be one of the necessary conditions for TLF oscillation. Deng et al. [7] were the first to notice this in our research group when simulating a low-speed rotor. This had later been observed in serval other compressor rotors, including Darmstadt Rotor 1 as shown in Figures 2 and 3 and NASA Rotor 67 in Figure 5. The top row of Figure 5 lists six instantaneous contours of static pressure coefficients with the last one almost exactly repeating the first one. The bottom row is the corresponding pressure coefficient distribution in the pressure side of the blade. The low-pressure spot on the pressure side is marked as A1, which is clearly casted by the low pressure core of the tip leakage vortex. At time 0/30T, where “T” represents the time interval to store data set by the time-accurate CFD solver, A1 is located at the leading edge. Because of it, the pressure difference across the tip clearance is low, which weakens the tip leakage vortex and pushes the TLF/MF interface downstream. As time goes on (e.g., at 20/30T and 30/30T), the low-pressure spot A1 moves towards the trailing edge, while the high pressure regains its control on the leading edge causing the first half of the TLF/MF interface to swing back. Close to the end of one period, the low-pressure spot A1 moves out of the blade chord, the TLF/MF interface returns to its starting position, and a new A1 emerges at the leading edge of the pressure side. The entire process repeats. The frequency of this oscillation is calculated as 0.586 BPF (Blade Passing Frequency) for this particular rotor in rotor relative frame.
The features of the aforementioned unsteady process were observed in compressor experiments. The first feature was that the location of the highest amplitude should be at the leading edge of the pressure side, neither at the spot where the shock/vortex interact nor at the shock itself. It is not even at the starting place of the tip leakage vortex. Figure 6 depicts the comparison of the root mean square (RMS) of static pressure distribution at casing. The left was the result of phase-locked RMS obtained from experiments of Darmstadt Rotor 1, while the right was the numerical simulation [12]. Both match well with each other.
The second feature is signature frequency bands of the casing static pressure measurements. Figure 7 demonstrates the relation between rotating relative frame and the casing stationary frame using numerical results. The top row shows the time series and its frequency spectrum for a probe fixed on the rotor relative frame. The first peak at 3056.68 Hz or 0.57 BPF is the main frequency and the second peak is its harmonic. The second row is the data taken from a probe at the same axial location but fixed on the casing. The frequency components are much richer due to the modulation of frequencies between the relative frame and the rotor rotating frame. Figure 8 compares the frequency spectra between the experiments and CFD simulations [12]. Note that the test results contain noises from turbulent, the variations of the blade geometry and/or assembly, and so forth; thus the frequency spectra appear as a few signature bands, not as pronounceable individual peaks.
(a) Time series
(b) Frequency spectrum
The comparison of the oscillatory TLF introduced here and Rotating Instability as described in [8] is of interest. Since there is no chance to study the compressors that were reported to possess Rotating Instability such as those in [8, 24], we are not able to justify whether or not these two are the same. However, there are at least two features that are common to both phenomena: (1) both embark a signature frequency band; (2) the frequency band is about 0.5 BPF. Regardless of whether this frequency is measured in the rotor rotating frame or in the casing stationary frame, as long as it is about 0.5 BPF, it would appear as the tip leakage vortices (TLV) alternate their trajectories between neighboring blade passages. That is, when one TLV is within the blade passage, its neighboring one would impinge the blade’s pressure surface. Both features were observed in our low-speed compressor.
4. The Role of Axial Momentum and the “Bell-Shaped” Curve
According to Vo’s hypothesis [2], it is necessary to be able to estimate the location of the interface between the incoming main flow (MF) and the tip leakage flow (TLF). Since such an MF/TLF interface exists in the rotor’s rotating reference frame, how it looks like on the casing stationary reference frame becomes crucial to experimentists. Cameron et al. [14] demonstrated that on casing the complex 3D curvy oscillatory surface of MF/TLF in the rotor rotating frame can be observed as a straight line on the casing stationary frame, as seen in Figure 9 [14]. This is because, when observed at casing, the spatial and temporal variations within the rotor frame are all naturally averaged. Figure 9(a) depicts the experimental result of the casing streaklines, which was taken from a transparent window when the compressor was operating at Point B. In Figure 9(b), on the left, it is the computational result of one instant showing the streaklines on casing, while on the right the resultant streaklines after spatial and temporal averaged display flow patterns qualitatively the same as those of experiments. The rotor under investigation was the transonic rotor in University of Notre Dame in USA, abbreviated as ND-TAC.
(a) Experimental result
(b) Numerical result
Cameron et al. [14] then proposed a simple model to estimate the location of the MF/TLF interface on casing, as illustrated in Figure 10. The axial momentum equation near the tip clearance was simplified as the balance between the pressure due to main flow and the axial momentum due to tip leakage jet. The details of derivation, which can be found in [14, 15], are omitted here due to the space limitation of this paper. Despite its simplicity, the trend of Xzs versus the compressor’s incoming flow coefficient fitted well with that of the throttling process, which suggests that this simple model did capture some physics embedded in the stalling mechanism.
Since the MF/TLF interface is a 3D surface, such as the one in Figure 11, it is not surprising that using Xzs alone is not sufficient to correlate it to stall. A control volume approach is thus proposed to include the 3D effect into a new model [16, 19]. Unfortunately, it is impossible to establish an analytical equation like the one for Xzs due to the complexity of the 3D unsteady flow. The new model is indeed a method of postdata processing based on 3D unsteady Reynolds-Averaged Navier-Stokes solutions. It starts with the same strategy as the Xzs equation: observing the flows within the rotor while sitting on the casing stationary reference frame. The details can be described below.
Consider a series of discrete control volumes installed at the tip region between the casing and the rotor passage tip as illustrated in Figures 12 and 13. If we integrate all the linear momentum in axial direction, according to the Newton’s second law, such an integral would be equal to the total axial force acting on the control volume. Because this is done on a fixed and stationary control volume, it is equivalent to pitch-wise smear the spatial variation. For unsteady cases, the total axial linear momentum at every time instant should first be averaged before spatial integration. The end result is a number representing the net axial force on the control volume, as measured by an observer on casing stationary frame. If this number is positive, it means that the net axial force on this control volume is pushing the fluid downstream. If it is negative, it means that the net force is pushing the fluid upstream. Therefore, the control volume approach proposed here makes the justification of the flow stability become simple. The key is where the net force is zero. The closer this location is to the leading edge, the easier the compressor would run into stall.
The rotor in Figure 12 is NASA Rotor 67. The smooth casing is placed on the left, while the casing treated with 6 grooves is on the right. The bottom row depicts the entropy contours on the suction surface of the blade, indicating the radial depth of influence by the tip leakage flow. This is how the depth of the control volumes is decided. For instance, in the case of Figure 12, the control volumes are taken as deep as up to 90% span. A volumetric illustration of a typical control volume is given in Figure 13. Considering only the axial direction, the momentum equation in a finite control volume form can be written as
Here, refers to axial direction. Each control volume covers one pitch. The two periodic surfaces of each control volume are ignored in the equation due to the periodic boundary conditions for the single passage simulation.
By adding the right-hand side of the above equation together and averaging it over the oscillation period, we will obtain the net axial momentum (or force) on the control volume. Plotting all the net forces for all the control volumes, we obtain a distribution curve of local net axial momentum (or force) (Figure 14(a)). Two curves in the plot are for peak efficiency (PE) and near stall (NS) points, respectively. The location where the net momentum is zero for PE is at 38% axial chord, while the one for the NS is at 14%, much closer to leading edge as expected. In Figure 14(b), the curves of cumulative axial momentum for both the smooth casing (SC) and the 6-groove casing treatment (CG6 as in Figure 12) at the NS point of SC are given (the comparison between them will be explained later). The cumulative axial momentum is equal to the cumulative distribution of the net force acting on the enlarged control volume from the inlet of the first control volume gradually to the outlet of the last control volume on the curve. This cumulative momentum curve reaches its peak at which the local momentum curve crosses the zero. After that, the local momentum becomes negative and the net momentum starts to decrease. Therefore, the curve forms a bell shape. It is thus called “bell-shaped curve” or simply “bell curve” for breviate.
(a) Local axial momentum
(b) Cumulative axial momentum
The peak of the bell curve divides the tip region into two. The one in the front (near the leading edge) is dominated by the incoming main flow, while the other one on the back is dominated by the tip leakage flow. The location of the peak is more valuable than the peak value itself because the peak value involves the blade force that is so complicated that it is hard to explain the meaning of the peak value in general. If needed, it has to be carefully examined case by case.
To demonstrate the usefulness of the bell curves, seven points were chosen on the same characteristic of Rotor 67 and seven bell curves were generated, which corresponded to the throttling process of the rotor [16]. Figure 15 depicts the results. One can see that the locations of the peak move exactly as expected towards the leading edge monotonically as the rotor is throttled to near stall.
Figure 14(b) also depicts the comparison of the two bell curves for the smooth casing and the treated casing, respectively. Both are at the same flow coefficient, that is, the near stall point of smooth casing. It is clear that the 6-groove casing treatment is able to move the peak location further downstream compared to the smooth casing. Therefore, one can predict that the CG6 would extend the stall margin of the rotor by comparing the bell curves at the near stall point with smooth casing, without even being numerically calculated to CG6’s own stall limit. The feature could greatly reduce the need of unsteady CFD simulations when it comes to design or optimize casing treatments.
Both Figures 14 and 15 are CFD results. In the next section, we will provide examples that have experimental results to validate the bell curves and demonstrate how they can be used in design or optimization of casing treatments and tip air injections. However, it must be noted that the idea of bell curves relies on the hypothesis that rotating stall would be triggered by the MF/TLF interface spillage out of the rotor leading edge. There is at least one counter example. Houghton et al. [25] showed that when testing the effect of single groove’s axial location on the stall margin improvement (SMI), there are two peaks on the curve of SMI versus axial chord. The hypothesis of MF/TLF spillage cannot explain what happens in this case. However, other than this special case, the method of bell curves works well with many stability enhancement techniques.
5. The Examples of Application
In this section, we will demonstrate how to apply the bell curves to design or optimize stability enhancement techniques. Every example had been experimentally validated.
5.1. Casing Grooves
The experimental data of the first example is taken from Ross [26], in which a set of casing grooves were tested. There are three grooves available on the especially decided casing (called variable casing treatment, or VCT). One can open or block any of the grooves to make 7 combinations of casing grooves, which are marked as 001, 010, and so forth on Figure 16. Number 1 means that the groove is open and number 0 means that the groove is closed. The seven SMIs together with that of the smooth casing were listed in Figure 16. The total eight bell curves were calculated and plotted together in one chart. The SMIs can be grouped into four groups, A, B, C, and D, each of which corresponds to two cases of casing grooves that exhibit similar SMIs. The order of SMI groups is A>B>C>D. The eight bell curves can also be grouped into four groups, each of which contains two bell curves whose peaks located almost at the same axial location. Furthermore, the order of the four groups of peaks along the axial chord corresponded EXACTLY to the order of SMI groups. The bell curves were plotted based on the CFD results at the same incoming flow coefficient as the NS point of smooth casing. At this point, the flow fields in most of the cases with grooves were still steady, which saved lots of computation time and expenses.
5.2. Skewed Axial Slots
Four skewed axial-slot casing treatments with different geometries were chosen to check the applicability of the bell curve (Table 1), which were named as CT-a, CT-b, CT-c, and CT-d for short. The data were taken from a large-scale low-speed compressor test rig [18]. In this study, each individual slot component covers exactly the tip chord of the blade and has a skew angle of 60° along the same direction of blade rotation. The depth of all slots is 10 mm. The comparison of the bell-shaped curves for these four CTs is shown in Figure 17. The experimental results are listed in Table 1. The trend of those peaks of the bell curves matches very well with the trend of SMIs.
5.3. Tip Air Injection with Self-Recirculation
Up to this point, the experiments were done prior to CFD and the bell curves, so the results were used for validation purposes. Hereby we are giving an example where the bell curves were done before the experiments and predicted the test results of SMI. A type of casing treatment (SELF-INJ), which recirculates the air from the trailing edge into the leading edge, was designed as shown in Figure 18. A high-fidelity time-resolved CFD simulation was done to predict its SMI by gradually rising the back pressure. The stall limit was recognized as the last point before the CFD collapsed. At the same time, the same CFD result at the near stall point of the smooth casing was used to construct the bell curve and compare it with those of the smooth casing, a double-grooved casing (generated separately for another paper [16]), a five-grooved casing, and a SAS casing. The resultant bell curves are depicted in Figure 19 [19], which predicted that the SMI of SELF-INJ should be similar to the double-grooved CG13 (not shown in the figure), much less than the other two casing treatments. The test results are given in Table 2 [19], which confirmed the prediction in Figure 19.
6. Summary and Future Work
Research on oscillatory TLF and its applications are reviewed in this paper. The oscillatory TLF is a phenomenon that has been observed and studied by many researchers. As one of the research teams worldwide, we offered both CFD and test results and explained them from a view angle of flow structure at near stall. Since the oscillation of TLF was considered as a consequence of the axial momentum imbalance, a novel control-volume-based method, the bell curve, is proposed to estimate the axial momentum balance between the main flow and the tip leakage flow. The bell curve was then applied to make comparison of SMIs for various casing treatment methods and other stability enhancement methods that involve tip leakage flows.
These research efforts have been centered on one single purpose, that is, enlarging the stall margin of axial compressors. However, towards this goal, there is still a long way to go. Four ideas are proposed here for future research:(1)Momentum transport in radial direction from PE to NS(2)Identification of critical stage in multistage environment(3)Integrated optimization of casing treatment and blade design for high-pressure compressors(4)Active control of heavily loaded low-pressure compressors
The first idea is to extend the momentum analysis from tip region to all radial section of the blade span through the throttling process from PE to NS. It will help to understand the blade loading transition as a compressor approaches to stall. The second is to extend the current research from single rotors to multistages. These two ideas are basic research. The third and fourth ideas are for high-pressure compressor and low-pressure compressor, respectively. They are the applications of the first two ideas, with a hope to design safer modern compressors.
Nomenclature
| : | Density |
| : | Root mean square of static pressure |
| : | Mass flow coefficient of compressor |
| : | Shear stress on solid surfaces |
| : | Area |
| : | Tip axial chord |
| : | Tip clearance |
| : | Axial momentum |
| : | Force |
| : | Unit vector normal to given surface |
| : | Static pressure |
| : | Time interval to store data in unsteady CFD |
| : | Absolute velocity |
| : | Relative velocity |
| : | Axial location of zero shear region |
| : | Axial direction |
| : | Axial location of bell curve peak. |
| : Blade tip | |
| : | Casing surface |
| : | Leakage jet through the tip clearance |
| : | Positive axial direction |
| : | Negative axial direction. |
| BPF: | Blade passing frequency |
| CFD: | Computational fluid dynamics |
| CT: | Casing treatment |
| DP: | Design point |
| IGTI: | International Gas Turbine Institute |
| MF: | Main flow |
| NS: | Near stall point |
| PE: | Peak efficiency point |
| RMS: | Root mean square |
| SAS: | Skewed axial slots |
| SC: | Smooth casing |
| SFB: | Signature frequency band |
| SMI: | Stall margin improvement |
| TLF: | Tip leakage flow |
| TLV: | Tip leakage vortex |
| URANS: | Unsteady Reynolds-averaged Navier-Stokes |
| UTLF: | Unsteady tip leakage flow |
| VCT: | Variable casing treatment. |
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
The support of the National Natural Science Foundation of China through Grants no. 51676184 and no. 51606189 is acknowledged. The authors thank many colleagues and former students at Chinese Academy of Sciences who have contributed and are still contributing their efforts to enrich this ongoing research experience.