Hugo Carrillo

Axel Osses

Sergio Uribe

Cristóbal Bertoglio

Optimal Dual-VENC (ODV) Unwrapping in Phase-Contrast MRI

Keywords: Phase-Contrast MRI, dual-VENC, unwrapping

Dual-VENC strategies have been proposed to improve the velocity-to-noise ratio in Phase-Contrast MRI. However, they are based on aliasing-free high-VENC data. The aim of this work is hence to propose a dual-VENC velocity estimation method allowing high-VENC aliased data. For this purpose, we reformulate the phase-contrast velocity as a least squares estimator, providing a natural framework for including multiple encoding gradient measurements. By analyzing the mathematical properties of both single-and dual-VENC problems, we can justify theoretically high/low-VENC ratios such that the aliasing velocity can be minimized. The resulting reconstruction algorithm was assessed using three types of data: numerical, experimental and volunteers. In clinical practice, this method would allow shorter examination times by avoiding tedious adaptation of VENC values by repeated scans.

I. INTRODUCTION

Velocity-encoded Phase-Contrast MRI (PC-MRI) is a wellestablished method for measuring flow velocities, with several applications to quantitative analysis of cardiovascular pathologies [START_REF] Srichai | Cardiovascular applications of phase-contrast mri[END_REF]. The velocity-encoding magnetic gradients are set by the choice of the velocity encoding parameter, or VENC [START_REF] Dyverfeldt | 4d flow cardiovascular magnetic resonance consensus statement[END_REF]. It is well known that the velocity-to-noise-ratio (VNR) worsens when increasing the VENC. However, if VENC is set lower than the true velocity (which is unknown prior to the scan), velocity aliasing occurs. Moreover, even for VENC values slightly larger than the true velocity, velocity aliasing may occur due to measurement noise. These restrictions promote in clinical practice to acquire images at different VENCs, obligating the MRI operator to manually select the image for one specific VENC, while the aliased images are ignored and the time spent is squandered.

Corresponding author. Email: c.a.bertoglio@rug.nl. Address: Nijenborgh 9, 9747 AG Groningen, The Netherlands * C.B. and S.U. are both last authors.

Velocity aliasing is one of the main limitations for measuring complex features of blood flows, particularly, when high and low velocities are present in the same image, such as in heart, valvular and vascular malformations.

Then, VENC has to be set high, but as a consequence, low VNR is present in low velocity regions, for instance in recirculation regions in aneurisma or false lumen in dissections, to name a few. This leads to important inaccuracies when further analysis of the flow is performed [START_REF] Callaghan | Use of multi-velocity encoding 4d flow mri to improve quantification of flow patterns in the aorta[END_REF]. Aliasing is also problematic in many PC-MRI techniques, like Tissue Phase Mapping [START_REF] Petersen | Myocardial tissue phase mapping with cine phase-contrast mr imaging: regional wall motion ana lysis in healthy volunteers[END_REF] and Elastography [START_REF] Hirsch | Magnetic resonance elastography: physical background and medical applications[END_REF], where the motions magnitude vary across the regions of interest.

In order to reduce aliasing artefacts, unwrapping algorithms have been developed by assuming that the velocity field is smooth in space and/or time, see e.g. [START_REF] Loecher | Phase unwrapping in 4d mr flow with a 4d single-step laplacian algorithm[END_REF] and references therein. Nevertheless, they often fail when the aliased regions are large. Therefore, voxelwise dual-VENC strategies have been proposed, i.e. without any assumption of smoothness of the flow [START_REF] Lee | Three-point phase-contrast velocity measurements with increased velocity-to-noise ratio[END_REF], [START_REF] Nett | Four-dimensional phase contrast mri with accelerated dual velocity encoding[END_REF], [START_REF] Ha | Multi-venc acquisition of four-dimensional phase-contrast mri to improve precision of velocity field measurement[END_REF], [START_REF] Callaghan | Use of multi-velocity encoding 4d flow mri to improve quantification of flow patterns in the aorta[END_REF], [START_REF] Schnell | Accelerated dual-venc 4d flow mri for neurovascular applications[END_REF]. They have been based on unwrapping low-VENC data by using the high-VENC reconstruction, which is assumed aliasing-free. While actual approaches allow to improve the VNR with respect to a single high-VENC acquisition, they fail when the high-VENC data is aliased. Also, there is a lack of mathematical support for choosing the low-and high-VENCs. All of these issues limits the applicability of dual-VENC techniques, particularly when the peak velocities are uncertain.

Therefore, the aim of this work is to provide a mathematical framework to obtain aliasing-free velocity estimations from dual-VENC data, even when the both VENC acquisitions are aliased. The key is the least-squares formulation of the PC-MRI problem, whose mathematical properties allow to propose optimal combinations of VENCs to achieve this goal. We also present a numerical algorithm for dual-VENC reconstructions, which is successfully applied to numerical, experimental and volunteer data sets.

II. THEORY A. Classical PC-MRI

Assuming a constant velocity field, the usual starting point of classical PC-MRI is the model for the phase of the transverse magnetization at the echo-time [START_REF] Lee | Three-point phase-contrast velocity measurements with increased velocity-to-noise ratio[END_REF]:

ϕ G = ϕ 0 + ϑ G (1)
with ϕ 0 ∈ [0, 2π) the reference phase, and

ϑ G = ϑ G (u) = γuM 1 (G) (2)
the velocity dependent phase. Here, u ∈ R the flow velocity component parallel to the velocity-encoding gradient G = G(t) ∈ R, with t the encoding time, and M 1 (G) ∈ R the firstorder moment of G(t). The constat γ > 0 is the giromagnetic ratio.

From now on, we deal with different gradients G i with different amplitudes. Assuming that we have measured two phases ϕ G0 and ϕ G1 with G 0 = G 1 , the phase-contrast velocity is estimated by

u pc := ϕ G0 -ϕ G1 π VENC(G 0 , G 1 ), (3)
with

VENC(G 0 , G 1 ) = π γ(M 1 (G 0 ) -M 1 (G 1 ))
.

In the case that the true velocity |u true | ≤ |VENC|, then u pc = u true . But if |u true | > |VENC|, the phase difference exceeds ±π and aliasing occurs, i.e. u true = u pc . However, increasing the VENC decreases the VNR. Therefore, choosing the VENC parameter is an iterative manual process trying to set it as small as possible to maximize VNR and at the same time large enough to avoid aliasing.

B. Dual-VENC approaches

It is well known that for any VENC value, u true belongs to the set of infinite but numerable solutions of type

u pc + 2kVENC(G 0 , G 1 ), k ∈ Z. (4)
Therefore, it is natural to extend the velocity estimation problem such that k can be also estimated using additional encoding gradient measurements.

Assuming that now three measurements with gradients G 0 = 0 < G 1 < G 2 are available, two velocities at different VENC values can be reconstructed: the phase-contrast velocity u 1 at VENC 1 = VENC(G 1 , 0) and a set of velocities

u 2 + 2kVENC 2 at VENC 2 = VENC(G 2 , 0), with VENC 1 > VENC 2 , k ∈ Z.
Standard dual-VENC unwrapping strategies, see e.g. [START_REF] Schnell | Accelerated dual-venc 4d flow mri for neurovascular applications[END_REF], [START_REF] Lee | Three-point phase-contrast velocity measurements with increased velocity-to-noise ratio[END_REF], aim to find the correct low-VENC velocity from an un-aliased high-VENC velocity u 1 . Hence, an improved VNR should be achieved. Here, we will compare our new dual-VENC approach against the one from [START_REF] Schnell | Accelerated dual-venc 4d flow mri for neurovascular applications[END_REF], which is defined as:

u SDV =        u 2 + 2 • VENC 2 if 1 < D < 2 u 2 -2 • VENC 2 if -2 < D < -1 u 2 + 4 • VENC 2 if 3 < D < 4 u 2 -4 • VENC 2 if -4 < D < -3 with D = u 1 -u 2 ; 1 = 1.6 • VENC 2 ; 2 = 2.4 • VENC 2 ; 3 = 3.2 • VENC 2 ; 4 = 4.8 • VENC 2 .
In the reminder of this article, we will denote it as standard dual-VENC (SDV).

Note that the SDV reconstruction will be aliased if |VENC 1 | < |u true |. The new dual-VENC method based on our analysis will overcome this issue by optimally choosing both VENC 1 and VENC 2 based on a reformulation of the phase-contrast problem presented next.

C. Least-squares formulation of the single-VENC problem

For a given velocity encoding gradient G let us denote the measured phase of transverse magnetization by φG .

Assume now that we have available two measurements: a reference one for G = 0, and another for G = 0. We formulate the velocity reconstruction as a standard maximum-likelihood estimation problem from the phase measurements, by means of the least-squares function

J G (u) = 1 2 cos( θG ) -cos(ϑ G (u)) 2 + 1 2 sin( θG ) -sin(ϑ G (u)) 2 (5) = 1 -cos θG -ϑ G (u) (6)
with θG = φG -φ0 the "measured" velocity dependent phase.

Least-squares formulations have also been recently applied in the context of unwrapping methods using the information of contiguous voxels for various types of single-and dual-VENC acquisitions [START_REF] Loecher | Velocity reconstruction with nonconvex optimization for low-velocity-encoding phase-contrast mri[END_REF]. However, no analysis of their properties or potential for optimizing the VENC combinations was reported.

Figure 1 shows examples of the functions J G (u), for different gradients represented by VENC(G, 0). The synthetic measurements were generated with a unitary magnitude and the phases from Equation (1) using ϕ 0 = γBt E with B = 1.5 T , γ = 267.513e3 rad/T /ms, t E = 5 ms, a velocity u true = 1 m/s. It can be appreciated that the functions are periodic, with the period depending on the VENC, and also that the true velocity is a local minimum independent on the VENC. The following propositions proof these observations. Proposition 1: J G (u) is a periodic function with period 2VENC(G, 0).

Proof: It suffices to see that the cosine and sine are 2πperiodic functions, and

ϑ G (u + 2VENC(G, 0)) = γ(u + 2VENC(G, 0))M 1 (G) = γuM 1 (G) + 2π = ϑ G (u) + 2π
Proposition 2: The critical points u k of J G (u) are

