Blade Element Method

The blade element method formulation is based on the methodology that is proposed by Rajagoplan [279] and others. It is primarily designed to model the effects of a helicopter rotor on the flow field without actually prescribing the exact geometry of rotor blades.

This approach consists of two parts. The first step is to identify the location of cells within the computational domain that belong to the virtual disk and to which the source term is applied. The second step involves the calculation of the momentum source term vector for each cell, which was identified and marked in the first step.

Step 1: Cell Marking and the Allocation to the Blade Elements/Buckets Based on the Interpolation Grid

The information for this step is user-specified and provided through the virtual disk nodes:

Disk Geometry
Describes the geometrical dimensions and the orientation of the disk
Virtual Disk Resolution
Specifies the size of the two-dimensional interpolation grid in polar coordinates, representing the actual blade elements. The blade properties such as chord, twist, and aerodynamic characteristics prevail throughout the whole element.

Source Term Distribution: Time Averaged

For the Time Averaged approch, the source term is calculated per blade element on the interpolation grid and then transferred to the finite volume mesh. Each element is associated with a number of cells on the finite volume mesh. For the coupling between the local flow conditions and the modeled rotor, the information is taken from the finite volume cell which is the closest to the center of a blade element on the interpolation grid. Thus, the specified size of the interpolation grid is the actual computational effort of the influence of the rotor on the flow field. The sketch below explains the flow field coupling between the finite volume mesh (red) and the virtual disk interpolation grid (black).

The marking procedure is itself a two-step procedure. First, all the cells that fall within the circular disk area are marked to receive the source term. Second, each of the marked cells is assigned a bucket/element number corresponding to its location with respect to the interpolation grid. So each of the cells that was marked during the first step belongs to one of the elements/buckets of the interpolation grid.

This approach requires that each element/bucket has at least one cell in its domain. If the finite volume mesh is very coarse and the interpolation grid is too fine then, you are asked to reduce the resolution. The following two figures show the interpolation grid for two different resolutions. For the fine interpolation grid in the second figure, no finite volume cell come under the domain of some for some of the elements/buckets at the innermost radius.

Source Term Distribution: Time Accurate

The Time Accurate approach tracks the motion of the blades and adds source terms only to the volume cells which correspond to the location of the blades. For a rotor with four blades, the cells are marked as shown below according to the blade shape and location specified within the Disk Geometry node.

To ensure successful coupling with the interpolation grid, the finite volume mesh must be sufficiently fine that at least 5 cells are marked across the blade width. A similar mesh resolution should be applied for the entire sweep area of the rotor. In the example above, 10 cells are marked.

With the default interpolation mode, Simcenter STAR-CCM+ supports overlapping of the marked cells of different blades. The source term is computed accordingly:

Step 2: Calculation of Momentum Source Term Strength

The second part of the blade element method calculates the strength of momentum source term for the elements on the interpolation grid. Based on the geometry of the rotor (that is orientation and origin), a disk local cartesian coordinate system $(X′, Y′, Z′)$ is created. This cartesian coordinate system has its origin at the center of the rotor and the $Z′$ -axis in the direction corresponding to the thrusting direction of the disk. The $Z′$ -axis is perpendicular to the plane of rotation, while the $X′$ and $Y′$ axes lie in the plane of the disk. The following sketch demonstrates the disk local coordinate system:

The rotational speed $ω$ is always oriented along the $Z′$ -direction. + $ω$ corresponds to a counter-clockwise wise (CCW) rotation and - $ω$ corresponds to a clock-wise (CW) rotation when viewed from the + $Z′$ direction.

To account for blade flapping, two additional coordinate systems are defined:

A rotor-based cylindrical system $(r, ϕ, z)$ with respect to the basis of the local coordinate system $(X′, Y′, Z′)$
A Cartesian coordinate system with respect to the basis of the cylindrical coordinate system. It is called the flap coordinate system.

