On the Mathematical Modelling of Adaptive Darrieus Wind Turbine ()

Palanisamy Mohan Kumar^{1,2*}, Kulkarni Rohan Ajit^{1}, Narasimalu Srikanth^{1}, Teik-Cheng Lim^{2}

^{1}Energy Research Institute @Nanyang Technological University, Singapore City, Singapore.

^{2}School of Science and Technology, Singapore University of Social Sciences, Singapore City, Singapore.

**DOI: **10.4236/jpee.2017.512015
PDF HTML XML
1,137
Downloads
3,142
Views
Citations

Darrieus wind turbines are experiencing a renewed interest for their application in decentralized power generation and urban installation. Much attention and research efforts have been dedicated in the past to develop as an efficient standalone Darrieus turbine. Despite these efforts, these vertical axis turbines are still low in efficiency compared to the horizontal axis counterparts. The current architecture of the turbine and their inherent characteristics limit their application in low wind speed areas as confirmed experimentally and computationally by past research. To enable and extend their operation for weak wind flows, a novel design of Adaptive Darrieus Wind Turbine (ADWT) is proposed. The hybrid Darrieus Savonius rotor with dynamically varying Savonius rotor diameter based on the wind speed enables the turbine to start, efficiently operate and stop the turbine at high winds. As the wake of Savonius rotor has a profound impact on the power performance of the combined rotor, the wake of two buckets Savonius rotor in open and closed configuration is reviewed. The current study aims to develop an analytical model to predict the power coefficient and the influence of other design parameters on the proposed design. The formulated analytical model is coded in python, and the results are obtained for the 10 kW rotor. Parametric analysis on the chord length and the diameter of the closed Savonius rotor is performed in search of an optimized diameter to maximize the annual energy output. Blade torque and the rotor torque are evaluated with respect to azimuthal angle and compared with conventional Darrieus rotor. The computed results show that peak power coefficient of ADWT is 13% lower than the conventional Darrieus rotor at the rated wind speed of 10 m/s.

Keywords

Wind Turbine, Low Wind, Analytical Model, Darrieus, Savonius, Adaptive

Share and Cite:

Kumar, P. , Ajit, K. , Srikanth, N. and Lim, T. (2017) On the Mathematical Modelling of Adaptive Darrieus Wind Turbine. *Journal of Power and Energy Engineering*, **5**, 133-158. doi: 10.4236/jpee.2017.512015.

1. Introduction

Renewable energy sources are increasingly popular for their emission-free power generation. The tremendous increase in the wind turbine installation is expected to continue in the future with the current worldwide installation of 468 GW [1] . The advancement in other renewable sources such as solar energy and inexpensive storage system led to the growth of decentralized power generation or distributed energy generation. Reduction in transmission losses and the immediate response to the local power demand are the notable advantages of distributed energy generation [2] . There is renewed interest in the development of Vertical Axis Wind Turbines (VAWTs) for their simple design, Omni directionality and east of maintenance [3] . Straight bladed Giromill or H-rotor turbines are preferred than other VAWTs, especially for their efficiency. Regardless of the above-said merits, the startup issues [4] and the lack of aerodynamic power regulation at high winds are immediate hurdles in their development. Multiple strategies, field tests, and computational studies spread all over the literature to address these two critical issues.

The past attempts contribute significantly for the enhancement of startup characteristics, yet one design has not been singled out as an implementable solution for both startup and over speed regulation without affecting the performance of the Darrieus turbine at higher Tip Speed Ratio (TSR). Rotors with cambered airfoils are found to generate higher starting torque at low wind speed compared to symmetric airfoils [5] . Albeit, the performance at high wind speed is curtailed due to increased drag in the downstream, the reduction in peak power coefficient is small compared with an overall increase in the annual energy output. New airfoils are designed to extend the fatigue life of the blade and to reduce the manufacturing cost of the blades [6] especially for urban turbines, where the flow is characterized by highly turbulent and multidirectional [7] . Conventional airfoils are modified by incorporating a cavity to enhance the lift at low Reynolds number (Re) [8] . The computational study on the trapped vortex airfoil shows significant improvement in starting torque, yet the performance deteriorates at higher TSR [9] . Recent proposals such as blade pitching, trailing edge flaps and morphing blades addresses both the starting issues and the over speed regulation, but the cost of precise mechanical actuators and the complex sensing elements limits its commercial application. Increasing the solidity by increasing the number of blades enhances the starting torque, but the performance at higher TSR will decrease due to the blade wake interaction from the preceding blade [10] . Blades with trailing edge flaps are proposed as a potential solution for low wind startup and to aerodynamically regulate the rpm of the rotor. The solution is practically complex and the improvement in stating characteristics is not attractive as they are not able to sustain the rotation. An elegant, low-cost solution is still in search of the question that is hovering around the development of Darrieus wind turbine for decades. The current study attempts to provide a solution by proposing a novel Darrieus rotor. The remaining part of the study is aimed at developing an analytical model for predicting the performance of the proposed rotor.

2. Working Principle of Adaptive Darrieus Wind Turbine