u k = φG -φ0 γM 1 (G) + kVENC(G, 0) , k ∈ Z (7)
Proof: From ( 6) we see that

∂J G ∂u = -γM 1 (G) sin( θG -ϑ G ). (8)
At the critical points we must then have:

sin( θG -ϑ G ) = 0 ⇐⇒ ϑ G (u k ) = θG + kπ, k ∈ Z. (9)
Finally, using Equation ( 2) we obtain

u k = φG -φ0 γM 1 (G) + k π γM 1 (G) . ( 10
)
Proposition 3: At the critical points of J G (u), the second derivatives are given by

∂ 2 J G ∂u 2 (u k ) = C • (-1) k , k ∈ Z, C > 0.
Proof: Taking the derivative in (8) we obtain:

∂ 2 J G ∂u 2 (u k ) = γ 2 M 1 (G) 2 cos( θG (u k )-ϑ G u ) = C •(-1) k (11)
where the last equality holds due to Equation [START_REF] Ha | Multi-venc acquisition of four-dimensional phase-contrast mri to improve precision of velocity field measurement[END_REF].

In conclusion, we have just proved that Equation (4) corresponds to the local minima of the cost function J G by taking k as an even number in Equation [START_REF] Loecher | Velocity reconstruction with nonconvex optimization for low-velocity-encoding phase-contrast mri[END_REF].

It is also straightforward to show that the true velocity u true belongs to the set of local minima of J G when the measurements are noise-free. Indeed, in that case φG = φ0 + γM 1 (G)u true + 2kπ, and if we choose ϕ G (u true ) = φ0 + γM 1 (G)u true , then J G (u true ) = 0 from Equation (6).

D. The dual-VENC least squares problem

We assume now that we have measured the magnetization vector with three encoding gradients G 0 = 0 < G 1 < G 2 . We can then define the dual-VENC least squares sum function as:

J Σ (u) = 2 i=1 1 -cos θGi -ϑ Gi (u)
Figure 2 shows the single-and dual-VENC least-squares functions for different VENC combinations VENC 1 > VENC 2 = βVENC 1 , 0 < β < 1. Hence, the VENCs can be set in terms of VENC 1 and β. Note that VENC 1 is set lower than u true and is kept fixed in all plots, while β is variable. We can first observe that in all cases local and global minima are present in the dual-VENC functions J Σ (u). However, the true velocity is always a global minimum since it is a local and global minimum for each VENC, as shown in the previous section.

Remarkably, the periodicity of J Σ is now the least common multiplier (lcm) between the periodicity of the single-VENC functions, i.e. L Σ := lcm(2VENC 1 , 2VENC 2 ). As a consequence, if β is carefully chosen, e.g. in Figure 2(a) and 2(b), J Σ has a larger period than the original single-VENC functions, namely L Σ > 2VENC 1 . Therefore, even though VENC 1 , VENC 2 < |u true |, we can still distinguish u true from the other global minima since they have larger absolute values.

However, if we do not choose β well, e.g. as in Figures 2(c) and 2(d), then L Σ = 2VENC 1 and the global minima with smallest absolute value will not be u true if VENC 1 < u true and velocity aliasing occurs.

A general method for computing the aliasing limit is: for β = α/α 0 , with α, α 0 ∈ N the smallest possible values, then it is easy to verify that the periodicity of J Σ is L Σ = α2VENC 1 , since

L Σ = k 1 2VENC 1 = k 2 2βVENC 1 , k 1 , k 2 ∈ Z leading to k 1 = α, k 2 = α 0 .
Then, aliasing will occur when

||u true | -L Σ /2| < |u true |, i.e. VENC 1 < |u true |/α.

E. Choice of β

As shown in Table I, in the case without any measurement noise, to maximize the periodicity of J Σ one should choose VENC 2 ≈ VENC 1 , making the aliasing velocity very small, or for instance β = 0.7 or β = 0.55 as indicated in Table I.

However, the presence of noise deforms the dual-VENC functions, see Figure 3, since the noise is independent for each VENC. Therefore, local minima from both single-VENC cost functions that are not necessarily u true can get close to each other. Hence, there is an increased risk for u true not being global minima when α is large. In order to maximize the robustness to noise, the local minima of both single-VENC functions should be separated as much as possible. As shown in Figure 2(b), this is indeed the case for β = 0.66. For β = 0.75 this separation is less pronounced, however β = 0.75 would allow to lower the aliasing velocity if noise is low. In general, the optimal choice of β should be optimized to the SNR of the specific MRI scanner, but β = 0.66 is always the most robust to noise due to the largest separation between minima.

In the experiments, we will use these two values, β = 0.66 and β = 0.75. Additionally, in the experiments with numerical data, we will show the poor performance of β = 0.7 when noise is present.

F. The optimal dual-VENC (ODV) algorithm

Based on the considerations above, we now detail the ODV velocity estimation algorithm. For the given user-defined parameters VENC 1 and VENC 2 = βVENC 1 , 0 < β < 1:

1) Measure phases φGi for three gradients:

G 0 = 0 and G 1 , G 2 such that VENC(G 1 , 0) = VENC 1 and VENC(G 2 , 0) = VENC 2 . 2) Find the global minima u * k , k ∈ Z: u * k = argmin u ∈ [-u max , u max ] J Σ (u),
with u max = lcm(2VENC 1 , 2VENC 2 )/2. The estimated dual-VENC velocity corresponds to u * k with smallest absolute value.

III. METHODS

This section summarizes setups with three types of data: synthetic, phantom and volunteer. In all cases we applied the formula (3) for single-VENC and dual-VENC with both standard [START_REF] Schnell | Accelerated dual-venc 4d flow mri for neurovascular applications[END_REF] (SDV) and new ODV methods. For the ODV algorithm, the global minima was found using a sampling of the cost function J Σ with uniform spacing of the velocity of VENC 2 • 10 -3 , which was found to be small enough to avoid numerical artefacts in the global optimization.

A. Synthetic data

The reference phase is defined as ϕ 0 = γB 0 t E with B 0 = 1.5 T , γ = 267.513e3 rad/T /ms, t E = 5 ms. For the phases of the non-zero flow encoding gradients, we consider ϕ G1,2 = ϕ 0 + u true π/VENC 1,2 , with u true = 1m/s. Using these phases, reference magnetization measurements were built assuming a unitary magnitude. The estimation is shown in terms of VENC 1 and VENC 2 = βVENC 1 , with β = {0.66, 0.7, 0.75}.

We also compute estimations using magnetization measurements perturbed with an additive Gaussian noise with zeromean and standard deviation of 20% of the magnitude. We express these results in terms of mean estimated velocity for 2000 realizations of the noise and twice the standard deviation.

B. Phantom data

In order to preliminary assess the ODV we used a flow phantom that consisted of a rigid straight hose of 15mm internal diameter, 25mm external diameter. The hose was connected to a MRI-compatible flow pump (CardioFlow 5000 MR, Shelley Medical Imaging Technologies, London, ON, Canada) with a constant flow rate of 200 mL/s. The system was filled with a blood-mimicking fluid (40% distilled H2O, 60% Glycerol) and the set up was similar as in [START_REF] Urbina | Realistic aortic phantom to study hemodynamics using mri and cardiac catheterization in normal and aortic coarctation conditions[END_REF], [START_REF] Montalba | Variability of 4d flow parameters when subjected to changes in mri acquisition parameters using a realistic thoracic aortic phantom[END_REF]. The MRI data sets were acquired on a clinical 1.5T Philips Achieva scanner (Philips, Best, The Netherlands). The protocol consisted of through-plane PC-MRI sequence with a single cardiac phase due to constant flow rate. The scan parameters were: in-plane resolution was 1x1 mm with a slice thickness of 8 mm, 1 prospective cardiac phase, FA = 12 o , TR=9.2 ms, TE=4.9 ms, matrix size = (256,256).. The data was using non-symmetric pairs of encoding gradients with VENC = 150, 100, 70 cm/s with one surface coil. The acquisitions were performed using single-VENC protocols and the dual-VENC reconstructions were computed using only one of the zero-encoding gradients of the corresponding dual-VENC pair.

C. Volunteer data

Eight healthy volunteers underwent MRI in the same 1.5T Achieva scanner using a 5 elements cardiac coil. The protocol consisted of through-plane PC-MRI sequence perpendicular to the ascending aorta just above the valsalva sinus. We used several VENC values: 33.3, 37.5, 50, 66.7, 75, 100 and 150cm/s. These choices allow to generate dual-VENC reconstructions with both values of β = 0.66 and β = 0.75. The raw data was obtained and the reconstruction of each bipolar gradient was performed offline using matlab. Data from the multiple coils were combined using the method proposed in [START_REF] Bernstein | Reconstructions of phase contrast, phased array multicoil data[END_REF]. The data was acquired using the following scan parameters: in-plane resolution was 1x1 mm with a slice thickness of 8 mm, 25 cardiac phases using prospective ECG triggering, FA = 15 o , TR=5.5 ms, TE=3.7 ms, matrix size = (320, 232). Temporal resolution depended on the heart rate of the patients, varying between 35ms to 48 ms.

As in the panthom, the acquisitions were performed using single-VENC protocols. One issue with this approach is that the TE may be different depending in the scan setting, particularly may increase for low VENCs [START_REF] Bernstein | Minimizing te in moment-nulled or flow-encoded two-and three-dimensional gradientecho imaging[END_REF]. Since we use only the reference phase of VENC 1 , the value of the reference phase used in the dual-VENC reconstructions for VENC 2 was scaled by T E the echo times given by acquisitions with VENC 1 and VENC 2 , respectively. This is justify simply by the knowledge about the reference phase being proportional to T E [START_REF] Brown | Magnetic Resonance Imaging: Physical Principles and Sequence Design[END_REF].

IV. RESULTS

A. Synthetic data

Figure 4 shows the estimated velocity against VENC 1 without noise, confirming the unwrapping properties of both dual-VENC approaches: for SDV aliasing occurs when VENC 1 < u true , and for ODV when VENC 1 < u true /2, VENC 1 < u true /7 and VENC 1 < u true /3 with β = 0.66, β = 0.7 and β = 0.75, respectively.