The basis of the flap system $(s, ϕ, n)$ is defined such that $e_{s}$ points in the spanwise direction of the blade. The basis vector $e_{ϕ}$ of the flap coordinate system points in the same direction as $e_{ϕ}$ of the rotor-based cylindrical system. The cylindrical system rotates with the blade and the flap system is tilted with respect to the cylindrical system. Given the azimuth location $ψ$ of the blade, the cylindrical system is rotated such that $e_{r}$ is aligned with the blade span. The flap system is tilted by $δ (ψ)$ , that is, by the flap angle of this azimuth location. All calculations are performed in the flap system. The resulting force vector is transformed back to the lab system and then added to the momentum equation as a source term.

The local flow field velocity in the laboratory coordinate system is:

图 1. EQUATION_DISPLAY

v = v_{x} \hat{i} + v_{y} \hat{j} + v_{z} \hat{k}

(4990)

The local flow velocity is transferred to the flap coordinate system using a sequence of transformation matrices $M_{1}, M_{2}, M_{3}$ :

图 2. EQUATION_DISPLAY

[\begin{matrix} v_{s} \\ v_{ϕ} \\ v_{n} \end{matrix}] = M_{3} M_{2} M_{1} [\begin{matrix} v_{x} \\ v_{y} \\ v_{z} \end{matrix}]

(4991)

where the transformation matrices relate the basis vectors of the various coordinate systems as:

图 3. EQUATION_DISPLAY

[\begin{matrix} e_{s} \\ e_{ϕ} \\ e_{n} \end{matrix}] = M_{3} [\begin{matrix} e_{r} \\ e_{ϕ} \\ e_{z} \end{matrix}] [\begin{matrix} e_{r} \\ e_{ϕ} \\ e_{z} \end{matrix}] = M_{2} [\begin{matrix} e_{x^{'}} \\ e_{y^{'}} \\ e_{z^{'}} \end{matrix}] [\begin{matrix} e_{x^{'}} \\ e_{y^{'}} \\ e_{z^{'}} \end{matrix}] = M_{1} [\begin{matrix} e_{x} \\ e_{y} \\ e_{z} \end{matrix}]

(4992)

The blade has a velocity due to rotation, which can be written in the flap coordinate system as:

$v^{b l} = M_{3} M_{2} M_{1} (ω \times r)$

where $r$ is the radius vector.

The velocity relative to the blade $v^{r e l} = v_{s} e_{s} + v_{ϕ} e_{ϕ} + v_{n} e_{n}$ is thus given by:

$\begin{matrix} v^{r e l} = v - v^{b l} \\ = M_{3} M_{2} M_{1} v - M_{3} M_{2} M_{1} (ω \times r) \end{matrix}$

To calculate the effective angle of attack

α

, the following blade angles are taken into account:

induced angle of attack (inflow angle) $β$
blade flap angle $δ (ψ)$
geometrical pitch angle $θ (ψ)$
geometrical twist angle $θ_{t w} (r)$

The blade flap angle is calculated as:

图 4. EQUATION_DISPLAY

δ = δ_{0} - δ_{1 c} \cos (ψ) - δ_{1 s} \sin (ψ)

(4993)

where

δ_{0}

is the coning angle,

δ_{1 c}

is the longitudinal flap angle, and

δ_{1 s}

is the lateral flap angle. The time derivative is given by:

图 5. EQUATION_DISPLAY

\dot{δ} = ω \frac{ⅆ δ}{ⅆ ψ} = ω (δ_{1 c} \sin (ψ) - δ_{1 s} \cos (ψ))

(4994)

and the flap velocity:

图 6. EQUATION_DISPLAY

v^{f l a p} = v_{s}^{f l a p} e_{s} + v_{ϕ}^{f l a p} e_{ϕ} + v_{n}^{f l a p} e_{n} = r \dot{δ}

(4995)

The induced angle of attack is now calculated as:

图 7. EQUATION_DISPLAY

β = \arctan [\frac{- (v_{n} + v_{n}^{f l a p})}{s i g n (v_{ϕ}) \sqrt{{(v_{ϕ})}^{2} + {(v_{s})}^{2}}}]

(4996)

where the direction of $v_{ϕ}$ in the reverse flow region for high advance ratios is taken into account by the $s i g n (v_{ϕ})$ .