Of several solutions that are discussed above, hybrid Darrieus and Savonius turbine is a potential candidate that be redesigned to improve the low wind speed performance and over speed regulation. The conventional hybrid Darrieus turbine integrates a Savonius rotor to a common shaft with the Darrieus rotor. The strategy is that the high torque generated by the Savonius rotor accelerates the rotor to higher TSR. Similar to other concepts, the hybrid Darrieus-Savonius also suffers from poor performance when the Darrieus rotor accelerates beyond 1. The optimum TSR for a two-bladed Darrieus rotor lies between 3 to 5 [11] and the optimum TSR for Savonius rotor is 1. The Savonius rotor tends to generate resistive torque and in fact energy must be expended to rotate Savonius rotor for the TSR above 1. The mismatch between the optimum TSR for the two rotors severely degrades the performance at higher TSR. Practically, a conventional hybrid Darrius-Savonius rotor will not accelerate beyond 1.5 resulting in iota of improvement in annual energy output. Hence a novel design has been put forward to minimize the influence of Savonius rotor beyond TSR 1. The strategy is to transform the Savonius rotor into a shape that leaves minimum wake downstream without any resistive torque. Two bucket Savonius rotor can be transformed into a nominal cylinder if they are able to slide. The wake behind the downstream is axisymmetric with minimum width compared to other shapes. The wake width and the kinetic energy imbibed dictate the performance of the Darrieus rotor. Two-stage two bucket Savonius rotor offset at 90^{0} can improve the directional starting. The three operating configurations are shown in the Figures 1(a)-(c). At low wind speed Darrieus rotor torque
$\left({M}_{d}=+ve\right)$ and the Savonius rotor torque
$\left({M}_{s}=+ve\right)$ are in the same direction when the Savonius buckets a and b are arranged as shown in the Figure 1(a). As the TSR reaches 1 the buckets slide towards the axis of rotation to form a cylinder without any torque generation
$\left({M}_{s}=0\right)$ as shown in the Figure 1(b). In the extreme wind conditions and above the rated rotor rpm, the Savonius buckets slides in opposite direction generating resistive torque
$\left({M}_{s}=-ve\right)$ decelerating the rotor as shown in the Figure 1(c). The number of blades on the Darrieus rotor, orientation of the rotor to the oncoming wind, angular offset between the Savonius buckets and the Darrieus rotor, ratio between the diameters of Darrieus rotor to the Savonius rotor are the crucial design parameters that determine the starting characteristics. The angular offset between Darrieus blades and the Savonius buckets will have minimal impact on the performance, as two stages are arranged offset at 90˚. Hence from the structural perspective, the Savonius buckets can slide on the Darrieus blade connecting struts eliminating the requirement of additional structures. Thus, the ADWT has the capability to start the turbine at low wind speed, let it operate with minimal effect on the Darrieus

Figure 1. (a) ADWT rotor in open configuration; (b) closed configuration; (c) at high wind braking configuration; (d) flow over single stage Savonius; (e) flow over double stage Savonius; (f) flow sequence on Single stage Savonius.

rotor and decelerate the rotor when it rotates beyond the rated rpm. The construction and the mechanical arrangement are less complex making this concept commercially implementable.

3. Analytical Model of ADWT in Open Configuration (Open Savonius)

The analytical model of the ADWT in open configuration is similar to the conventional hybrid Darrieus-Savonius rotor. In order to simplify the development of analytical model, the Savonius buckets are arranged in line with the blades of Darrieus rotor. The simplified configuration for ADWT is with two Savonius buckets without an overlap mounted on the same axis with the Darrieus rotor with two blades. The velocity diagram of the ADWT is shown in the Figure 2.

Figure 2. Velocity diagram of ADWT rotor.

3.1. Wake of Savonius Rotor in Open Configuration (Conventional Two Bucket Savonius Rotor)

Before proceeding with the analytical solution it is indispensable to predict the wake flow pattern of a two-bladed Savonius rotor. Countless studies focused on the flow over the Savonius rotor rather than on the downstream wake. A detailed flow investigation and adequate knowledge is a must to generalize the wake structure on the downstream. It is possible to deduce the wake flow pattern on the downstream within the flight path of the Darrieus blades by examining the flow leaving the Savonius buckets. Typical flow of a two bucket Savonius rotor is presented in the Figures 1(d)-(f). The flow is characterized by a sequence in which flow is first attached to the convex side of advancing blade. As the rotation advances, the flow transfers from the advancing blade’s convex surface to returning blade’s concave side. The flow simultaneously enters in-between center shaft space, followed on the upstream flow on the convex side of the returning blade. The vortex shedding from the returning and advancing blade tips contributes to increased wake width. The wake pattern is highly turbulent due to alternating suction and pressure zone occurrence on the concave and convex side of the blade. The flow on a single stage Savonius rotor is different from the double stage rotor [12] . The flow tends to jump from the high pressure zone from one stage into the low pressure zone of the adjoin stage. A spanwise flow is induced due to the flow movement from one stage to another and widening the wake width when it leaves the rotor. The unsteady and highly turbulent wake when interacts with the Darrieus blades flight path would have dispersed to a large wake width from smaller diameter Savonius rotor. As the Savonius rotor diameter increases, the flow in fact has high energy and highly turbulent. Hence a pragmatic way is to assume the Savonius rotor as a bluff body or as an actuator disk absorbing a part of kinetic energy leaving behind low energy wake. The wake width of the Savonius rotor on the Darrieus flight path is assumed to be equal to the Savonius rotor diameter. In reality some portion of the flow may be accelerated due to vortices from the bucket tips. The Darrieus blades may not be able to extract energy from the vortices due to large Angle of Attack (AoA).

3.2. Mathematical Model

${V}_{\infty}=u{V}_{\infty i}$ (1)

And the equilibrium induced velocity is

${V}_{e}=\left(2u-1\right){V}_{\infty i}$ (2)

With ${V}_{e}$ as the input velocity for the downstream half-cycle of the rotor the induced velocity at the end of the streamtube is

${V}^{\prime}={u}^{\prime}\left(2u-1\right){V}_{\infty i}$ (3)

The relative velocity for the for the upstream half-cycle of the rotor, $-\text{\pi}/\text{2}<\theta <\text{\pi}/\text{2}$ , is given by the expression

${W}^{2}={V}^{2}\left[{\left(X-\mathrm{sin}\theta \right)}^{2}+{\mathrm{cos}}^{2}\theta {\mathrm{cos}}^{2}\delta \right]$ (4)

where $X=\omega r/V$ represents the local tip speed ratio. The general expression for the angle of attack is

$\alpha ={\mathrm{sin}}^{-1}\left[\frac{\mathrm{cos}\theta \mathrm{cos}\delta}{\sqrt{{\left(X-\mathrm{sin}\theta \right)}^{2}+{\mathrm{cos}}^{2}\theta {\mathrm{cos}}^{2}\delta}}\right]$ (5)

By equating the blade element theory and the momentum equation for each stream-tube

${f}_{up}{\left(\frac{V}{{V}_{\infty}}\right)}^{2}=\text{\pi}\eta \left(\frac{V}{{V}_{\infty}}\right)\left[\left(\frac{{V}_{\infty i}}{{V}_{\infty}}\right)-\left(\frac{V}{{V}_{\infty}}\right)\right]$ (6)

${f}_{up}u=\text{\pi}\eta \left(1-u\right)$ (7)

$\eta =\frac{r}{{R}_{D}}$ (8)

where ${f}_{up}$ is the function that characterizes the upwind conditions