Similar results for noisy measurements are presented in Figure 5, now including the aforementioned confidence interval. As one expects, the spread of the estimations are lower for β = 0.66. Moreover, in the single-VENC cases we confirm that aliasing starts even before the theoretical value due to the noise. This is also evident for SDV, while ODV is clearly more robust. We can also see that for ODV and β = 0.7 the confidence interval does not decrease uniformly with VENC 1 due to the nondesirable effect of overlapping of the single-VENC least squares functions mentioned in Section II-E. A similar, but less pronounced effect, occurs with β = 0.75. Therefore, in the real data acquisitions we continue using only β = 0.66 and β = 0.75.

B. Phantom data

The results for the phantom experiments are presented in Figure 6. The peak velocity in the tube is about 120 cm/s, what can be inferred from the single-VENC image with VENC 1 = 150. The wall of the tube can be distinguish as the noise ring separating the flow and the surrounding zerovelocity fluid. We first show the single-VENC PC-MRI, where aliasing for the two smaller VENCs can be clearly appreciated. We also confirm that SDV cannot handle the aliasing when both VENC values are lower than the true velocity, while ODV is able to sucessfully reconstruct un-aliased images from two aliased ones. C. Volunteers data Figure 7 presents the velocity profiles on the descending aorta for the different VENC combinations and different reconstruction methods for Volunteer 6. The figures for all the volunteers can be found in the Supplementary Material.

In all volunteers it is confirmed that ODV is the most robust method when decreasing the VENC, allowing to reconstruct velocities using lower VENCs than the true velocity, in contrast to SDV. Moreover, the theory is verified: aliasing is practically inexistent for (VENC 1 , VENC 2 ) = (50, 37.5) (β = 0.75), while aliasing always occurs at (50, 33.3) (β = 0.66). Indeed, the peak velocity is approximately 130 cm/s, for β = 0.75, VENC 1 = 50 > 130/3 ≈ 40, hence no aliasing appears. For β = 0.66, VENC 1 = 50 < 130/2 ≈ 65, hence aliasing appears. The actual noise level of the acquisition seems to not affect the performance of the ODV.

Figure 8 summarises the ODV results for all volunteers when varying the VENC. The error is computed in terms of the 2 -norm for the voxels inside the lumen, relative to the 2 -norm of the reference image (average of VENC 150 cm/s with 3 repetitions).

V. DISCUSSION In this work, we present a method for reconstructing velocities using dual-VENC images, for the first time in the literature when both single-VENC images are aliased. The main advantage of the method is that the true velocity does not need be know exactly in advance, since alias is allowed for both VENCs. All previous works have proposed to unwrap low-VENC images using high-VENC images without aliasing [START_REF] Lee | Three-point phase-contrast velocity measurements with increased velocity-to-noise ratio[END_REF], [START_REF] Nett | Four-dimensional phase contrast mri with accelerated dual velocity encoding[END_REF], [START_REF] Ha | Multi-venc acquisition of four-dimensional phase-contrast mri to improve precision of velocity field measurement[END_REF], [START_REF] Callaghan | Use of multi-velocity encoding 4d flow mri to improve quantification of flow patterns in the aorta[END_REF], [START_REF] Schnell | Accelerated dual-venc 4d flow mri for neurovascular applications[END_REF]. The theoretical findings are confirmed in real data sets from an experimental phantom and volunteers.

The choice of the VENC's ratio β = 0.66 is the most robust to noise, independent on the MRI scanner settings. However, for the volunteers scanned here, β = 0.75 works satisfactory and therefore it allows lower aliasing limits for the ODV estimations than β = 0.66, as given by the theory. Let us recall that β can be kept fixed (for instance, optimized once for typical scan settings), while the scanner user only needs to choose VENC 1 as in a single-VENC acquisition.

Note that unwrapping methods using contiguous voxelslike the ones from [START_REF] Loecher | Velocity reconstruction with nonconvex optimization for low-velocity-encoding phase-contrast mri[END_REF] -can be still applied after the estimation with ODV. The unwrapping would then probably perform better due to the larger periods of the candidate solutions, e.g. L Σ = 6VENC 1 for β = 0.75 and L Σ = 4VENC 1 for β = 0.667.

Concerning the limitations of our study, the method was not assessed in patients, only in volunteers. It is well known that dual-venc approaches (as any other cardiovascular MRI sequences) are challenging due to variabilities during the experiment (not only measurment noise) [START_REF] Lee | Three-point phase-contrast velocity measurements with increased velocity-to-noise ratio[END_REF], such as cardiac rhythm changes and subjects' motion. However, this variability will impact in similar manner the standard dual-VENC approach as well as the method proposed here. Another limitation is that data acquisition was performed for the two VENCs in a serial fashion, and therefore MRI scan protocols tailored to the ODV reconstructions have to be developed yet. This could be also done by including k-space undersampling techniques as in [START_REF] Nett | Four-dimensional phase contrast mri with accelerated dual velocity encoding[END_REF], what would allow dual-VENC protocols comparable in scan time to single-VENC ones, what is of high interest for the application of ODV to 4Dflow. Moreover, as in standard PC-MRI, there is the implicit assumption that the velocity is constant in space and time and therefore, neither the single-nor the dual-VENC approaches count for effects like dephasing of spins and turbulence.

VI. CONCLUSION

We present a robust method for estimating velocities from dual-VENC data in PC-MRI. The main contribution of this work is that both a theoretical and an extensive empirical analysis was carried out, turning out that there are highand low-VENC combinations that can considerably reduce the aliasing issues. For example, in the volunteer data the ODV allows to choose the high-VENC up to a third of the maximal velocity. In clinical practice, the scanner operator has only to choose a single expected velocity, as for standard single-VENC PC-MRI. Then, the low-VENC value can be automatically fixed by the scanner in terms of the high-VENC. Moreover, the reconstruction method is simple enough to be implemented directly in the MRI scanner. Next steps are to assess the ODV in cases with high velocity variability, like stenotic vessels or valves, and 4Dflow, and application to other phase-contrast techniques, like elastography.

  Fig. 1. cost functions J G (u) for utrue and two VENC values.

Fig. 2 .

 2 Fig. 2. Cost functions J G (u) and J Σ (u) for different VENC 1 , VENC 2 = βVENC 1 .

Fig. 3 .

 3 Fig. 3. Cost functions J G (u) for different pair of values of VENC and the sum cost function J Σ with noisy magnetization measurements (standard deviation 20% of magnitude).

Fig. 4 .

 4 Fig. 4. Synthetic data (noise-free): single-and dual-VENC.

Fig. 5 .

 5 Fig. 5. Synthetic data (20% noise): single-and dual-VENC.

Fig. 6 .

 6 Fig. 6. Phantom data: single-(PC) and dual-VENC.

  ACKNOWLEDGMENT H.C. was supported by CONICYT-Doctorado nacional 21151645. C.B and A.O. acknowledge the funding of Conicyt