It can be seen from Eqn. (4996) that the blade flap angle affects the inflow angle (induced angle of attack) through the blade flap velocity, and thus the effective angle of attack $α$ .

The pitch angle $θ$ is defined as:

图 8. EQUATION_DISPLAY

θ = θ_{0} - θ_{1 c} \cos ψ - θ_{1 s} \sin ψ

(4997)

where $ψ$ is the azimuth angle, $θ_{0}$ is the collective pitch angle, $θ_{1 c}$ is the lateral cyclic pitch angle, and $θ_{1 s}$ is the longitudinal cyclic pitch angle.

The final effective angle of attack is then calculated as:

图 9. EQUATION_DISPLAY

α = θ + θ_{t w} - β

(4998)

Knowing the angle of attack $α$ and the blade relative velocity $v^{r e l}$ experienced by the airfoil, either the Mach number or the Reynolds number are calculated. The aerodynamic force coefficients are also used to calculate the Mach and the Reynolds number. You specify these coefficients in the form of tables in the simulation tree. Which number is calculated depends on how you specify the aerodynamic force coefficients—as a function of Mach number $C_{l}$ , $C_{d} = f (α, M a) or$ or as a function of Reynolds number $f (α, R e)$ .

The lift and drag forces are then computed as:

图 10. EQUATION_DISPLAY

\begin{matrix} L′ = \frac{1}{2} \cdot ρ \cdot {| v^{r e l} |}^{2} \cdot C_{l} \cdot c \cdot d r \\ D′ = \frac{1}{2} \cdot ρ \cdot {| v^{r e l} |}^{2} \cdot C_{d} \cdot c \cdot d r \end{matrix}

(4999)

where $c$ is the blade chord.

To account for tip-loss due to trailing vortices at the tip of a rotor blade, the lift force $\begin{matrix} L′ \end{matrix}$ , the angle of attack $α$ or lift and drag $L' a n d D'$ are corrected depending on the chosen tip correction method.

The forces are then converted into the normal and tangential force (transformation from the wind axis to the rotor flap system):

$\begin{matrix} {F′}_{n} = L′ \cos β - D′ \sin β \\ {F′}_{t} = L′ \sin β + D′ \cos β \end{matrix}$

As there are no aerodynamic forces along the blade span, $F_{s} = 0$ . The resultant aerodynamic force vector can be assembled on the blade element as $(F_{s}, F_{ϕ}, F_{n})$ in the flap system.

Finally, the force vector is transformed into the laboratory coordinate system:

图 11. EQUATION_DISPLAY

(\begin{matrix} F_{x} \\ F_{y} \\ F_{z} \end{matrix}) = M_{1}^{T} M_{2}^{T} M_{3}^{T} (\begin{matrix} F_{s} \\ F_{ϕ} \\ F_{n} \end{matrix})

(5000)

The instantaneous force acting on the fluid element is then $- F$ .

The time-averaged source terms $S = (S_{x}, S_{y}, S_{z})$ are added to the discretized momentum conservation equations. Since the blade only spends a finite fraction of its total revolution time passing through the control volume, the time-averaged source terms are scaled by the time fraction:

图 12. EQUATION_DISPLAY

S = \frac{b Δ φ}{2 π} (- F)

(5001)

where $b$ is the number of blades and $Δ φ$ is the angular distance through which the blade traverses in passing through a blade element.

The time-accurate source terms do not require any scaling. The specific source term from a given blade element is distributed over the respective cells in the finite volume mesh, proportional to their cell volume.

Numerical Trim Algorithm

In the context of rotor crafts or helicopters, the term trim means the correct adjustment of the aircraft controls and the overall orientation of the rotorcraft in order to obtain a desired steady flight condition. For example, for the hover condition, the objective of trimming the rotor is to achieve an equilibrium, that is, that the net sum of forces and the net sum of moments about the center of gravity are zero. Simcenter STAR-CCM+ uses the isolated trim (wind tunnel trim) algorithm. The isolated trim algorithm evaluates the rotor pitch controls, which are the collective, the lateral cyclic, and the longitudinal cyclic pitch.