${f}_{up}=\frac{Nc}{8\text{\pi}{R}_{D}}{\displaystyle {\int}_{-\text{\pi}/\text{2}}^{\text{\pi}/\text{2}}\left({C}_{N}\frac{\mathrm{cos}\theta}{\left|\mathrm{cos}\theta \right|}-{C}_{T}\frac{\mathrm{sin}\theta}{\left|\mathrm{cos}\theta \right|\mathrm{cos}\delta}\right){\left(\frac{W}{V}\right)}^{2}\text{d}\theta}$ (9)

${C}_{N}={C}_{L}\mathrm{cos}\alpha +{C}_{D}\mathrm{sin}\alpha $ (10)

${C}_{T}={C}_{L}\mathrm{sin}\alpha -{C}_{D}\mathrm{cos}\alpha $ (11)

Airfoil coefficients ${C}_{L}$ and $C$ are obtained from the wind tunnel test or from literature and interpolating for local Reynolds number and the local angle of attack.

Defining the blades local Reynolds number as $R{e}_{b}$ for local conditions given by

$R{e}_{b}=Wc/{V}_{\infty}$ (12)

The turbine Reynolds number will be

$R{e}_{b}=\left(R{e}_{t}n/X\right)\sqrt{{\left(X-\mathrm{sin}\theta \right)}^{2}+{\mathrm{cos}}^{2}\theta {\mathrm{cos}}^{2}\delta}$ (13)

For each blade in the upstream position, the non-dimensional force coefficients as functions of the azimuthal angle θ are given by

${F}_{N}\left(\theta \right)=\frac{cH}{S}{\displaystyle {\int}_{-1}^{1}{C}_{N}{\left(\frac{W}{{V}_{\infty}}\right)}^{2}\left(\frac{\eta}{\mathrm{cos}\delta}\right)\text{d}\zeta}$ (14)

${F}_{T}\left(\theta \right)=\frac{cH}{S}{\displaystyle {\int}_{-1}^{1}{C}_{T}{\left(\frac{W}{{V}_{\infty}}\right)}^{2}\left(\frac{\eta}{\mathrm{cos}\delta}\right)\text{d}\zeta}$ (15)

By integrating for the entire blade

${T}_{up}\left(\theta \right)=\frac{1}{2}{\rho}_{\infty}c{R}_{D}H{\displaystyle {\int}_{-1}^{1}{C}_{T}{W}^{2}\left(\frac{\eta}{\mathrm{cos}\delta}\right)\text{d}\zeta}$ (16)

The average half cycle of the rotor torque produced by N/2 of the N blades is given by:

$\stackrel{\xaf}{{T}_{up}}=\frac{N/2}{\text{\pi}}{\displaystyle {\int}_{-\text{\pi}/2}^{\text{\pi}/2}{T}_{up}\left(\theta \right)\text{d}\theta}$ (17)

The average torque coefficient will be:

$\stackrel{\xaf}{{C}_{Q1}}=\frac{NcH}{2\text{\pi}S}{\displaystyle {\int}_{-\text{\pi}/2}^{\text{\pi}/2}{\displaystyle {\int}_{-1}^{+1}{C}_{T}\left(\frac{\eta}{cos\delta}\right){\left(\frac{W}{{W}_{\infty}}\right)}^{2}\text{d}\theta \text{d}\zeta}}$ (18)

Thus, the power coefficient for the upstream half can be written as

${C}_{P1}=\frac{\omega {R}_{D}}{{V}_{\infty}}\stackrel{\xaf}{{C}_{Q1}}$ (19)

Similarly for the downstream half cycle. The relative velocity for the for the downstream half-cycle of the rotor, $\text{\pi}/\text{2}<\theta <\text{3\pi}/\text{2}$ , is given by the expression

${{W}^{\prime}}^{2}={{V}^{\prime}}^{2}{\left({X}^{\prime}-\mathrm{sin}\theta \right)}^{2}+{\mathrm{cos}}^{2}\theta {\mathrm{cos}}^{2}\delta $ (20)

where ${X}^{\prime}=\omega r/{V}^{\prime}$ represents the local tip speed ratio. The general expression for the angle of attack is

${\alpha}^{\prime}={\mathrm{sin}}^{-1}\left[\frac{\mathrm{cos}\theta \mathrm{cos}\delta}{\sqrt{{\left({X}^{\prime}-\mathrm{sin}\theta \right)}^{2}+{\mathrm{cos}}^{2}\theta {\mathrm{cos}}^{2}\delta}}\right]$ (21)

By equating the blade element theory and the momentum equation for each stream-tube

${f}_{dw}{\left(\frac{{V}^{\prime}}{{V}_{\infty}}\right)}^{2}=\text{\pi}\eta \left(\frac{{V}^{\prime}}{{V}_{\infty}}\right)\left[\left(\frac{{V}_{\infty i}}{{V}_{\infty}}\right)-\left(\frac{{V}^{\prime}}{{V}_{\infty}}\right)\right]$ (22)

${f}_{dw}{u}^{\prime}=\text{\pi}\eta \left(1-{u}^{\prime}\right)$ (23)

${f}_{dw}=\frac{Nc}{8\text{\pi}R}{\displaystyle {\int}_{\text{\pi}/2}^{\text{3\pi}/2}\left({{C}^{\prime}}_{N}\frac{\mathrm{cos}\theta}{\left|\mathrm{cos}\theta \right|}-{{C}^{\prime}}_{T}\frac{\mathrm{sin}\theta}{\left|\mathrm{cos}\theta \right|\mathrm{cos}\delta}\right){\left(\frac{{W}^{\prime}}{{V}^{\prime}}\right)}^{2}\text{d}\theta}$ (24)

${{F}^{\prime}}_{N}\left(\theta \right)=\frac{cH}{S}{\displaystyle {\int}_{-1}^{1}{{C}^{\prime}}_{N}{\left(\frac{{W}^{\prime}}{{V}_{\infty}}\right)}^{2}\left(\frac{\eta}{\mathrm{cos}\delta}\right)\text{d}\zeta}$ (25)

${{F}^{\prime}}_{T}\left(\theta \right)=\frac{cH}{S}{\displaystyle {\int}_{-1}^{1}{{C}^{\prime}}_{T}{\left(\frac{{W}^{\prime}}{{V}_{\infty}}\right)}^{2}\left(\frac{\eta}{\mathrm{cos}\delta}\right)\text{d}\zeta}$ (26)

${T}_{dw}\left(\theta \right)=\frac{1}{2}{\rho}_{\infty}cRH{\displaystyle {\int}_{-1}^{1}{{C}^{\prime}}_{T}{{W}^{\prime}}^{2}\left(\frac{\eta}{\mathrm{cos}\delta}\right)\text{d}\zeta}$ (27)