(

Fig. 7 .Fig. 8 .

 78 Fig. 7. Volunteer 6. First row: single-VENC PC-MRI. Second row: SDV. Third row: ODV. Velocities are coloured as in Figure 6.

  Table I gives examples of VENC 1 such that the global minimum of J Σ with lowest magnitude corresponds to u true depending on β.

	β	0.95	0.9	0.75	0.7	0.66	0.55	0.5
	α	19	9	3	7	2	11	1

TABLE I .

 I EXAMPLES OF ALIASING LIMITS FOR DECREASING VALUES OF β. ODV METHOD ALLOWS ALIASING-FREE ESTIMATION IF VENC

1 > |utrue|/α.

Basal Program PFB-03, and A.O. Fondecyt 1151512 and Fondap CR2-1511009. A.O and S.U. acknowledge funding from Millenium Science Initiative of the Ministry of Economy, Development and Tourism, grant Nucleus for Cardiovascular Magnetic Resonance. S.U. also thanks to Fondecyt 1181057.


David Marten
email: david.marten@tu-berlin.de

George Pechlivanoglou

Christian Nayeri

Christian Oliver Paschereit

Nonlinear Lifting Line Theory Applied To Vertical Axis Wind Turbines: Development of a Practical Design Tool

Keywords:

Recently a new interest in vertical axis wind turbine (VAWT) technology is fueled by research on floating support structures for large scale offshore wind energy application. For the application on floating structures at multi megawatt size, the VAWT concept may offer distinct advantages over the conventional horizontal axis wind turbine (HAWT) design. As an example VAWT turbines are better suited for upscaling and, at multi megawatt size, the problem of periodic fatigue cycles reduces significantly due to a very low rotational speed. Additionally, the possibility to store the transmission and electricity generation system at the bottom, compared to the tower top as in a HAWT, can lead to a considerable reduction of material logistics costs. However, as most VAWT research stalled in the mid 90's [1], no established and sophisticated tools to investigate this concept further exist today. Due to the complex interaction between unsteady aerodynamics and movement of the floating structure fully coupled simulation tools, modelling both aero-and structural dynamics are needed. A nonlinear Lifting Line Free Vortex Wake code was recently integrated into the open source wind turbine simulation suite QBlade. This paper describes some of the necessary adaptions of the algorithm, which differentiates it from the usual application in HAWT simulations. A focus is set on achieving a high robustness and computational efficiency. A short validation study compares simulation results with those of a U-RANS and a Double Multiple Streamtube (DMS) simulation.

INTRODUCTION

The aerodynamics and the flow field around a VAWT are very complex. During one revolution of the rotor the angle of attack, and thus the blade loads, are constantly changing. Additionally, at low tip speed ratios (TSR), the blades experience dynamic stall at certain azimuthal positions. This inherent unsteadiness during operation is the main reason why no practical, efficient and robust design tool exists for VAWT, which is comparable to the very popular blade element momentum (BEM) method, widely used for HAWT design and certification. Paraschivoiu [START_REF] Paraschivoiu | Wind Turbine Design -With Emphasis on Darrieus Concept[END_REF] developed the double multiple stream tube algorithm (DMS), as an extension of the BEM method, for VAWT. In contrast to the BEM the DMS algorithm employs two actuator discs, one for the upstream and one for the downstream half circle of rotation. This method is used to predict the performance of VAWT in the early design stage, however for high blade loading, due to high solidities or a high TSR [START_REF] Mcintosh | Convergence Failure and Stall Hysteresis in Actuator-Disk Momentum Models Applied to Vertical Axis Wind Turbines[END_REF] convergence problems arise. One problem is that momentum balance based codes assume a steady state flow field and this assumptions is even more violated for VAWT than for HAWT rotors. On the other hand Reynoldsaveraged Navier-Stokes (RANS) or Large Eddy (LES) CFD simulations are still too computationally expensive to be used as a design tool in coupled dynamics simulations. Being of intermediate complexity and physical accuracy Lifting Line Theory (LLT) based codes, combined with a free vortex wake model, are a promising candidate for a capable VAWT design tool. The Lifting Line Theory captures the unsteady blade loading, whereas the complex wake structure (see Figure 1) is modeled through free evolving vortex elements. LLT based codes have been successfully applied to VAWT simulation (as by Ferreira [START_REF] Simão | The near wake of the VAWT, 2D and 3D views of the VAWT aerodynamics[END_REF]), but up to date no such code, dedicated to VAWT, is commercially or publicly available. This paper describes a newly developed, efficient and robust LLT based design and simulation tool for VAWT. An unsteady LLT simulation module for VAWT is integrated with the open source wind turbine simulation tool QBlade [START_REF] Marten | Implementation, Optimization and Validation of a Nonlinear Lifting Line Free Vortex Wake Module Within the Wind Turbine Simulation Code QBlade[END_REF]. As a result of this integration, advanced simulations can be setup using the airfoil, polar, turbulent wind field and blade database of QBlade. Additional aerodynamic models for the tower deficit and ground effects have been implemented (details in [START_REF] Marten | QBlade Guidelines v0.9[END_REF]). Furthermore the turbine model is capable of prescribed 6 degree of freedom movement enabling an easy coupling to other dynamics simulation modules. The focus of this paper is on the specific differences and treatments that distinguish the implementation of a LLT code for VAWT from an implementation for HAWT rotors and how these affect its robustness, efficiency and accuracy. Moreover the performance of the code will be demonstrated in a validation study, comparing the obtained results to published data from a CFD and a DMS simulation. The resulting software is freely distributed under an open source license. The code can be found under: sourceforge.net/projects/qblade.

LLT MODIFICATIONS FOR VAWT SIMULATIONS

Recently a nonlinear lifting line free vortex wake algorithm has been implemented in QBlade. The LLT is a very general simulation method that can be used to model both HAWT and VAWT rotors with an almost identical formulation. Blade forces are calculated using the lifting line formulation where the circulation is modeled through distributed vortices, located at the quarter chord positions of the blade (see Figure 2). The circulation is calculated from tabulated airfoil data and the blade inflow angles that are computed under the influence of the free wake induction. The wake is described by freely floating vortex elements that are shed from the trailing edge during every time step.

Compared to BEM based methods which model the wake with a momentum balance across the rotor disc this is a much more physically sound approach that relies on far less assumptions. As the history of the flow is always contained in the wake, in the form of previously shed wake elements, this method is especially suited to model transient turbine behavior, such as start-or ramp up cases or any situation involving relative motion of the whole rotor. Especially in these transient cases a significant improvement, compared to BEM based models, can be expected. An additional advantage is that, due to the explicit wake modeling, the simulation results also contain the velocity distribution around the rotor, allowing to extract wake deficit profiles or induced velocity distributions around the rotor (Figure 1). Because the LLT is also based on far less assumptions than momentum balance based methods, generally a far lower number of empirical corrections are required. As an example: the effect of blade tip loss is already included in the LLT formulation whereas momentum balance based codes rely on a tip loss correction. Furthermore, the convergence problems that momentum balance based methods are often facing are not an issue with LLT based codes as they can even model situations of flow reversal within the rotor. This robustness, combined with a reasonable computational cost (compared to CFD methods), makes the LLT a very promising candidate as a universal design method for VAWT.

More general details about the implemented LLT method can be found in the report of van Garrel [START_REF] Van Garrel | Development of a wind turbine aerodynamics simulation module[END_REF], details about the implementation in QBlade have been published in [START_REF] Marten | Implementation, Optimization and Validation of a Nonlinear Lifting Line Free Vortex Wake Module Within the Wind Turbine Simulation Code QBlade[END_REF]. The aim of this paper is to give some details about the specifics that differentiate the implementation for VAWT simulation from that for HAWT simulation.

Calculation of relative blade velocities and azimuthal discretization

During every time step an angle of attack (AoA) needs to be calculated at every blade section. As this AoA determines the lift and drag coefficients, stored in data tables, and thus the forces acting on the respective blade section, it has a very large effect on the predicted performance of the rotor. The AoA is evaluated at the quarter chord position of each blade panel from the total velocity that the blade experiences. The total blade velocity is a superposition of inflow velocity, induced velocity, due to wake and rotor induction and a relative velocity component due to the blades' movement.

In most implementations of the LLT for HAWT applications the velocity component due to blade movement is computed via a simple derivative of the blades' position in time [START_REF] Van Garrel | Development of a wind turbine aerodynamics simulation module[END_REF][START_REF] Kloosterman | Development of the Near Wake behind a Horizontal Axis Wind Turbine[END_REF]. This practice is usually accurate enough for HAWT, where a 10° azimuthal discretization is usually recommended [START_REF] Abedi | Development of Vortex Filament Method for Aerodynamic Loads on Rotor Blades[END_REF]. The error between the actual tangential velocity and the approximation through the derivation does not have a very large effect on evaluated AoA as the circular path of the blade is usually perpendicular to the inflow. For VAWT however the inflow velocity is orthogonal to the axis of rotation and thus the tangential velocity vector due to blade movement lies in plane with the inflow. In this case the error that is introduced due to the first order derivative has a very large effect on the computed AoA (Figure 3) and using the derivation procedure makes the simulation very sensitive to the chosen azimuthal discretization. The tangential component is constructed via a unit tangent vector and the known angular velocity magnitude. The component describing all other blade movements is then calculated by derivation. From this component the part of the blades' velocity due to rotation of the rotor has to be subtracted, as it is already included in the, now separate, tangential component. This is done by reconstructing the blades' angular position of the last time step for the current rotor axis and using this position in the derivative. This treatment makes the simulation results very insensitive to the chosen azimuthal discretization and reasonable results can be obtained with an azimuthal discretization of up to 20° (see Figure 4). The advantage of this velocity decomposition, when compared to higher order time discretization, is the significantly lower computational cost.

Wake treatment

The computational cost of the LLT method is situated between the cost of momentum balances based methods and those of CFD. To be used in certification test like scenarios, where several million aerodynamic time steps are required, a high efficiency is of utmost importance. The driving factor for the computational cost of the LLT method is the number of free vortex elements, contained in the wake. To simulate the free evolution of the wake, each wake element has to be convected under the influence of induction from all other wake elements thus the computational time needed to run a simulation scales roughly linearly with the squared number of free vortex elements. As new wake elements are created at the blades trailing edge during every time step the computational cost of the LLT grows exponentially over time when no measures are taken. To enable the implementation of several strategies that reduce the number of wake elements an object orient attachment-detachment scheme has been implemented. This is an important feature of the code as it allows to remove wake elements from the simulation domain, based on arbitrary conditions, without breaking the connectivity of the complex wake mesh.

Wake truncation

One practical approach to limit the number of free wake elements to a fixed size is to simply remove all wake elements from the domain after they have reached a certain age (example in Figure 5). For this treatment the vortex age is made dimensionless with the time the rotor needs for a full revolution, thus the wake age, as it is used here, represents the age of a vortex in full revolutions of the rotor. To demonstrate the effect of this wake truncation on the resulting accuracy of the simulation, several test cases for exemplary tip speed ratios and a range of the parameter governing the maximum vortex age have been conducted. Figure 6 shows the relative error of the resulting averaged power coefficient compared to a baseline simulation, over different values of the wake age. The three values for tip speed ratio that have been chosen represent low (TSR: 2), optimal (TSR: 3.2) and high (TSR: 4) induction cases. The comparison shows that a wake length of around 8 rotor revolutions is an acceptable value, even for tip speed ratios, as the relative error due to the wake truncation drops below 1% for all simulated cases. The effect of wake truncation on the simulation accuracy

Wake reduction

When free vortex wake methods are applied to simulate HAWT it is common practice to reduce the amount of free vortices by lumping the trailing vortices into concentrated tip vortices and removing the shed vorticity from the wake. This usually takes place at some point sufficiently far away downstream from the rotor disc (as described in [START_REF] Montgomerie | Vortex Model for Wind Turbine Loads and Performance Evaluation[END_REF]). As the trailing vorticity usually is dominant in the wake of a HAWT, the removal of the shed vorticity only has very little effect on the simulation accuracy. Additionally, because the trailing vorticity is conserved, as it is lumped, this is a very effective means of reducing the number of wake elements drastically while still preserving a detailed wake shape several rotor diameters downstream of the rotor plane.

Due to the ever changing angle of attack of a VAWT, and the resulting change in blade circulation, neither the trailing nor the shed vorticity is dominant in the flow field. Both types of vorticity play an equal part in the description of the flow field and thus the strategy as applied to HAWT cannot be applied to VAWT simulations. The strategy that is applied instead is based on the assumption that in every time step there are vortices being shed from the trailing edge that have a more or less negligible circulation. Usually, for straight blade geometries, this is the trailing vorticity, shed from the inner blade regions, or the shed vorticity being shed when the blade is at an azimuthal position with very light loading. Based on this assumption a percentage of vortices is removed right after it is being shed. The amount of vorticity to be removed is a user input. Based on this input the vortices with the lowest absolute circulation, created at the last time step, are removed until the number of vortices that was specified in the user input has been met. To demonstrate the effect that this wake reduction has on the accuracy of the simulation the relative error of the computed power coefficient (averaged over one rotor revolution) is plotted over the percentage of removed wake elements in Figure 7. It can be seen that, regardless of the simulated tip speed ratio, up to 50% (effectively reducing simulation time by a factor of 4) of the vortices can be removed in this way without introducing an error above 1%. Figure 8 shows an exemplary visualization of a wake where 30% of the vortices have been removed.

VALIDATION / CODE COMPARISON

This sections compares the quasi steady state performance predictions of the LLT method with those of a DMST model, already implemented in QBlade [START_REF] Marten | Development and application of a simulation tool for vertical and horizontal axis wind turbines[END_REF], and those of a 2D U-RANS simulation carried out at the University of Florence (described in [12; 13]). The rotor used in this comparison is a straight bladed H-Darrieus rotor, simulated in a one bladed configuration (Table 1, Figure 8). To compare the global results of the simulation methods graphs of the averaged power coefficient over the tip speed ratio are shown in Figure 10. Furthermore the local azimuthal distributions of the generated blade torque at three tip speed ratios, for a low, high and optimal induction case, are compared. The LLT simulation obtains the power curve via a ramp up simulation, varying the rotational speed at a constant 8m/s inflow. Adaptive time stepping is used to fix the azimuthal discretization for all simulated cases to and angle of 5°.The main simulation parameters for the LLT are listed in Table 2. The lift and drag polars used in the DMS and the LLT simulations are XFOIL [START_REF] Drela | Viscous-Inviscid Analysis of Transonic and Low Reynolds Number Airfoils[END_REF] generated polars (NCRIT = 5, forced transition at the leading edge) in the Reynolds number range between 5e4 and 1e6 (see Figure 9). Both codes interpolate within this Reynolds number range at every time step. The polars were extrapolated with the Viterna method, with a CD90 value of 1.85. Dynamic stall effects are not modelled in the LLT & DMS simulations. Figure 11 shows the azimuthal variation of the torque coefficient at a relatively low TSR of 1.67. The distributions differ greatly between URANS and the two engineering methods. Due to the low TSR of 1.67, the AoA varies between +-35° and the blade experiences deep stall during every rotation. Because no dynamic stall model is implemented in neither the LLT nor the DMS code these differences are to be expected. Again, it can be seen that the DMS tip loss model has a large effect on the predicted torque. Bearing in mind that the set of polars used is plain XFoil data created without any special modifications or a dynamic stall model, the agreement for the torque coefficient in Figure 12 is surprisingly good between the URANS and the LLT code.

The DMS code over-and under predicts the URANS results and also predicts the torque maxima at slightly lower azimuthal angles. The simulated TSR of 3.34 is near the maximum CP value, at this operational point the AoA varies between +15° and -10°, thus dynamic stall effects are expected to play a less significant role. Torque coefficient over azimuthal angle @ 300rpm / TSR: 3.34, 1 bladed configuration Figure 13 compares the torque profiles for a high induction case above the optimal TSR. The AoA at this operational point varies between +8° and -7°and dynamic stall effects play a minor role now. As a result the LLT prediction is slightly improved. The DMS results show the same trend for over and under prediction of torque dependent on the use of the tip loss model.

Figure 13. Torque coefficient over azimuthal angle @ 400rpm / TSR: 4.45, 1 bladed configuration

Runtime Information

As an estimate for the real life performance of the LLT: With the settings from Table 2, on an Intel® Xeon® E3-1230V3, with 3.3GHz and 4 physical (+4 virtual) cores the simulation takes 137 seconds for 10 revolutions of the rotor (or 720 time steps). At these settings the wake is made up of roughly 7000 free wake elements. If, for the same settings the azimuthal discretization is increased to 10°, the simulation only takes 22 seconds (10 revolutions or 360 time steps) with roughly 3500 free wake elements present in the domain.

Recently, after these simulations have been performed, the OpenCL [START_REF] Stone | OpenCL: A Parallel Programming Standard for Heterogeneous Computing Systems[END_REF] framework was implemented to evaluate the Biot-Savart equation in parallel on GPUs. This massive parallelization decreases the needed computing time by a factor of up to 30, as demonstrated in [START_REF] Marten | Integration of an Unsteady Nonlinear Lifting Line Free Vortex Wake Algorithm in a Wind Turbine Design Framework[END_REF].

CONCLUSION AND OUTLOOK

A Lifting Line Free Vortex Wake code, optimized for VAWT simulations, was implemented in the wind turbine design suite QBlade. During the implementation special emphasis was placed on achieving a high robustness and computational efficiency of the code and some of the measures taken have been discussed in the previous sections. A validation and comparison to a URANS and DMS simulations shows that the LLT code generates reasonable results for "steady state" performance simulations. However it was also found that the code is lacking a dynamic stall model (by now a Beddoes-Leishman type dynamic stall model has been implemented and validated, see [17; 18]). Furthermore the impact of polar data quality and different polar extrapolation techniques on the simulation results needs to be further investigated to increase the confidence of simulation results. However the present LLT model is a promising starting point towards a fully coupled aero-elastic simulation environment for VAWT.

Figure 1 .

 1 Figure 1. Left: Detail of free evolving VAWT wake, right: velocity distributions for high induction / high tip speed ratio case

Figure 2 .

 2 Figure 2. Geometry of a blade panel, position of lifting line, shed and trailing vortex elements

Figure 3 .

 3 Figure 3. Time derivative of position vector vs tangent vector for a VAWT

Figure 4 .

 4 Figure 4. Power coefficient over one rotation for six different azimuthal discretizations 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 time [s] Power Coefficient Cp

Figure 5 .

 5 Figure 5. Top view of single blade wake at a TSR of 4 with three different values for wake truncation. From top to bottom: 4, 7, and 10

Figure 6 .

 6 Figure 6. The effect of wake truncation on the simulation accuracy

Figure 7 .Figure 8 .

 78 Figure 7. The effect of wake reduction on the simulation accuracy

Figure 9 .

 9 Figure 9. Lift polars used in the DMS & LLT simulations

Figure 10 .

 10 Figure 10. CP over TSR curves, 1 bladed configuration

Figure 11 .

 11 Figure 11. Torque coefficient over azimuthal angle @ 150rpm / TSR: 1.67, 1 bladed configuration

Figure 12 .

 12 Figure 12. Torque coefficient over azimuthal angle @ 300rpm / TSR: 3.34, 1 bladed configuration

[ 1 ]

 1 H. J. Sutherland, D. E. Berg and T. D. Ashwill. A Retrospective of VAWT Technology. Sandia Report,

Table 1 .

 1 H-Darrieus reference rotor

	Blade Number	1
	Airfoil	NACA0018
	Radius	0.85m
	Chord	0.246
	Solidity	0.144

Table 2 .

 2 LLT Simulation parameters

	Fixed Azimuthal discretization	5°
	Wake truncation	8 revolutions
	Wake reduction	0.3
	Number of blade elements	8 (sinusoidal spacing)
	Inflow velocity	8 m/s


Irfan Ahmed
email: irfan.ahmed@uni-kassel.de

Martin Lawerenz

Development and Numerical Investigation of an Actuator Concept for Shape-Adaptive Airfoil Profile for Horizontal Axis Wind Turbine Rotor

Keywords:

To cope with the increase in demand on renewable energy, the modern horizontal axis wind turbines have grown in dimension and have already reached a rotor diameter of 150 meter. With the introduction of offshore wind parks, the door to extract multi-megawatt power from wind is opened and has set a trend of continuing increase in wind turbine rotor diameter. This imposes a trade-off between the energy production and the controllability of future wind turbines. The increase in rotor diameter results in higher torque requirement on the active pitch/stall-control motor. Moreover, it is difficult to achieve optimal inflow condition over the entire span length for rigid blades. This results in high structural loading, which in turn reduces the life span of the blades. In light of this situation, smart rotor blades with local airfoil form-adaption promises to be an alternative and effective control scheme for future wind turbines. The authors present the results of investigations on active airfoil form-adaption concept for wind turbine rotor blades. A novel actuator concept is developed for incorporating airfoil trailing edge form-adaption in wind turbine application. Structural analysis of the shape morphing system is carried out and the aerodynamic performance of the shape-adaptive airfoil is presented.

NOMENCLATURE

S XX normal stress [M L -1 T -2 ] S Y Y normal stress [M L -1 T -2 ] S ZZ normal stress [M L -1 T -2 ] S XY shear stress [M L -1 T -2 ] S Y Z shear stress [M L -1 T -2 ] S ZX shear stress [M L -1 T -2 ] {f } global force vector [M L 2 T -2 ] [M ] global mass matrix [M ] [K] global stiffness matrix [M T -2 ] U local displacement [L]

INTRODUCTION

Modern wind turbines run mostly at variable rotational speed and use collective pitch control to optimize the energy yield and blade loading. The larger swept area in future wind turbines means the blades would be more susceptive to changing inflow conditions and wind shear. In light of this situation, collective or even individual pitch control alone would not be able to fully tame the radial and azimuthal gradient of the inflow. To compensate, future mega-watt rated wind turbines would require a faster yawing and pitching mechanism which would lead to increased loading on the entire structure. An alternative solution is the introduction of shape-morphing airfoil sections which would allow faster and more detailed flow control over different span-wise locations of the rotor blades.

As the idea of employing shape-adaption technique in wind turbine blades is still in its early stage, the bulk of investigations [START_REF] Wang | Preliminary Investigation of the Active Flow Control Benefits on Wind Turbine Blades[END_REF][START_REF] Debiasi | Shape Change of the Upper Surface of an Airfoil by Macro Fiber Composite Actuators[END_REF][START_REF] Ragheb | Multi-Element Airfoil Configurations for Wind Turbines[END_REF][START_REF] Peter Bjørn Andersen | Advanced Load Alleviation for Wind Turbines using Adaptive Trailing Edge Flaps: Sensoring and Control[END_REF][START_REF] Hoeijmakers | Plasma Actuators for Active Flow Control on Wind Turbine Blades[END_REF]] are aimed at understanding the effectiveness and the performance of the flow-control mechanisms for wind turbine applications. The development of dedicated actuator concepts for shape-adaptive wind turbine blades is still in its early stage. A review of the actuator concepts and the state of the art for morphing profiles in wind turbine applications can be found in [START_REF] Barlas | Closed-loop Control Wind Tunnel Tests on an Adaptive Wind Turbine Blade for Load Reduction[END_REF]. Recently [START_REF] Meldgaard | Wind Turbine Blade With Deflectable Flaps[END_REF] have developed pneumatic actuator for wind turbine blade with deflectable flap. [START_REF] Rebsdorf | Wind Turbine Blades with Cambering Flaps[END_REF] developed a pneumatic actuator for flaps to change the airfoil camber. A recent work by [START_REF] Wang | Preliminary Investigation of the Active Flow Control Benefits on Wind Turbine Blades[END_REF] describes the advantages in using active flow control techniques in wind turbine applications. Plasma actuators are tested by [START_REF] Hoeijmakers | Plasma Actuators for Active Flow Control on Wind Turbine Blades[END_REF] as active flow-control mechanism for the wind turbines.