Main Rotor Controls

Generally, a helicopter has four independent controls: three for the main rotor and one for the tail rotor. In the Simcenter STAR-CCM+ implementation, the rotor is treated in isolation and only the rotor controls are considered.

The main rotor controls are:

Collective pitch:
By using this control, a pilot can change the pitch angle of all the blades of the main rotor collectively, that is, by the same amount. This change effectively influences the total thrust generated by the main rotor.
Cyclic pitch:
By using this control, a pilot can change the pitch angle once per revolution. The cyclic pitch is categorized as follows:
- Lateral cyclic pitch $θ_{1 c}$ :
  Varies as per $\cos (ψ)$ . The lateral cyclic pitch is at a maximum or minimum when a blade is at 0° and 180°. The rotor response, which is phased by 90°, leads to a control effect in the lateral direction (y-axis) and effectively influences the roll moment.
- Longitudinal cyclic pitch $θ_{1 s}$
  Varies as per $\sin (ψ)$ . The longitudinal cyclic pitch is at a maximum or minimum when a blade is at 90° and 270°. The rotor response, which is phased by 90°, leads to a control effect in the longitudinal direction (x-axis) and effectively influences the pitch moment.

The total pitch angle therefore is a function of the azimuthal angle

ψ

and is given by:

图 13. EQUATION_DISPLAY

θ = θ_{0} - θ_{1 c} \cos (ψ) - θ_{1 s} \sin (ψ)

(5002)

Algorithm

The rotor aerodynamic response is a highly coupled function of its control inputs and its aerodynamic environment. Although the thrust generated by rotor is mainly influenced by collective pitch, it is also influenced by the cyclic pitch inputs. Likewise, the roll and pitch moments are directly linked to cyclic pitch inputs but also produce a change in thrust. Defined in the local coordinate system of the virtual disk:

图 14. EQUATION_DISPLAY

\begin{matrix} T & = & T (θ_{0}, θ_{1 c}, θ_{1 s}) \\ M_{x} & = & M_{x} (θ_{0}, θ_{1 c}, θ_{1 s}) \\ M_{y} & = & M_{y} (θ_{0}, θ_{1 c}, θ_{1 s}) \end{matrix}

(5003)

where $T$ is the thrust, $M_{x}$ is the roll moment, and $M_{y}$ is the pitch moment.

This coupled system of equations can be interpreted as how the control input vector $x$ is related to the rotor response vector $y$ where:

x = (\begin{matrix} θ_{0} \\ θ_{1 c} \\ θ_{1 s} \end{matrix}) and y = (\begin{matrix} T \\ M_{x} \\ M_{y} \end{matrix})

Each flight condition corresponds to a state of the rotor response vector $y$ , which requires a particular state of the control input vector $x$ , which in turn affects the aerodynamic characteristics of a helicopter.

To investigate the aerodynamic performance of a helicopter, you first choose a flight condition. This means that the response vector $y$ is known, but not the control input vector $x$ . You start the iterative procedure with an initial guess of the input vector $x$ . $x$ is then iteratively updated until the solution converges towards the target response vector $y$ .

For the iterative procedure to update

x

, the coupled system of equations is linearized. Then the Newton-Raphson approach is applied. To linearize, the Taylor series expansion for

y

about

x

is applied:

图 15. EQUATION_DISPLAY

y (x + Δ x) = y (x) + J Δ x + {O (Δ x)}^{2}

(5004)

Considering only the first-order terms, Eqn. (5004) is rearranged to read:

图 16. EQUATION_DISPLAY

Δ x = J^{- 1} {y (x + Δ x) - y (x)}

(5005)

where $J$ is the Jacobian matrix of the dependent quantities, that is, of the response vector $y$ with respect to the control inputs $x$ .

Eqn. (5005) is rewritten to:

图 17. EQUATION_DISPLAY