$\stackrel{\xaf}{{T}_{dw}}=\frac{\frac{N}{2}}{\text{\pi}}{\displaystyle {\int}_{\text{\pi}/2}^{\text{3\pi}/2}{T}_{dw}\left(\theta \right)\text{d}\theta}$ (28)

The downstream torque coefficient is given by

${C}_{Q2}=\frac{\stackrel{\xaf}{{T}_{dw}}}{1/2{\rho}_{\infty}{V}_{\infty}^{2}S{R}_{D}}$ (29)

${C}_{Q2}=\frac{NcH}{2\text{\pi}S}{\displaystyle {\int}_{\text{\pi}/2}^{\text{3\pi}/2}{\displaystyle {\int}_{-1}^{+1}{{C}^{\prime}}_{T}\left(\frac{\eta}{cos\delta}\right){\left(\frac{W}{{W}_{\infty}}\right)}^{2}\text{d}\theta \text{d}\zeta}}$ (30)

${C}_{P2}=\left(\frac{\omega R}{{V}_{\infty}}\right)\stackrel{\xaf}{{C}_{Q2}}$ (31)

${C}_{P-D}={C}_{P1}+{C}_{P2}$ (32)

For the Savonius rotor as shown in the Figure 2, the power coefficient can deduced as follows

Suppose pressure difference on retreating side is $\Delta {P}_{r}$

$\delta {Q}_{p}=\Delta {P}_{r}\mathrm{sin}\theta \cdot 2r\mathrm{cos}\theta \ast \text{Area}$ (33)

Assuming advancing side contributes to negative torque through drag

$\delta Q=\left(\Delta {P}_{r}\mathrm{sin}\theta \cdot 2r\mathrm{cos}\theta -\delta {Q}_{p}-\Delta {P}_{a}\mathrm{sin}\theta \cdot 2r\mathrm{cos}\theta \right)4rh$ (34)

$Q={\displaystyle {\int}_{0}^{\text{\pi}/2}\Delta \left({P}_{r}-\Delta {P}_{a}\right)\mathrm{sin}\theta \mathrm{cos}\theta \ast 4{r}^{2}h\text{d}\theta}$ (35)

$Q=2{r}_{s}^{2}{H}_{s}{\displaystyle {\int}_{0}^{\text{\pi}/2}\Delta \left({P}_{r}-\Delta {P}_{a}\right){\mathrm{sin}}^{2}\theta}$ (36)

Average power p is obtained by integrating torque from 0 to $\text{\pi}$

${P}_{s}=\omega \cdot Q=\frac{\omega}{\text{\pi}}{\displaystyle {\int}_{0}^{\text{\pi}}Q\text{d}\alpha}$ (37)

Normalized power coefficient can be given as

${C}_{p-s}=\frac{P}{0.5\rho {V}_{e}^{3}\left(4{r}_{s}^{2}{H}_{s}\right)}$ (38)

Let’s assume

$\frac{\Delta {P}_{r}}{0.5\rho {v}_{e-s}^{2}}={S}_{1}$ (39)

${V}_{e-s}$ is equivalent velocity or relative velocity, ${v}_{is}$ is absolute velocity.

${V}_{e-s}={v}_{es}-{v}_{is}$ (40)

Consider retreating side: $\Delta {P}_{p}$ value is known. Solving the integral,

${Q}_{p}=2{r}^{2}{H}_{s}\frac{1}{2}\rho {s}_{1}\left({v}_{es}^{2}\right)$ (41)

${Q}_{p}={r}^{2}{H}_{s}\rho {s}_{1}{\left({v}_{\infty}-{v}_{is}\right)}^{2}$ (42)

${Q}_{M}={r}^{2}{H}_{s}\rho {s}_{1}\left({v}_{\infty}^{2}+2{r}^{2}{w}^{2}-\frac{\text{\pi}}{2}{v}_{\infty}r\omega \mathrm{cos}\alpha -2{v}_{\infty}r\omega \mathrm{cos}\alpha -2{v}_{\infty}r\omega \mathrm{sin}\alpha \right)$ (43))

Similarly, for the advancing side:

$\Delta {P}_{a}/0.5\rho {v}_{\infty}^{2}={S}_{2}$ (44)

Again, resolving the component will yield the equation of the torque for advancing side.

${Q}_{D}=2\rho {r}^{2}{\alpha}_{s}^{2}{H}_{s}{S}_{2}$ (45)

Combining all the parts together the final equation for the ${C}_{p}$ is given as

${C}_{P-S}=\left[\frac{{S}_{1}\phi}{4}-\frac{{S}_{1}{\phi}^{2}}{2\text{\pi}}+\frac{{S}_{1}-{S}_{2}}{{S}_{1}}\left(\frac{{\phi}^{3}}{4}\right)\right]$ (46)

The power coefficient of ADWT rotor in open configuration

${C}_{P}={C}_{P-D}+{C}_{P-S}$ (47)

4. Analytical Model of ADWT in Closed Configuration (Cylinder)

The power coefficient under the influence of Savonius rotor in closed condition can be derived by treating it as a nominal cylinder placed in a steady and homogenous wind flow. Even this simplified assumption gives rise to complex flow wake structures. The objective of the model is to predict the wake width and the velocity deficit due to the presence of cylinder. The stream tubes that are influenced by the wake width are identified and the input velocity is modified with the velocity deficit calculated during the iteration. However there will be axisymmetric flow acceleration on either side of the wake which is not accounted. The cylinder can be considered as non-rotating as the cylinder TSR is low, and the flow field displayed by both rotating and non-rotating flow for low rpm of 80 ~ 100 is similar. The current approaches and the assumptions will lead to the development of an analytical model that can be well integrated into the existing subroutine coded in python.

4.1. Wake of Savonius Rotor in Closed Configuration (Cylinder)

A wake boundary occurs between two fluids carrying different momentum along with their flow. The momentum deficit in the particular region of unidirectional flow is highly unstable and gives rise to the zone of turbulence mixing layer downstream at a point where the two streams meet for the first time. Though at far downstream the static pressure tends to equalize within the flow and wake, the velocity or the momentum deficit continues to travel along with the flow.