The results showed an improvement of lift and delay of flow separation for low Reynolds numbered flow. [START_REF] Debiasi | Shape Change of the Upper Surface of an Airfoil by Macro Fiber Composite Actuators[END_REF] used macro-fiber composite actuators to adapt the airfoil upper surface during operation. [START_REF] Ragheb | Multi-Element Airfoil Configurations for Wind Turbines[END_REF] have investigated the performance optimization of the thick airfoil used near the rotor hub using various combinations of aerodynamic elements. [START_REF] Peter Bjørn Andersen | Advanced Load Alleviation for Wind Turbines using Adaptive Trailing Edge Flaps: Sensoring and Control[END_REF] investigated the use of deformable trailing edge in reducing the blade loading. [START_REF] Zuteck | Adaptive Blade Concept Assessment: Curved Planform Induced Twist Investigation[END_REF] investigated the concept and material requirement of passive twist for wind turbine blades. In scope of this investigation, active trailing-edge camber adaption of airfoil is selected as the shape adaption strategy to develop the actuator system. This allows the performance optimization of the airfoil profile for fixed incidence angle as well as performance evaluation of the morphed profile for varying inflow conditions.

ACTUATOR CONCEPT

This work focuses on the development of an actuator to actively adapt the shape of the airfoil trailing edge section of horizontal axis wind turbine rotor. In the following sections, the strategy and the design constraints of the shape adaption and the actuation concept are presented.