{(\begin{matrix} Δ θ_{0} \\ Δ θ_{1 c} \\ Δ θ_{1 s} \end{matrix})}_{(i)} = {(\begin{matrix} \frac{\partial T}{\partial θ_{0}} & \frac{\partial T}{\partial θ_{1 c}} & \frac{\partial T}{\partial θ_{1 s}} \\ \frac{\partial M_{x}}{\partial θ_{0}} & \frac{\partial M_{x}}{\partial θ_{1 c}} & \frac{\partial M_{x}}{\partial θ_{1 s}} \\ \frac{\partial M_{y}}{\partial θ_{0}} & \frac{\partial M_{y}}{\partial θ_{1 c}} & \frac{\partial M_{y}}{\partial θ_{1 s}} \end{matrix})}_{(i)}^{- 1} \times {{(\begin{matrix} T \\ M_{x} \\ M_{y} \end{matrix})}_{(d)} - {(\begin{matrix} T \\ M_{x} \\ M_{y} \end{matrix})}_{(i)}}

(5006)

In Eqn. (5006), the subscript (i) corresponds to the i^th trim iteration and d corresponds to the desired target state of the rotor response vector $y$ .

The initial guess for $x$ , yields an initial $y$ , which is usually not equal to the target state of $y_{d}$ . To proceed with the trim iteration, the Jacobian matrix must be determined using a frozen flow field. To achieve this, each component of the control input vector is independently perturbed and the rotor response is recomputed. The partial derivatives are then computed using 2^nd-order central differences. With $J$ available, the correction to the control input vector $Δ x$ is determined by solving Eqn. Eqn. (5005). The computed correction is under-relaxed to achieve numerical damping and is then used to update the control input vector. The procedure is repeated until the solution converges towards the target $y$ within a specified tolerance.

In time-average mode, the control input correction is applied repeatedly. The convergence of the source term triggers the application of the correction.

In time-accurate mode, the rotor runs at initial condition for a few revolutions until flow is periodic. Before applying the control input correction, the trim solver collects the data of $y$ for one revolution to build average values. The event of application is steered by either the convergence of the source term or by the number of revolutions between two corrections, which are specified in the Trim Option properties.

Overview of the Source Term Update

The following flowchart provides an overview of the blade element source term update process. This procedure is triggered when the main flow solver is ready to accept source term contributions from optional models into the momentum equation. When the main flow related fields such as pressure and velocity are available, the blade element source term is evaluated.

Next, Simcenter STAR-CCM+ checks if the blade element source term is converged up to the Source Tolerance (see User Guide > Modeling Physics > Virtual Disk Model > Blade Element Method > Setting the Source Convergence Control). Furthermore, Simcenter STAR-CCM+ checks if the ramping of the rotation rate is finished. If either of the checks fails, the blade element source term is under-relaxed by the Source Under Relaxation Factor (see see User Guide > Modeling Physics > Virtual Disk Model > Blade Element Method > Setting the Source Convergence Control) and is added to the source term of the momentum equation. The solution process then continues with the next model. This step corresponds to the no branch of the first decision block in the flow chart above.

When the first decision block evaluates to yes, Simcenter STAR-CCM+ checks if the main rotor response equals the target values. This check corresponds to the second decision block in the flow chart above and is interpreted as the convergence of the trim state. The trim state is converged if the user-specified Trim Tolerance is reached (see User Guide > Modeling Physics > Virtual Disk Model > Blade Element Method > Setting the Trim Convergence Control). If the trim state is converged, the solution process continues with the next model. If the trim state is not converged, the rotor trim angles are updated according to the numerical trim algorithm. The computed trim angles are under-relaxed by the Trim Under Relaxation Factor (see User Guide > Modeling Physics > Virtual Disk Model > Blade Element Method > Setting the Trim Convergence Control).

The source term computation is performed at each iteration. Therefore, any changes in the user setup resulting in a different flow field, influence the blade element method source term vector. It is important to note that the trim algorithm is executed only after the blade element source term has converged to the specified level. If the external flow field is of highly unsteady nature, you are recommended to relax the convergence tolerance on the source term so that the trim angle update can spring into action. For extreme case, where the main rotor is operating in a wake and the desired flight condition is hover, you are advised to deactivate the source term convergence check by specifying a tolerance level of 100%. This setting ensures that the algorithm is executed at each iteration.