Detailed studies by the previous researchers have emitted a very important conclusion which relates the cylinder wake, its growth to the ratio of rotational to rectilinear speed ratio. The experimental studies by Prandtl [13] , Dfaz [14] gave the most vital results which suggested that the eddies behind the cylinder in terms of Karman vortex street disappears at higher rotational to rectilinear speed ratio. One of the most dominating reasons for this is, for the low values of rotational speed ratio (TSR), the vortices are shed in the alternate fashion behind the cylinder. The cylinder Re dictates the wake pattern behind the cylinder and the wake structure can be established by defining the Re. Alternating eddies are formed and moves along the direction of rotation progressively decrease until the TSR reaches 1. After that it starts dissipating completely transferring the turbulent kinetic energy to large eddy structures leaving behind constantly growing wake. In the case of VAWTs in combination with the central shaft the similar kind of pattern is observed. The flow structure behind the cylindrical shaft rotating at TSR greater than 1, creates a uniformly growing wake without breaking into vortex structures. So for the mathematical modeling, it is safe to assume that the wake structure behind the central shaft is in accordance with the two-dimensional wake behind a single body as suggested by H. Schlichting [15]

4.2. Mathematical Model

$\frac{\partial U}{\partial X}+V\frac{\partial U}{\partial Y}=-\frac{1}{\rho}\frac{\partial P}{\partial X}+\gamma \frac{{\partial}^{2}U}{\partial {Y}^{2}}$ (48)

The above equation is the governing partial differential equation

$-\frac{\partial P}{\partial y}=0$ (49)

$\frac{\partial U}{\partial X}+\frac{\partial V}{\partial Y}=0$ (50)

For cylinder wake which is relatively thick we cannot consider

$\frac{\partial P}{\partial X}=0$ (51)

Assuming, X is considerably large

$\frac{U}{{V}_{\alpha}}=1$ (52)

$\frac{V}{{V}_{\alpha}}=0$ (53)

${U}_{d}$ the difference between the velocities $\left({V}_{\alpha}-U\right)$ should be an even function (symmetric wake)

$V$ must be an odd function (asymmetric) in order to get the solution for governing differential equation.

Drag prediction from Boundary layer assumptions

$D={\rho}_{\propto}U{\displaystyle \underset{-\propto}{\overset{+\propto}{\int}}\left({V}_{\propto}-U\right)\text{d}y}$ (54)

Ud2 can be neglected from

$D={\rho}_{\infty}V{\displaystyle \underset{-\infty}{\overset{\infty}{\int}}\left({V}_{\infty}-U\right)\text{d}y}$ (55)

$D={\rho}_{\infty}V{\displaystyle \underset{-\infty}{\overset{\infty}{\int}}{U}_{d}\text{d}y}$ (56)

Comparing it with conventional drag equation

$\frac{1}{2}{\rho}_{\infty}{V}_{\infty}\cdot d\cdot {C}_{d}={\rho}_{\infty}{V}_{\infty}{\displaystyle \underset{-\infty}{\overset{\infty}{\int}}{U}_{d}\text{d}y}$ (57)