Constraints of Shape Adaption

In a first step, a number of criteria is defined for the selection of the actuator for wind turbine application. These boundary conditions are specified by the different aspects of the wind turbine life cycle, spanning from the production process to the aerodynamic and structural performance requirements of the components. The boundary conditions are given in table 1 . These boundary conditions predefine the concept of the shape adaption. The aerodynamic requirements define the shape of the morphed airfoil, whereas the structural requirements define the material of the shape adaptive airfoil section and the actuator system.

Aerodynamic Shape of the Adaptive Airfoil

In the next step, the target slope of the camber line at the trailing edge for the shape-adaptive airfoil is selected.

Figure 1 shows schematic of flow-turning over an airfoil along with the velocity and angle definition used in this work. The initial form of the airfoil is depicted in blue and the final bent-shape is shown in green. The Euler turbine equation states

a = (c c,2 -c c,1 )u (1) a = (tanβ 2 -tanβ 1 )c m,1 u 1 (2)
For a constant-incident flow, the specific work a changes with the change in the angle ∆β 2 , which corresponds to the change in slope of the camber line at the trailing edge for the shape-adaptive airfoil. The specific work done by a deformed airfoil is given as

a * = (tanβ * 2 -tanβ 1 )c m,1 u 1 (3) a * = (tan(β 2 + ∆β 2 ) -tanβ 1 )c m,1 u 1 (4) α j β j β 1 β 2 c j c m,j c c,j ∆β 2 i u j w j w 1 w 2 Figure 1.
Flow Turning over Shape-Adaptive Airfoil Hence, a target bending of 8 • of the camber line is selected for the implementation in wind turbine application.

Actuation Concept

The actuation concept for incorporating a bending on the airfoil profile skin is depicted in the figure 2. The underlying concept is to use inward protruding beam elements on the inner side of the airfoil profile skin. The concept is similar to the stringer principal [START_REF] Peter | Erzielung einer optimierten Tragflügelverwolbung durch den Einsatz formvariabler Klappenstrukturen[END_REF]. Bending moment applied on the beams will deform the airfoil skins. The skins of the shape-adaptive airfoil section are modeled to be made of elastic material.

METHODICS

In the following sections, the core strategy of the design and analysis of the shape-adaptive airfoil profile are presented.

FSI Simulation

The interaction of fluid and structure is a considerable phenomenon in wind turbine flow. The modern wind turbine blades with sub 100 m diameter experience substantial loading owing to the aerodynamic forces. The design aspects viz. size, material, inflow velocity, design tip speed ratio, cone angle, control schemes each plays a significant part in the blade loading. The aeroelasticity of the wind turbine blade and tower is acknowledged in the scientific community and design tools are dedicated to model the aeroelastic behavior of the wind turbine flow [START_REF] Simms | NREL unsteady aerodynamics experiment in the NASA-Ames wind tunnel: a comparison of predictions to measurements[END_REF][START_REF] Rasmussen | Present status of aeroelasticity of wind turbines[END_REF][START_REF] Otto | State of the art in wind turbine aerodynamics and aeroelasticity[END_REF][START_REF] Hansen | Aeroelastic stability analysis of wind turbines using an eigenvalue approach[END_REF].

The majority of these codes resolves the aeroelasticity in a global manner where the basic beam equation is use to model the structure under the aerodynamic loading. From numerical point of view, the modeling of the fluid-structure interaction can be subdivided in the subtasks, a) resolving the fluid field variables, b) resolving the structural response to the body forces acting owing to the fluid flow, c) transferring the deformation of the structural mesh in the fluid model to incorporate the movement of the wetted body. A fully resolved monolithic approach to employ a system of aggregated equations to solve the coupled fluid-structure system for the wind turbine would be numerically challenging and too cumbersome to apply. Recent attempts have been taken by [START_REF] Guerri | Fluid structure interaction of a parked wind turbine airfoil[END_REF] to simulate the fluid structure interaction in a strongly coupled manner. The general trend for turbomachinery is to rather use the staggered approach, i.e. using well established codes for the fluid and the structural solver separately and coupling them at a predefined step. Recently, attempts [START_REF] Dobrev | Fluidstructure interaction in the case of a wind turbine rotor[END_REF] have been taken to employ the staggered fluid-structure simulation to resolve the wind turbine structural and flow fields.

Flow Simulation

The flow simulations around the airfoil shapes are carried out with the help of the panel method code XFOIL [START_REF] Drela | An analysis and design system for low Reynolds number airfoils[END_REF]. The XFOIL code employs linear-vorticity stream function panel method coupled with viscous formulation of the flow-field via integral boundary layer method [START_REF] Drela | Viscousinviscid analysis of transonic and low Reynolds number airfoils[END_REF]. The required inputs for the XFOIL code are the airfoil shapes in plain text format and the aerodynamic boundary conditions.

Structural Simulation

The structural simulation is carried out with the help of the open source software CalculiX [START_REF] Dhondt | The Finite Element Method for Three-Dimensional Thermomechanical Applications[END_REF]. The general equation of motion for the linear static system is as follows [START_REF] Dhondt | The Finite Element Method for Three-Dimensional Thermomechanical Applications[END_REF]

[M ]{ Ü } + [K]{U } = {f } (5)
With

{U } = U (t) (6)
{U (0)} = U 0 (7) Ü (0) = Ü (0) (8)
Here [M ] is the global mass matrix, {f } is the global force vector, and [K] is the global stiffness matrix. The time derivatives U and Ü of the displacement field U represents the velocity and the displacement of the volume respectively. The global force vector {f } is time-dependent, and contains both inner and outside forces. The displacement field {U } is assumed to be a continuous function of the finite elements in the field, and is described in [START_REF] Dhondt | The Finite Element Method for Three-Dimensional Thermomechanical Applications[END_REF] as

U (X) = N i=1 φ i (η, ζ, ξ) U i (9)
The function φ i is the shape function and is defined by the mesh element type. In this case, the 10-node tetrahedral element is used as shown in figure 3.

FSI Coupling Scheme

For transient aeroelastic problems, conventional staggered approaches are taken in both serial [START_REF] Batina | Unsteady Euler Airfoil Solutions Using Unstructured Dynamic Meshes[END_REF] and parallel modes [START_REF] Weeratunga | Direct Coupled Aeroelastic Analysis Through Concurrent Implicit Time Integration on a Parallel Computer[END_REF]. Farhat and Lesoinne [START_REF] Farhat | Two Efficient Staggered Algorithms for the Serial and Parallel Solution of Three-Dimensional Nonlinear Transient Aeroelastic Problems[END_REF] have recently proposed an improved serial staggered algorithm for time dependent coupled fluid-structure model. Hence a similar strategy is applied for an assumed steady state behavior. A steady-state staggered fluidstructure interaction scheme is used to numerically simulate the shape-adaptive system. The structural simulation is performed using finite element simulation software CalculiX CrunchiX [START_REF] Dhondt | The Finite Element Method for Three-Dimensional Thermomechanical Applications[END_REF]. Figure 4 depicts the steady-state staggered FSI coupling scheme. At first, the 2D viscous flow field around the airfoil is simulated with the help of XFOIL. In the next step, structural simulation of the morphing structure under aerodynamic and actuation forces is performed using CrunchiX [START_REF] Dhondt | The Finite Element Method for Three-Dimensional Thermomechanical Applications[END_REF]. The deformed shape of the airfoil profile is then used for a second aerodynamic simulation via XFOIL. The resulting aerodynamic forces are then incorporated in further structural simulations to recalculate the required actuation forces. A parametric study is conducted to determine the shape and the position of the cam profile, the dimensions and relative positions of the beam elements, the required actuation torque, and the material characteristics of the airfoil skin to achieve optimal profile pressure distribution for the morphed airfoil shape.

FSI Coupling Constraints