$\underset{-\infty}{\overset{\infty}{\int}}{U}_{d}\text{d}y}=\frac{1}{2}{V}_{\infty}\cdot d\cdot {C}_{d$ (58)

$\frac{{U}_{d}}{{V}_{\infty}}\propto \frac{{C}_{D}d}{2b}$ (59)

From Prandtl’s mixing length theory we know that

${v}^{\prime}={u}^{\prime}\cdot \text{const}=\text{const}\cdot l\cdot \text{d}u/\text{d}y$ (60)

${v}^{\prime}+{u}^{\prime}$ are turbulent velocity components. Also rate of increase of width ‘b’ of mixing zone is proportional to transverse velocity ${v}^{\prime}$

$\frac{Db}{Dt}=\text{const}\cdot \frac{l}{b}{u}_{d}=\text{const}\cdot {U}_{\mathrm{max}}$ (61)

$\frac{Db}{Dt}\propto {v}^{\prime}=\frac{Db}{Dt}\propto l\frac{\partial u}{\partial Y}$ (62)

Average value of $\frac{\partial u}{\partial Y}$ i.e.; velocity deficit along Y direction is considered to be proportional to transverse velocity $\frac{{U}_{\mathrm{max}}}{b}$

$\frac{Db}{Dt}=\text{const}\cdot \frac{l}{b}{u}_{\mathrm{max}}$ (63)

Let

$\frac{l}{b}=\Delta =\frac{Db}{Dt}=\text{const}\cdot \Delta \cdot {U}_{\mathrm{max}}$ (64)

For two dimensional

$\frac{Db}{Dt}=\text{const}\cdot \frac{l}{b}{u}_{d}=\text{const}\cdot \Delta \cdot {U}_{d}$ (65)

Equating above expressions

$\frac{Db}{Dt}=\text{const}\frac{\text{d}b}{\text{d}x}$ (66)

$\frac{\text{d}b}{\text{d}x}\propto \Delta \frac{u}{{V}_{\infty}}$ (67)

Introducing equation of wake width variation into drag equation

We derived previously

$2b\frac{\text{d}b}{\text{d}x}\propto \Delta {C}_{D}d$ (68)

Using variable differential form

$b\propto {\left(\Delta x{C}_{D}d\right)}^{1/2}$ (69)

We know from 2-D incompressible flow equations

$\frac{\partial U}{\partial t}+U\frac{\partial U}{\partial x}+V\frac{\partial U}{\partial Y}=\frac{1}{\rho}\frac{\partial C}{\partial Y}$ (70)

From Prandtl mixing length theory

$-{V}_{\infty}\frac{\partial {V}_{d}}{\partial x}=2{l}^{2}\frac{\partial {U}_{d}}{\partial Y}\frac{{\partial}^{2}{U}_{d}}{\partial {Y}^{2}}$ (71)

We introduce $K=\frac{Y}{b}$ as independent variable. (m=constant)

$b=m\cdot \sqrt{\Delta}{\left({C}_{D}dx\right)}^{1/2}$ (72)

${U}_{d}={V}_{\infty}{\left(\frac{x}{{C}_{D}d}\right)}^{\frac{1}{2}}f\left(k\right)$ (73)

Put this function again in differential equation to get value of $f(\; x\; )$

$f\left(k\right)=\frac{m\sqrt{\Delta}}{18{\Delta}^{2}}{\left(1-{k}^{3/2}\right)}^{2}$ (74)

And from momentum equation of experimental data the constant $m\sqrt{\Delta}$ yields the value equivalent to $10\sqrt{\Delta}$

$b=\sqrt{10}\Delta {\left(x{C}_{D}D\right)}^{1/2}$ (75)

$\frac{{U}_{d}}{{V}_{\infty}}=\frac{\sqrt{10}}{18\Delta}{\left(\frac{x}{{C}_{D}d}\right)}^{1/2}{\left(1-{\left(\frac{y}{b}\right)}^{3/2}\right)}^{2}$ (76)

According to experiments by H. Reichardt

$\frac{{b}^{\prime}}{2}=\frac{1}{4}{\left(x{C}_{D}d\right)}^{1/2}$ (77)

Substituting $y=0$ in the ${U}_{d}$ equation then the velocity deficit is at central region which practically explains the maximum deficit velocity since the distribution is close to Gaussian distribution.

From Boundary layer theory and momentum equation after neglecting small terms we get

$\frac{{C}_{DT}}{4}=\frac{\Delta {V}_{m}}{{V}_{e}}\cdot \frac{L}{{D}_{T}}$ (78)

${\left({V}_{e}/\Delta {V}_{m}\right)}^{2}=\left(x+{x}_{0}\right)/\left(b{D}_{T}\right)$ (79)

For small Reynolds number of 10^{4}, the following result holds true

${\left(L/{D}_{T}\right)}^{2}=c\left(x+{x}_{1}\right)/{D}_{T}$ (80)

$\frac{\Delta {V}_{e}}{{V}_{e}}=\frac{\Delta {V}_{m}\cdot 2L}{{V}_{e}{N}_{T}{R}_{T}f{\theta}_{0}}$ (81)

5. Results and Discussion

The developed mathematical model is applied to 10 kW ADWT and conventional straight bladed Darrieus turbine. The dimensional details of the configurations are listed in Table 1, and the various configurations are shown schematically in Figure 3. The results are evaluated for the starting characteristics, blade and rotor torque and power coefficient variation. Though the starting characteristics does not reflect the actual starting conditions, it will provide insight into the low Re behavior of cylinder. Another possible configuration that may be of interest for the current study is the height of the Savonius rotor with respect to the height of the Darrieus. It is envisaged that half-length of the Savonius rotor may have minimal influence on the performance of the Darrieus rotor. On the negative aspect, the half-length Savonius may induce uneven loading on the

Table 1. Dimensional details of the investigated rotor.

Darrieus blade on the downstream half when it enters the wake zone of the cylinder, yet the advantage that can be reaped in the high TSR is attractive to incite an investigation. Blade torque and the rotor torque values are evaluated for the above-said architectures and compared with the conventional Darrieus turbine.

5.1. Parametric Study

A parametric analysis is conducted to evaluate the effect of solidity on the power and torque coefficient on a conventional Darrieus rotor. The solidity is varied by changing the chord length of the blades. The chosen chord for the study is 20 mm, 40 mm and 70 mm. The results are plotted as a function of TSR as shown in Figure 4. It is apparent from the result that, the peak power coefficient increases with a decrease in solidity but the operating range of TSR decreases. A maximum Cp of 0.46 is achieved at 3.8 TSR without accounting the effect of struts. For 70 mm chord length, the Cp decreases by 12.8% to 0.4. Though a low solidity rotor is preferred from the aerodynamic perspective, a thicker blade with increased chord length is favorable from structural perspective and desirable for starting characteristics. A longer blade chord with high thickness will enable the blade to withstand alternating fatigue stress and thus extending the operating life time with acceptable loss in peak power coefficient. Starting characteristics are compared for the conventional Darrieus rotor and ADWT rotor at low a wind speed of 2 - 4 m/s. The outcome of the prediction is plotted as Cp vs TSR as shown in the Figure 5(a). The conventional Darrieus rotor displays earlier starting than ADWT in closed condition. Low wind speed startup study is to shed additional light on the performance of conventional Darrieus rotor, since

Figure 3. (a.1) Wake pattern behind a Savonius rotor at 90˚ and 0˚; (a.2) wake behind a Darrieus turbine at high TSR; (b) velocity deficit of ADWT rotor in closed condition; (c.1) full length Savonius configurations; (c.2) half-length Savonius configuration; (c.3) ADWT in closed configuration without end plates.

Figure 4. (a) Rotor torque for half length Savonius; (b) Cp vs TSR for various chord lengths; (c) Cp vs TSR at low wind speed.

Figure 5. (a) Blade forces vs Azimuthal angle for half-length Savonius; (b) blade forces vs Azimuthal angle for full-length Savonius; (c) rotor torque for half-length Savonius.

the starting torque will be generated by Savonius rotor if the ADWT is in open condition.

5.2. Blade Torque and Rotor Torque

The normal and the tangential forces are plotted as a function of azimuthal angles shown in the Figure 5. The well-established pattern of double peak is displayed by both conventional and ADWT rotor. The non-dimensional force coefficients ${F}_{T}$ and ${F}_{N}$ are computed from the predicted normal and tangential force. The normal force coefficients are lesser in the downstream half than the upstream half. The tangential force and normal forces are predicted for both full and half-length Savonius rotor. From the figure it is evident that the presence of cylinder increases the AoA resulting in the decrease of lift and increase of drag. The maximum tangential force is obtained at 0˚ as 900 N. The tangential force difference between conventional Darrieus and the ADWT is lesser than the difference in normal force. The dynamic stall has significant influence on the normal force.

5.3. Power Coefficient and Torque Coefficient

The power coefficients are evaluated for various diameters and the conventional Darrieus rotor. The diameters are investigated for the full length and half-length Savonius rotor including the blade tip loss as shown in Figure 6(a) and Figure 6(b) respectively. The power coefficient curve follows the same trend for all the diameters investigated and the difference in their magnitudes is diminutive. A maximum reduction in the power coefficient of 5% is reported for the cylinder of 2000 mm compared to the conventional Darrieus rotor. The peak power coefficient achieved by the conventional Darrieus rotor is 0.42 at TSR 2.5. Also, the difference between the half-length and the full-length Savonius rotor is almost negligible. Hence full-length Savonius rotor is preferred aerodynamically from the starting perspective. The power coefficient is anticipated to reduce drastically as explained in past literature, but in reality, the reduction is negligible. The potential reason is a strong correlation between the wake width, AoA with respect to azimuthal position. The same trend is followed in the torque coefficient prediction for half and full length Savonius rotor. A maximum ${C}_{t}$ of 0.17 is achieved at 2.2 TSR for conventional Darrieus whereas for the configuration with cylinder the maximum ${C}_{t}$ is less than 0.16. Hence it can be concluded from the power coefficient and torque coefficient prediction that the diametrical ratio between the Darrieus rotor and the Savonius rotor in closed configuration can be as high as 1:0.5 with acceptable loss in power coefficient.

6. Experimental Verification of Analytical Results

To validate the derived analytical model, the computed results are compared against the experimental results. An analytical model of the closed condition can be predicted close to the reality, experimentation has been carried out with dif-

Figure 6. (a) Power coefficient comparison with half-length Savonius; (b) power coefficient comparison with full-length Savonius; (c) torque coefficient comparison.

Figure 7. Closed Savonius diameters (cylinders) under investigation.

ferent diameter cylinders that will represent the closed Savonius rotor integrated with Darrieus rotor. The performance of the Darrieus rotor with 80 mm, 115 mm, 130 mm, 150 mm is compared with the conventional Darrieus rotor. The rotors are investigated for the Reynolds number of 268,737 corresponding to the wind speed of 9 m/s, as unsteady and flow separation behavior at low Re is more pronounced at low speeds. The Darrieus rotor consists of two blades of NACA 0018. The cylinders under investigation are made in two halves, so that the cylinders can be interchanged without disturbing the Darrieus rotor arrangement. The Darrieus rotor diameter is 400 mm and height is 300 mm. The various diameters employed in the wind tunnel test along with the Darrieus rotor are shown in the Figure 7.

Subsonic open circuit wind tunnel used for this study has a square cross-section of 700 mm × 700 mm. The wind tunnel is powered by 900 mm diameter axial flow fan (Multi-Wing) with the power capacity of 6 kW. A variable frequency drive regulates the wind speed from 0 - 15 m/s. The settling chamber with flow straightener streamlines the air flow. The flow velocity in the test section is uniform within 0.1%. The airfoils can be interchanged in the test setup without disturbing the center shaft or end plates to minimize the errors induced by shaft misalignment. The test turbine is mounted on the aluminum shaft with deep groove ball bearings on bottom end and spherical bearings on top end as shown in Figure 8 to accommodate for shaft misalignment. The rpm of the rotor was measured by proximity sensor. Hotwire anemometer (KANOMAX) with the accuracy of 0.1 m/s was mounted at a distance of 300 mm from the bell mouth exit downstream. The braking torque was applied by magnetic particle brake (Placid Industries). The required tip speed ratio was achieved by varying the

Figure 8. 150 mm diameter cylinder mounted on Darrieus rotor.

amperage to the magnetic particle brake. The difference between the rotor torque and the braking torque was measured by the torque sensor (LORENZ Transducers). The output data from the above sensors are logged by DEWE 43V data acquisition system. The test rotor has a frontal swept area of 400 mm × 300 mm and occupies nearly 39% of the test section area, which is more than the allowable blockage. The total blockage including solid and wake blockage for the non -standard shapes are accounted. The experimental values are blockage corrected before comparison with the analytical outcome.

The analytical computation and the experimental results are compared as shown in Figure 9(a) & Figure 9(b). For conventional Darrieus (bare turbine), the analytical model overpredicts by 68%, whereas for the 80 mm diameter cylinder, the difference in prediction is 53.8%. As the cylinder diameter increases the difference in prediction decreases, as noted from 115 mm, the difference diminishes to 36.6%. Conclusively, the analytical model overpredicts the power coefficient. The experimental values may tend to be lower due to unavoidable errors in the experiment. The operating Re is still in the lower range and the aerodynamic coefficients obtained thru Xfoil is relatively higher than the actual value.

7. Conclusion

An analytical model has been developed for the proposed novel Adaptive Darrieus Wind Turbine and the results are obtained for 10 kW rotor. The mathematical model has been developed for ADWT in open configuration and closed configuration. For simplification, the Savonius rotor in the open configuration is assumed to be a bluff body with the wake width equal to the diameter of the

Figure 9. Analytical and experimental comparison at 9 m/s.

bluff body, whereas in the closed configuration, it is treated as a nominal cylinder. The open configuration of ADWT is highly complex to model unless both the rotors are assumed to be aerodynamically independent which is not in reality. In the open configuration, the Savonius rotor dominates during the starting and low TSR with negligible torque from Darrieus rotor. Hence the current study emphasizes on the performance of the ADWT in closed configuration. Since the rpm of the Darrieus rotor is low in closed configuration, the cylinder is treated as non-rotating. The velocity deficit has been calculated for various diameter cylinder and the stream tubes that is encompassed by the wake width are modified with deficit velocity. The power coefficient is predicted against various TSR and the blades, rotor torque are calculated as a function of azimuthal angle. All the parameters that are investigated are compared for full and half-length Savonius rotor. The results indicate that the power loss between full and half-length Savonius is negligible. On various rotor diameter studies, 12.8% loss is reported for 2000 mm Savonius rotor. Hence an optimum configuration can be with 70 mm chord and full-length Savonius rotor of diameter 1200 mm from aerodynamic and structural perspective. The future work continues with scaling the proposed ADWT to 10 kW to monitor the performance in the field conditions. The dynamic variation of Savonius rotor will highlight the structural issues arising out, which has not been considered so far in the study. The 10 kW system is expected to shed light on the required actuation systems and the cost.

Acknowledgements

This research was supported by the National Research Foundation, Prime Minister’s Office, Singapore under its Energy Innovation Research Programme (EIRP Award No. NRF2013EWT-EIRP003-032: Efficient Low Flow Wind Turbine).

Nomenclature

$c$ Blade chord, m

${C}_{D}$ Drag coefficient

$C{F}_{N},C{{F}^{\prime}}_{N}$ Upwind and downwind elemental blade normal force coefficients

${C}_{L}$ Lift coefficient

${C}_{N}$ Normal force coefficient

${C}_{P}$ Darrieus rotor power coefficient

${C}_{P-D}$ ADWT rotor power coefficient

$C{F}_{T},C{{F}^{\prime}}_{T}$ Upwind and downwind elemental blade tangential force coefficients

C_{Q} Rotor torque coefficient

${C}_{T}$ Tangential force coefficient

$f$ Free-vortex frequency

${f}_{dw},{f}_{up}$ Downwind and upwind functions

${F}_{N},{{F}^{\prime}}_{N}$ Upwind and downwind blade normal force coefficients

${F}_{T},{{F}^{\prime}}_{T}$ Upwind and downwind blade tangential force coefficients

$H$ Half-height of the rotor, m

$l$ Blade length, m

$N$ Number of blades

$R,{R}_{D}$ Darrieus rotor radius at the equator, m

$R{e}_{t}$ Turbine Reynolds number

$R{e}_{b}$ Blade Reynolds number

$S$ Rotor swept area, m^{2}

$s$ Rotor solidity

$u,{u}^{\prime}$ Upwind and downwind interference factors

$V,{V}^{\prime}$ Upwind and downwind induced velocities, m/s

${V}_{e}$ Induced equilibrium velocity, m/s

${V}_{\infty}$ Wind velocity at the equator level, m/s

$W,{W}^{\prime}$ Upwind and downwind relative inflow velocity, m/s

$X,{X}^{\prime}$ Upwind and downwind local tip-speed ratio

${X}_{EO}$ Tip-speed ratio at the equator

$z$ Local turbine height, m

$\alpha ,{\alpha}^{\prime}$ Upwind and downwind local angle of attack, degree

$\beta $ Rotor maximum diameter/height ratio

$\delta $ Angle between the blade normal and the equatorial plane, degree

$\theta $ Azimuthal angle, degree

${\rho}_{\infty}$ Freestream density, kg/m^{3}

$\omega $ Turbine rotational speed, s^{-1}

$\eta $ r/R, Non-dimensional Cartesian coordinate

$\zeta $ z/H, Non-dimensional Cartesian coordinate

${A}_{s}$ Savonius turbine swept area, m^{2}

${C}_{P-S}$ Savonius power coefficient

${r}^{\prime}$ Savonius bucket radius, m

${H}_{s}$ Savonius rotor height, m

${N}_{s}$ Number of Savonius buckets

$Q$ Savonius turbine torque, N・m

$\phi $ Savonius Turbine tip-speed ratio

${P}_{r}$ Pressure on the returning Savonius bucket, Pa

${P}_{a}$ Pressure on the advancing Savonius bucket, Pa

${v}_{e-s}$ Equivalent velocity of Savonius rotor, m/s

${v}_{is}$ Absolute velocity of Savonius rotor, m/s

${v}_{\infty}$ Freestream wind speed, m/s

Conflicts of Interest

The authors declare no conflicts of interest.

[1] | Lu, X. and McElroy, M.B. (2017) Global Potential for Wind-Generated Electricity, Wind Energy Engineering. Elsevier, Amsterdam, 51-73. |

[2] |
Che, L., Zhang, X., Shahidehpour, M., Alabdulwahab, A. and Abusorrah, A. (2017) Optimal Interconnection Planning of Community Microgrids with Renewable Energy Sources. IEEE Transactions on Smart Grid, 8, 1054-1063. https://doi.org/10.1109/TSG.2015.2456834 |

[3] |
Li, Q., Maeda, T., Kamada, Y., Murata, J., Furukawa, K. and Yamamoto, M. (2015) Effect of Number of Blades on Aerodynamic Forces on a Straight-Bladed Vertical Axis Wind Turbine. Energy, 90, 784-795. https://doi.org/10.1016/j.energy.2015.07.115 |

[4] |
Worasinchai, S., Ingram, G.L. and Dominy, R.G. (2016) The Physics of H-Darrieus Turbine Starting Behavior. Journal of Engineering for Gas Turbines and Power, 138, Article ID: 062605. https://doi.org/10.1115/1.4031870 |

[5] | Kirke, B.K. (1998) Evaluation of Self-Starting Vertical Axis Wind Turbines for Stand-Alone Applications. Griffith University, Logan. |

[6] | Kumar, M., Surya, M.M.R., Sin, N.P. and Srikanth, N. (2017) Design and Experimental Investigation of Airfoil for Extruded Blades. International Journal of Advances in Agricultural and Environmental Engineering, 3, 2349-1523. |

[7] |
Karthikeya, B.R., Negi, P.S. and Srikanth, N. (2016) Wind Resource Assessment for Urban Renewable Energy Application in Singapore. Renewable Energy, 87, 403-414. https://doi.org/10.1016/j.renene.2015.10.010 |

[8] | Kumar, M., Surya, M.M.R. and Srikanth, N. (2017) On the Improvement of Starting Torque of Darrieus Wind Turbine with Trapped Vortex Airfoil. International Conference on Smart Grid and Smart Cities, Singapore, 23-26 July 2017, 120-125. |

[9] | Kumar, M., Surya, M.M.R. and Srikanth, N. (2017) Comparative CFD Analysis of Darrieus Wind Turbine with NTU-20-V and NACA0018 Airfoils. International Conference on Smart Grid and Smart Cities, Singapore, 23-26 July 2017, 108-114. |

[10] |
El-Samanoudy, M., Ghorab, A.A.E. and Youssef, S.Z. (2010) Effect of Some Design Parameters on the Performance of a Giromill Vertical Axis Wind Turbine. Ain Shams Engineering Journal, 1, 85-95. https://doi.org/10.1016/j.asej.2010.09.012 |

[11] |
Fujisawa, N. and Shibuya, S. (2001) Observations of Dynamic Stall on Darrieus Wind Turbine Blades. Journal of Wind Engineering and Industrial Aerodynamics, 89, 201-214. https://doi.org/10.1016/S0167-6105(00)00062-3 |

[12] |
Iio, S., Nakajima, M. and Ikeda, T. (2008) Performance of Double-Step Savonius Rotor for Environmentally Friendly Hydraulic Turbine. Journal of Fluid Science and Technology, 3, 410-419. https://doi.org/10.1299/jfst.3.410 |

[13] | Prandtl, L. (1935) The Mechanics of Viscous Fluids. Aerodynamic Theory, 3, 208. |

[14] |
Díaz, F., Gavaldà, J., Kawall, J.G., Keffer, J.F. and Giralt, F. (1983) Vortex Shedding from a Spinning Cylinder. The Physics of Fluids, 26, 3454-3460. https://doi.org/10.1063/1.864127 |

[15] | Schlichting, H., Gersten, K., Krause, E. and Oertel, H. (1955) Boundary-Layer Theory. Springer, Berlin. |

Journals Menu

Contact us

customer@scirp.org | |

+86 18163351462(WhatsApp) | |

1655362766 | |

Paper Publishing WeChat |

Copyright © 2023 by authors and Scientific Research Publishing Inc.

This work and the related PDF file are licensed under a Creative Commons Attribution 4.0 International License.