The numerical scheme presented in this work is based on two dedicated softwares for the simulation of the flow field and the structural response coupled with a mapping algorithm in between the individual software. The mapping algorithm is developed with the assumption of a homogeneous aerodynamic loading along the span-wise direction of the blade element. This is valid in the case of the wind turbine application only when no rotational effect is considered and is similar in strategy to the classical blade element momentum theory [START_REF] Kulunk | Aerodynamics of Wind Turbines, Fundamental and Advanced Topics in Wind Power[END_REF]. In the actual case, the flow around the wind turbine blade would be 3D in nature and the effects could be incorporated using either 3D Navier-Stokes simulation [START_REF] Herráez | Insight into Rotational Effects on aWind Turbine Blade Using Navier-Stokes Computations[END_REF] or from experimental data [START_REF] Bak | Three-Dimensional Corrections of Airfoil Characteristics Based on Pressure Distributions[END_REF]. If the 3D effects are neglected, the aerodynamic field around the wind turbine blade can be resolved with the help of 2D fluid solver. The structural model, however, requires a 3D representation of the body to calculate the structural response for the given mass and stiffness matrix. Hence a mapping strategy is to be taken to transfer the field variables between the wetted surfaces. The mapping strategy is shown in figure 5. The 2D (z, x) fluid mesh is assumed to be symmetric in y direction. The aerodynamic pressure distribution is thus mirrored along the y axis. The aerodynamic loads are interpolated over (z, x), and transferred to the free surface of the 3D (y 1 , y 2 , y 3 ) structure mesh. A spline inter-polation scheme is used to interpolate the aerodynamic forces on the nodes of the C3D10 tetrahedral elements. A similar strategy is adopted for the coupling of the displacement matrix of the structure mesh to the fluid mesh. Hence the displacement field of the nodal elements of the structural mesh are translated to the deformed shape of the aerodynamic body. The mesh of the aerodynamic solver is then updated and the respective sweep of aerodynamic simulation is carried out. Benchmark tests of the flow solver XFOIL [START_REF] Selig | Low Reynolds Number Airfoil Design Lecture Notes[END_REF] and of the structural solver CalculiX CrunchiX [START_REF] Starnini | Calculix-Theory benchmarks: linear elastic[END_REF] are available and show reliable performance. To validate the predictions made by the individual softwares, actual experimental investigations are required on prototype blade section.

NUMERICAL SIMULATION

The boundary conditions, preparations and execution of the numerical simulations are presented in the following sections.

Geometry and Boundary Condition

For the numerical simulation the following geometric and boundary conditions are applied. Figure 6 depicts the base geometry of the shape-adaptive airfoil. The aerodynamic shape of the base airfoil is taken to be that of a NACA 0012 airfoil. The chord length of the airfoil is l c = 200 mm, the span length is l s = 198 mm. The airfoil is divided into two distinctive parts, the first 40% of the chord length is modeled as to be rigid, and the later 60% of the chord length is modeled to be shape adaptive. The shape adaptive part is divided into two slender beam sections forming the respective surface profiles of the upper and lower surfaces of NACA 0012 airfoil. The thickness of the beam elements are taken as to be t = 2 mm. To provide for the transfer of the actuator forces, two inward protruding beam sections each of 3 mm thickness are incorporated in the adaptive sections of the airfoil. The distance of the upper and lower inward-protruding beam sections from the fixed part are l upper = 22 mm and l lower = 37 mm respectively. The heights of the beams are h upper = 14 mm, and h lower = 8.7 mm respectively.

Discretization of the Aerodynamic Flow-field

The discretization of the airfoil for the aerodynamic flow simulation is directly related to the airfoil profile coordinate point distribution. Hence, numerical routine is used to generate the base airfoil and the deformed airfoil shapes. The deformation of the shape-adaptive airfoil is obtained from the CalculiX CrunchiX simulation result and is updated for subsequent aerodynamic simulation. Batch scripts are prepared to load the deformed airfoil profile in XFOIL and the geometric design routine of XFOIL is used to check and smooth the paneling of the airfoil profile. The airfoils are discretized in 226 panels.

Discretization of the Structure

The structural mesh is generated with the help of open source mesh generation software NETGEN [START_REF] Schöberl | [END_REF][START_REF] Schöberl | NETGEN -An advancing Front 2D/3D-Mesh Generator based on Abstract Rules[END_REF]. The CAD model of a blade section is loaded in NETGEN as step (ISO 10103) format. The blade-geometry contains filleted surface to eliminate structurally weak sections. In a first step, the geometric model is loaded in the NET-GEN software and analyzed for topological failure and healed if any found. Through a selection of parameters, the mesh granularity, local and global mesh-size grading, discretization of curved geometry and the number of surface and volume mesh optimization are set, and the final volume mesh is generated. The software automatically detects the edges present in the CAD model and divides the 3D mesh in several faces, with a boundary number assigned to each face. This information is used for later definition of the boundary conditions in the structural solver. The automatic mesh generator discretizes the geometry in 183451 nodal elements. A total of 57360 surface elements and 92510 volume elements are constructed from these nodal points.

FSI Sweep

In the first step, Python routine is used to generate the CalculiX GraphiX batch data to (*.fbl) to read the mesh generated by NETGEN and save the volume and surface meshes for the structural simulation. In the following step, the CalculiX GraphiX routine is run in batch mode and the element and node indices along with the boundary condition files are saved. Numerical routines using the open source software GNU Octave are developed to read the CalculiX GraphiX output file and the aerodynamic loading files obtained from the XFOIL simulation result. In the next step, the FSI mapping algorithm is implemented and the input files for the structural solver Cal-culiX CrunchiX are prepared. In the following step, GNU Octave codes for interpreting the deformed shape of the structure are run and the deformed airfoil shape is passed to the XFOIL code for aerodynamic simulation.

RESULTS

In the following sections the results of the FSI simulations are presented.

FSI Simulation First Sweep

The first FSI sweep starts with the aerodynamic flow simulation of the base airfoil with XFOIL. The resulting aerodynamic loadings are used for the structural solver. The actuator forces required for deforming the base airfoil with the target deformation angle are obtained via repeated iterations of the structural simulation.

XFOIL Simulation of Base airfoil

In a first step, the NACA 0012 airfoil is simulated for an incidence angle i = 0 • . The resulting pressure distribution around the airfoil surface profile is depicted in figure 7. For the symmetric airfoil with 0 • incidence angle, the aerodynamic loadings are identical for the upper and lower surfaces of the airfoil.

CalculiX Simulation

In the next step, the FSI coupling routine is implemented, first by reading the volume and surface meshes generated by the CalculiX GraphiX software. The aerodynamic loadings are fetched from the XFOIL simulation result of the base airfoil. Actuator forces are defined in the CalculiX CrunchiX input file, and structural simulations are carried out to determine the actuator loads for achieving the target bending of the airfoil. The structural simulations are carried out for static analysis with the help of SPOOLES solver [START_REF] Ashcraft | The User Manual for SPOOLES, Release 2.2: An Object Oriented Software Library for Solving Sparse Linear Systems of Equations[END_REF]. The resulting deflection of the airfoil profile are depicted in figure 8. The deflection D1 and D3 are along the airfoil thickness and chord directions respectively. The trailing edge of the airfoil is deflected by D1 = 28mm. The inward-protruding beam elements are rotated in the x-z plane owing to the bending of the airfoil upper and lower surfaces. These rotation is translated in a maximum shift of D3 = -4.34 mm for the upper inwardprotruding beam, and that of D3 = 3.21 mm for the lower inward-protruding beam.

Resulting Airfoil Shape

The resulting airfoil shape after the first FSI sweep is shown in figure 9. The camber line of the shape-adaptive airfoil is deformed by 8 • total, which translates to a bending of about 28 mm at the trailing edge.

FSI Simulation Second Sweep

In the second FSI sweep, the aerodynamic loading on the airfoil is updated from the aerodynamic simulation result from XFOIL. The actuator forces are kept unchanged to those of the first FSI sweep.

XFOIL Simulation of Shape-Adapted Airfoil

The aerodynamic loading from the XFOIL simulation for the shape-adapted airfoil is presented in figure 10.

The pressure distributions show rough profile which result the waviness in the airfoil panels. This is a direct result of the mapping algorithm for the FSI simulations.

CalculiX Simulation with Recalculated Aerodynamic Forces

In the next step, structural simulation are carried out with the updated aerodynamic loadings. The actuator forces are kept unchanged for this stage to investigate the effect of the changed aerodynamic loadings on the deformation of the airfoil surfaces.

Resulting Airfoil Shape

The resulting shape after the second sweep are shown in figure 11. The shapes of the airfoil profile are changed owing to the updated aerodynamic loadings. The pressure side has a total of 1.12 mm less deflection and the suction side has 0.36 mm more deflection at the trailing edge. The difference in the deflections are 0.56% and 0.18% of the chord length respectively.

CalculiX Simulation with Recalculated Actuator Forces

In the next step, structural simulations are carried out with varying actuator forces to achieve the target camber deformation of the airfoil. Hence 6 further iterations are carried out to obtain the end result. The resulting trailing edge deflection for the airfoil upper and lower surfaces are shown in figure 12. the The final actuator forces on the upper and lower beam sections are increased 4.23% and 1.51% respectively to achieve the camber line bending of 8 • .

Aerodynamic Analysis

In the following step, the aerodynamic analysis of the shape adapted airfoil is carried out. The shape-adapted airfoil has a lift coefficient of c l = 1.2867 and a drag coefficient of c d = 0.02128. The airfoil profile coefficient of pressure distributions are depicted in figure 13. Results of the first FSI sweep and those of the final FSI sweep are presented along with the coefficient of pressure distribution for the symmetric NACA 0012 airfoil. The rough distribution of the c p values are observed at sections where the deformed airfoil profile shows wavy nature. The cause of this lies on small changes in deflection values interpolated by the FSI mapping algorithm. The kinematic shape parameter over the airfoil profiles are shown in figure 14. Interestingly the values of h k for both shape adapted airfoil differ from that of the base NACA 0012 airfoil in the fixed section upstream to the shape-adaptive section. Whereas the change in along z axis. The material is stressed at the most along the z direction owing to the tensile and compressive loads from the shape adaption. The critical areas are the sections from the beginning of the shape-adaptive part to the inward-protruding beam sections where the actuator loads are applied. The upper section of the shape-adaptive part of the airfoil undergoes the highest stress. Results from initial design iterations suggest that, a shift in the position of the inward-protruding beams further towards the trailing edge would improve the situation. Hence the positioning of beams are selected at the furthest possible station for the current airfoil dimension. In the next step, the von Mises stress distribution is observed and is depicted in figure 18. The highest von Mises stress is observed in the suction side near the beginning of the shape-adaptive part of the airfoil. However, with a maximum value of 106 N/mm 2 the stresses remain well within the safety region of the selected material.

DAT1: STRESS TIME:1.000000 ENTITY:MISES max: 1.06e+03 min: 0.00e+00

Fatigue Analysis

In order to check the feasibility of the proposed design, the system should be checked for the fatigue loading. In this section, the system is checked for the stresses resulting from the shape deformation. The following assumptions and boundary conditions are considered for the prediction of the fatigue life cycle for the adaptive system. The loading is taken to be unidirectional with a stress ratio R = 0, the maximum stress is taken to be of the von Mises stress calculated by the structural solver σ max = 106 MPa. Now from empirical database presented in [START_REF]Metallic Materials and Elements for Aerospace Vehicle Structures MIL-HDBK-5J. 31[END_REF] the fatigue life cycle can be predicted for the fatigue stress with the help of the following equations σ eq = σ max (1 -R) 0.52 [START_REF] Otto | State of the art in wind turbine aerodynamics and aeroelasticity[END_REF] logN f = 20.83 -9.09log(σ eq )

It is to be noted that the equations are given for stresses in kilopound per square inch (ksi) [START_REF]Metallic Materials and Elements for Aerospace Vehicle Structures MIL-HDBK-5J. 31[END_REF]. The resulting fatigue life cycle is N f = 10 10 . Hence it is to be noted that in the actual case, the stress ratio would take a positive value and the fatigue life cycle would increase. The choice of R = 0 is here appropriate as the zero deformation load of the wind turbine blade life cycle is an unknown value at this moment.

DESIGN

The deflection D3 from the the beam elements give a pretext for the design of the actuation system. The actuator is modeled as a kinematic module able to rotate around its axis and has a cam profile which determines the actuation forces exerted on the different beam sections. The actuation force is achieved through implementing a torque which translates in a coupled bending moment on the profile upper and lower skins.

DESIGN OF THE CAM PROFILE

The design of the cam profile for the actuation is shown in figure 19. In this case a double cam profile is used to simultaneously exert the actuating forces on the upper

SUMMARY

In scope of this work, multidisciplinary design strategy of shape-adaptive airfoil for wind turbine application is presented. A low order fluid structure interaction strategy is developed and implemented for a design case scenario for actively adaptive airfoil profile. This strategy allows one to plan, design and pre-investigate potential design-flaws with very little effort. A design scenario is presented and detailed aerodynamic and structural characteristics are investigated for the prototype. A novel actuation system is proposed to implement in airfoil shape adaption.

The FSI mapping algorithm should be extended to incorporate 3D mapping. This would allow one to investigate the effect of 3D aerodynamic phenomenon on the structural response and vice versa. The initial design does not consider the unsteady nature of the fluid structure interaction. For the further design of the cam profile, the unsteady effect of the shape change has to be incorporated. For that purpose, the flow field has to be resolved for unsteady effects. Hence it calls for the inclusion of CFD simulation of the flow field around the airfoil shape and detailed experimental investigations of the aerodynamic and structural behaviors.

α

  absolute angle [rad, • ] β relative angle [rad, • ] δ boundary layer thickness [L] δ * boundary layer displacement thickness [L] ǫ modulus of elasticity [M L -1 T -2 ] η local coordinate of structure mesh [-] ρ density[M L -3 ] σ eq equivalent strength [M L -1 T -2 ] σ max maximum strength [M L -1 T -2 ] σ ts ultimate tensile strength [M L -1 T -2 ] σ y yield strength [M L -1 T -2 ] θ boundary layer momentum thickness [L] ξ local coordinate of structure mesh [-] ζ local coordinate of structure mesh [-] a specific work done by an airfoil section [L 2 T -1 ] c absolute flow velocity [LT -1 ] c d coefficient of drag [-] c l coefficient of lift [-] c p coefficient of pressure [-] h k boundary layer kinematic shape parameter [-] w relative flow velocity [LT -1 ] i incidence angle [rad, • ] l length [L] l c chord length [L] l s span length [L] p pressure [M L -1 T -2 ]r radius [L] t thickness [L] u blade speed [LT -1 ] D1 deflection along x axis [L] D2 deflection along y axis [L] D3 deflection along z axis [L] M a Mach number [-] N f fatigue life cycle [-] R stress ratio [-] ReReynolds number[-]

Figure 2 .

 2 Figure 2. Schematic of Bending of Airfoil Skin

Figure 3 .

 3 Figure 3. 10-Node Tetrahedral Element C3D10 [3]

Figure 4 .

 4 Figure 4. FSI Steady-state Staggered Coupling Scheme

Figure 5 .

 5 Figure 5. FSI Mapping Scheme

Figure 6 .

 6 Figure 6. Base Geometry of Shape-Adaptive Airfoil

Figure 7 .

 7 Figure 7. Surface Pressure Distribution for Base Airfoil

Figure 8 .

 8 Figure 8. CalculiX: Deformations (a) D1 (b) D3

Figure 9 .

 9 Figure 9. Shape-adapted Airfoil with Base Airfoil Profile

Figure 10 .

 10 Figure 10. Updated Aerodynamic Loading

Figure 11 .

 11 Figure 11. Shape-adapted Airfoils for Subsequent FSI Sweeps

Figure 12 .

 12 Figure 12. Iterations of Sweep 2

Figure 13 .

 13 Figure 13. Airfoil Profile Coefficient of Pressure Distributions

Figure 17 .

 17 Figure 17. CalculiX: Distributions of Principal Stresses (a) S XX (b)S Y Y (c)S ZZ

Figure 18 .

 18 Figure 18. CalculiX: von Mises Stress Distributions

  and lower protruding beam elements. The dimensions of the cam profile are taken to incorporate the required bending shape of the airfoil. Hence, the deflections D3 of the upper and lower beams have to be realized by actuating the cam profile. The cam is designed to have a base diameter of d b = 12 mm. The upper and lower eccentricity are taken as l 1 = 3 mm and l 2 = 2 mm respectively. These values correspond to the deflection D3 of the beams at the mid-plane. The cam has to be rotated counter-clockwise to achieve the required shape adaption of the base airfoil. A gradual slope is incorporated to hinder sudden shape change. The values of the angles are selected as γ 1 = 62 • , γ 2 = 114 • ,γ 3 = 74 • and γ 4 = 110 • . The material for the cam is taken as stainless steel (DIN X2CrNiMoCuWN25-7-4).

4 Figure 19 .

 419 Figure 19. Design of the Cam Profile

Table 1 .

 1 Boundary Conditions for Form-adaptive Wind Turbine Blade Concept.

	Goal	Subtask
		Increased startup lift
		Controllability of operational loading
	Aerodynamic	Gradual attaining of stall
		Good post-stall characteristics
		Reduced structural loading
	Structural	Increased life cycle
		Reduced weight
		Compatible with existing wind turbine
	Construction	blade structure
		Low complexity of production
	Miscellaneous	Ease of Maintenance
		Cost

ACKNOWLEDGMENTS

The authors would like to acknowledge Mr. Mark Drela for his airfoil panel method code XFOIL, Mr. Guido Dhondt and Mr. Klaus Wittig for the structural solver CalculiX, and Mr. Jochaim Schöberl for his mesh generation software NETGEN.

the pressure side occur close to the shape-adaptive section, the suction side distributions of h k exhibit an early peak near the nose section. For the shape-adaptive airfoils from both FSI sweeps, distinctive jump in h k is observed upstream of the shape-adaptive section. The values of h k reach 9 for both the cases. An investigation of the boundary layer shows, however, no trace of flow separation. The boundary layer displacement thickness over the airfoil profiles are shown in figure 15. Also in this case, a shift in the change in boundary layer upstream of the shape-adaptive sections are observed for the pressure sides of the shape-adaptive airfoil. For the suction side, the values of δ * increase starting from the shapeadaptive section. The boundary layer momentum thickness over the airfoil profiles are shown in figure 16. For the shape-adaptive airfoils, the suction side boundary layer momentum thickness θ increases way upstream of the shape-adaptive section. This, in combination of the relatively unchanged boundary layer displacement thickness distribution (figure 15) results in the difference in the suction side kinematic shape factor change relative to the NACA 0012 airfoil observed in (figure 14). Also noticeable is the change of the boundary layer momentum thickness distribution for the pressure sides upstream of the shape-adaptive sections. Comparing the onset of these changes to those of the boundary layer displacement thickness (figure 15), the shift of the kinematic shape factors in (figure 14) for the pressure sides are explained.

Structural Analysis

The resulting shape-adapted airfoil is then analyzed for structural properties. The material of the airfoil is taken as Aluminium alloy (EN AW 2024-T4 [AlCu4Mg1]). The alloy has a density of ρ = 2.77 g/cm 3 an ultimate tensile strength σ ts = 469 MPa, yield strength of σ y = 324 MPa, and a modulus of elasticity of ǫ = 73.1 GPa [START_REF] Hesse | Aluminium-Schlüssel:= Key to aluminium alloys 7[END_REF].

In a first step, stress distributions are investigated to check for design weak points. The principal stresses of the shape-adapted airfoil are presented in figure 17. The principal stress along x axis is presented in figure 17 (a). The bases of the upper and lower inward-protruding beams undergo the maximum stresses in this direction. The values remain within the safety limit for the given material. Figure 17 (b) depicts the principal stress along y axis. Although the value remains within the safety limit, the presence of stresses along y direction suggests compound bending of the airfoil surfaces along the spanwise direction. Figure 17 (c) depicts the principal stress