High Level Overview

At a high level, this method uses a dynamic inflow model that is capable of quickly computing the induced velocity of a rotor, to advect vortex elements downstream. These released vortex elements are then used to compute the inflow at the rotor. This allows us to quickly compute a dynamically reacting wake without the cost of computing the influence of each vortex element with all other vortex elements. The dynamic inflow model is based on potential flow theory, which allows the use of linear superposition to create a combined multirotor flow field.

Solution Procedure

The present model is inspired by Van der Wall’s [vdW12] generalization of Beddoes’ inflow prescribed wake model [Bed85]. A high-level flow chart outlining the solution procedure is shown in figure Fig. 1. The convection of the wake is loosely coupled to a nonuniformly loaded actuator disk dynamic inflow model. The coupled models are stepped until a quasi-steady state is achieved. The simulation is initialized with no trailed wake and uses a simple harmonic solution to equations (8) and (9) to initialize the dynamic inflow model states. At the start of each timestep the inflow on the blades is computed from the Biot-Savart integration of the wake (this is initially zero) and blade loads are computed. With the blade loads known the circulation model can then be updated. The wake model then releases new vortex elements at the blade tips. The dynamic inflow model will then be used to compute the induced velocity at each vortex element. The blade loads are then integrated into the pressure coefficients needed by the dynamic inflow model. Finally equations (8), (9), and (48) are updated to the next timestep.

Fig. 1 High level flow chart of the model

Dynamic Inflow Model

The dynamic inflow model used for this work is an implementation of the model developed in the following references Huang et al. [HPP15] with extensions provided by Gunar et al. [GPP21]. It provides the inflow caused by an actuator disk anywhere in space through a blending of several different inflow solutions depending on where they are best suited. As part of the Peters dynamic inflow model lineage, this model solves the incompressible potential flow equations along the free-steam line. By the nature of the potential flow equations being a form of Laplace’s equation, the system can be solved in ellipsoidal coordinates for the pressure discontinuity on the disk using a separation of variables solution with the pressure potential being one of either:

(5)\[ \begin{align}\begin{aligned}\Phi^m_n &= \frac{1}{\nu}P^m_n(\nu)Q^m_n(i\eta)\cos(m\psi)\\\Phi^m_n &= P^m_n(\nu)Q^m_n(i\eta)\cos(m\psi)\end{aligned}\end{align} \]

The first form is used with Nowak-He (NH) state variables and the second with Morillo-Duffy (MD) state variables. The NH form converges well on the disk while the MD form converges well everywhere else.

The derivation of the model will largely be skipped here for brevity, but an overview of the implementation details follows. There are two base velocity formulations and their adjoints that the rest of the solution builds upon, based on the two different potential forms:

\begin{align} V_{MD} &= \sum_{m,n}^\infty a^m_n P^m_n(\nu)Q^m_n(i\eta)\cos(m\psi) \\ V_{NH} &= \sum_{m+n=odd}^\infty \alpha^m_n \frac{1}{\nu}P^m_n(\nu)Q^m_n(i\eta)\cos(m\psi) \nonumber \\ &+ \sum_{m+n=even}^\infty a^m_n P^m_n(\nu)Q^m_n(i\eta)\cos(m\psi) \\ V_{MD}^* &= \sum_{m,n}^\infty \Delta^m_n P^m_n(\nu)Q^m_n(i\eta)\cos(m\psi) \\ V_{NH}^* &= \sum_{m+n=odd}^\infty \Lambda^m_n \frac{1}{\nu}P^m_n(\nu)Q^m_n(i\eta)\cos(m\psi) \nonumber \\ &+ \sum_{m+n=even}^\infty \Delta^m_n P^m_n(\nu)Q^m_n(i\eta)\cos(m\psi) \end{align}

where \(a^m_n\) are the MD state variables and \(\alpha^m_n\) are the NH state variables and \(\Delta^m_n\) and \(\Lambda^m_n\) are their respective adjoint states. These state variables are transformations of one-another through the following relation:

\begin{align} &\{a^m_j\} = [A^m_{nj}]^{-1}\{\alpha^m_n\} \\ &\{\Delta^m_j\} = [A^m_{nj}]^{-1}\{\Lambda^m_n\} \\ &[A^m_{nj}] = \frac{(-1)^\frac{n + j - 2m}{2}2\sqrt{2n + 1}\sqrt{2j + 1}}{\sqrt{H^m_n}\sqrt{H^m_j}(n + j)(n + j + 2)\left[(n + j)^2 - 1 \right]} \\ &H^m_n = \frac{(n + m - 1)!!(n - m - 1)!!}{(n + m)!!(n - m)!!} \end{align}

The state variables are computed through the following system:

(8)\[\boldsymbol{M}_{NH}\{\dot\alpha^m_n\} + \boldsymbol{D}\boldsymbol{V}\boldsymbol{L}^{-1}\boldsymbol{M}_{NH}\{\alpha^m_n\} = \boldsymbol{D}\{\tau^m_n\}\]

(9)\[\begin{split}-\boldsymbol{M}_{NH}\{\dot\Lambda^m_n\} + \boldsymbol{D}\boldsymbol{V}\boldsymbol{L}^{-1}\boldsymbol{M}_{NH}\{\Lambda^m_n\} = \boldsymbol{D} \begin{bmatrix}\ddots & & \\ & (-1)^{j + 1} & \\ & & \ddots\end{bmatrix}\{\tau^m_n\}\end{split}\]

Where \(\boldsymbol{V}\) is the diagonal mass flow matrix:

(10)\[\begin{split}\boldsymbol{V} = \begin{bmatrix} V_T & & & \\ & V & &\\ && V & \\ &&& \ddots \end{bmatrix}\end{split}\]

Where \(V_T\) is the total average flow at the rotor and \(V\) is the mass-flow parameter:

(11)\[ \begin{align}\begin{aligned}&V_T = \sqrt{V^2_\infty\sin^2{\chi} + (V_\infty\cos{\chi} + \bar v_z)^2}\\&V = \frac{V^2_\infty\sin^2{\chi} + (V_\infty\cos{\chi} + \bar v_z)(V_\infty\cos{\chi} + 2\bar v_z)}{V_T}\end{aligned}\end{align} \]

\(\bar v_z\) is the average induced velocity at the rotor:

(12)\[\bar v_z = \sqrt{3}\begin{bmatrix} 1 & 0 & \hdots \end{bmatrix}\boldsymbol{L}^{-1}\boldsymbol{M}_{NH}\{\alpha\}\]

The equations for the rest of the matrices can be found in reference [MP02].

With the system in place to time-march the required states, the final induced velocity computation can take place. The final induced velocity contains several different blends to ensure an accurate solution can be found both above and below the rotor as well as upstream and down stream. The first blended solution, called \(V_{BL}\), blends together the NH solution on the disk, and the MD off the disk, i.e.,

(13)\[ \begin{align}\begin{aligned}V_{BL} = \frac{1}{1 + bh}V_{NH} + \frac{bh}{1 + bh}V_{MD}\\V_{BL}^* = \frac{1}{1 + bh}V_{NH}^* + \frac{bh}{1 + bh}V_{MD}^*\end{aligned}\end{align} \]

where

(14)\[\begin{split}h = \begin{cases} 0 & \text{if } \eta < \varepsilon \\ \eta - \varepsilon & \text{if } \eta \geq \varepsilon \end{cases}\end{split}\]

and

(15)\[\begin{split}b = \begin{cases} 20\left[1 - \frac{y^2\sin{\chi}}{1 + \eta^2}\right] & x \leq 0, |y| \leq 1 \\ & \\ 20\left[1 - \frac{y^2\sin{\chi}}{1 + \eta^2 + 0.615(y^2 - 1)}\right] & x \leq 0, |y| > 1 \\ & \\ 20\left[1 - \frac{(x^2 + y^2)\sin{\chi}}{1 + \eta^2}\right] & x > 0, |y| \leq 1 \\ &\\ 20\left[1 - \frac{(x^2 + y^2)\sin{\chi}}{1 + \eta^2 + 0.615(y^2 - 1)}\right] & x > 0, |y| > 1 \end{cases}\end{split}\]

A value of \(\varepsilon = 0.01\) was used in this work.

Next, the downstream solution is constructed to account for the fact that the trailed rotor wake does not immediately decay to zero in edgewise flight. Leveraging the adjoint theorem and the blended solution the downstream velocity and its adjoint are given by

(16)\[\begin{split}V_{DS}(x, y, z, y) &= V_{BL}(-s_0, y, z, t - \sigma\sin{\chi}) \\ & + V_{BL}^*(+s_0, -y, z, t - \sigma\sin{\chi}) \\ & - V_{BL}^*(\sigma+s_0, -y, z, t)\end{split}\]

(17)\[\begin{split}V_{DS}^*(x, y, z, y) &= V_{BL}^*(-s_0, y, z, t + \sigma\sin{\chi}) \\ & + V_{BL}(+s_0, -y, z, t + \sigma\sin{\chi}) \\ & - V_{BL}(\sigma+s_0, -y, z, t)\end{split}\]

where \(s_0\) is the x-distance downstream where the blended solution is converged:

(18)\[\begin{split}s_0 = \begin{cases} \sqrt{1 - y^2 - z^2} & \text{if } y^2 + z^2 < 1 \\ 0 & \text{if } y^2 + z^2 \geq 1 \end{cases}\end{split}\]

and \(\sigma\) is the distance downstream to the desired point from \(s_0\), such that \(\sigma = -x - s_0\). The final velocity and its adjoint are then given by

(19)\[ \begin{align}\begin{aligned}V_F = V_{BL}[1 - f(s_0, \chi, y)] + V_{DS}[f(s_0, \chi, y)]\\V_F^* = V_{BL}^*[1 - f(s_0, \chi, y)] + V_{DS}^*[f(s_0, \chi, y)]\end{aligned}\end{align} \]

where

(20)\[\begin{split}f(\sigma, \chi, y) = \begin{cases} \frac{\sin^2{\chi}}{\sin^2{\chi} + \sigma g(\chi)} & \text{if } |y| \leq 1 \\ \frac{\sin^2{\chi}}{\sin^2{\chi} + (\sigma + 1.5\sqrt{y^2 - 1})g(\chi)} & \text{if } |y| > 1 \end{cases}\end{split}\]

and

(21)\[g(\chi) = 1.84\cos^{1/2}{\chi} - 4.06\cos{\chi} + 11.84\cos^{3/2}{\chi}\]

For the flow above the rotor disk, equation (19) is all that is needed. Below the disk one last application of the adjoint theorem is required:

\begin{align} V(x, y, z, t) &= V(z, r, \psi, t) \nonumber \\ & = V_F(t - s, 0, r, \psi) \nonumber \\ & + V_F^*(t - s, 0, r, \psi + \pi) \nonumber \\ & - V_F^*(t, -z, r, \psi + \pi) \end{align}

where \(s\) is the distance from the desired point to the point on the rotor along the free-stream line.

As seen in equations (16), (17), and (22), data from previous times is required to compute the velocity both below and downstream of the rotor disk. Further, equation (9) is unstable in forward integration due to the nagative sign. To address these issues an extension developed by Guner et al. [GPP21]. In this extension employs a quasi-steady approximation to compute the adjoint states at the current time, as well as a scheme to avoid reverse integration to compute past data. This is done by dropping the derivative terms of equations (8) and (9). Then solving equation (9) for the adjoint states yields:

(23)\[\begin{split}\{\Lambda^m_n\} = \boldsymbol{M}_{NH}^{-1}\boldsymbol{L}\boldsymbol{V}^{-1} \begin{bmatrix}\ddots & & \\ & (-1)^{j + 1} & \\ & & \ddots\end{bmatrix}\{\tau^m_n\}\end{split}\]

The solving the quasi-steady equation (8) for \(\{\tau^m_n\}\) and plugging into the previous equation yields:

(24)\[\begin{split}\{\Lambda^m_n\} = \boldsymbol{M}_{NH}^{-1}\boldsymbol{L}\boldsymbol{V}^{-1} \begin{bmatrix}\ddots & & \\ & (-1)^{j + 1} & \\ & & \ddots\end{bmatrix}\boldsymbol{V}\boldsymbol{L}^{-1}\boldsymbol{M}_{NH}\{\tau^m_n\}\end{split}\]

To avoid reverse integration past states are stored and an interpolation is performed to the desired time. Then by equation (24), the past adjoint states are also computed.

There are also instances where data from a future time are required. These instances arise when computing the downstream velocity very close to the disk. This largely avoided due to the downstream blending function, but there is a small region where the blending function hasn’t fully removed the downstream solution. For this another quasi-steady approximation is used. When data for a future time are required, they are instead approximated using the states from one full rotor revolution in the past.

Lift Model

Blade loads are computed using blade element theory. The lift at each blade element is computed as

(25)\[\mathrm{d}L = \frac{1}{2}\rho ||{\boldsymbol U}||^2 c C_l \mathrm{d}y\]

The local blade element flow velocity \(\boldsymbol U\) is defined as

(26)\[\boldsymbol U = u_t \hat x + u_p\hat z\]

The component \(u_t\) is the local flow velocity tangent to the rotor tip-path plane (TPP) and normal to the blade element. The component \(u_p\) is the flow velocity perpendicular to the TPP. The tangential velocity component is computed as

(27)\[u_t = \Omega y + V_\infty\cos(\alpha_s)\cos(\psi)\]

The perpendicular velocity is computed through the Biot-Savart integration of both the shed and trailed wake. This is discussed in detail below. The non-circulatory unsteady portion of the lift is taken care of through the shed-wake vortex elements allowing the use of the steady lift equation here.

The inflow angle into the blade element is computed as

(28)\[\phi = \tan^{-1}\left(\frac{u_p}{u_t}\right)\]

which then allows the computation of the local angle of attack

(29)\[\alpha = \theta - \phi\]

where \(\theta\) is the local blade twist angle. This angle of attack is then fed into the circulation model to compute the blade bound vortex strength.

Circulation Model

To compute the bound circulation on the blade the Weissinger-L model [Wei47] is used with extensions provided by Wickenheiser and Garcia [WG07]. The Weissinger-L model treats a tapered and swept blade as a lifting line with a known angle of attack at each blade section. The extensions provided by Wickenheiser and Garcia allow for arbitrarily curved blade planforms. The blade is divided into \(N\) colocation points with a cosine spacing such that:

(30)\[ \begin{align}\begin{aligned}\phi_\nu &= \frac{\nu\pi}{N + 1}\\r &= \cos(\phi_\nu)\end{aligned}\end{align} \]

where \(\nu = 0,...,N\) and :math:` -1 < r < 1`. The local angle of attack at the blade is then related to the normalized circulation strength \(G\) and blade geometry through the following integral:

(31)\[\alpha(\phi_\nu) = \frac{1}{2\pi}\int_0^\pi \frac{G'(\phi)}{\cos(\phi) - \cos(\phi_\nu)}\mathrm{d}\phi - \frac{1}{4\pi}\int_0^\pi P(\phi_\nu, \phi) G'(\phi)\mathrm{d}\phi + \frac{R}{4\pi c(\phi_\nu)}\int_0^\pi P(\phi_\nu, \phi) G'(\phi)\mathrm{d}\phi\]

Where \(G'\) is the derivative of the circulation strength of the bound vortex. It has the form:

(32)\[ \begin{align}\begin{aligned}G'(\phi) &= \sum_{n=1}^N G_n h_n(\phi)\\G(\phi) &= \sum_{n=1}^N G_n f_n(\phi)\end{aligned}\end{align} \]

where

(33)\[ \begin{align}\begin{aligned}h_n(\phi) &= \frac{2}{1 + N}\sum_{k=1}^N k\sin(k\phi_n)\cos(k\phi)\\f_n(\phi) &= \frac{2}{1 + N}\sum_{k=1}^N \sin(k\phi_n)\sin(k\phi)\end{aligned}\end{align} \]

The functions \(P\) and \(R\) are functions of the blade geometry

(34)\[ \begin{align}\begin{aligned}P(\phi_\nu, \phi) &= \frac{1}{\cos(\phi_\nu) - \cos(\phi)}\left[\frac{\xi(\phi_\nu) - \xi(\phi) + \frac{1}{2}}{\sqrt{\left[\xi(\phi_\nu) - \xi(\phi) + \frac{1}{2}\right]^2 + \left[R/c(\phi_\nu)\right]^2[\cos(\phi_\nu) - \cos(\phi)]^2}} - 1\right]\\R(\phi_\nu, \phi) &= \frac{\xi(\phi_\nu) - \xi(\phi) + \frac{1}{2} + \xi'(\phi)(\cos(\phi) - \cos(\phi_\nu))}{\left(\left[\xi(\phi_\nu) - \xi(\phi) + \frac{1}{2}\right]^2 + \left[R/c(\phi_\nu)\right]^2[\cos(\phi_\nu) - \cos(\phi)]^2\right)}\end{aligned}\end{align} \]

where \(\xi(\phi) = x(\phi)/\bar c\) is the \(x\)-coordinate position of the blade quarter chord line normalized by the average blade chord and \(\xi'(\phi)\) is its first spatial derivative.

A system of equations can be assembled by plugging in equation (32) into equation (31) and applying trapezoidal integration to the second two integrals. The first integral can be solved analytically. This results in the following equation:

(35)\[\begin{split}\begin{split} \alpha(\phi_\nu) =& \frac{1}{1 + N}\sum_{n=1}^N G_n \sum_{k=1}^N \frac{k\sin(k\phi_n)\sin(k\phi_\nu)}{\sin(\phi_\nu)} \\ &- \frac{1}{4(M + 1)}\sum_{n=1}^N G_n \left[\frac{P(\phi_\nu, 0)h_n(0) - P(\phi_\nu, \pi)h_n(\pi)}{2} + \sum_{\mu=1}^{M}P(\phi_\nu, \phi_\mu)h_n(\phi_\mu)\right] \\ &+ \frac{1}{4(M + 1)}\left(\frac{R}{c(\phi_\nu)}\right)^2\sum_{n=1}^N G_n \sum_{\mu=1}^M R(\phi_\nu, \phi_\mu) f_n(\phi_\mu)\sin(\phi_\mu) \end{split}\end{split}\]

where \(M\) is the number of integral points to compute and

(36)\[\phi_\mu = \frac{\mu\pi}{M + 1}\]

This can now be assembled as a system of linearly independent equations in the form

(37)\[\boldsymbol\alpha = \boldsymbol A \boldsymbol G\]

Where the coefficient matrix \(\boldsymbol A\) is built as such:

(38)\[\begin{split}\begin{split} \boldsymbol A =& \frac{1}{1 + N} \sum_{k=1}^N \frac{k\sin(k\phi_n)\sin(k\phi_\nu)}{\sin(\phi_\nu)} \\ &- \frac{1}{4(M + 1)}\left[\frac{P(\phi_\nu, 0)h_n(0) - P(\phi_\nu, \pi)h_n(\pi)}{2} + \sum_{\mu=1}^{M}P(\phi_\nu, \phi_\mu)h_n(\phi_\mu)\right] \\ &+ \frac{1}{4(M + 1)}\left(\frac{R}{c(\phi_\nu)}\right)^2\sum_{\mu=1}^M R(\phi_\nu, \phi_\mu) f_n(\phi_\mu)\sin(\phi_\mu) \end{split}\end{split}\]

After solving the system the vector \(\boldsymbol G = G_n\) contains the normalized circulation strength at the blade colocation points. Once the circulation on the blade is known, the tip vortex strength is estimated to be the maximum \(G_n\):

(39)\[\Gamma_v = u_t R \max(G_n)\]

The strength of the shed wake can also be computed by taking the difference of the current blade bound circulation strength and the previous timestep’s circulation strength. This ensures that the total circulation is conserved. Further it accounts for the unsteady loading on the blade when the induced inflow contribution from the shed wake is computed.

Vortex Core Model

The vortex core model used in this work is based on the model developed by Ramasamy and Leishman [RL07]. This is a semi-empirical model that takes into account the stretching a vortex filament undergoes as it evolves in time. As presented by Ramasamy and Leishman, it takes the following form

(40)\[r_c = \sqrt{r_0^2 + 4\alpha_L \nu \delta \int_{\zeta_0}^{\zeta} \frac{1}{\Omega(1 + \varepsilon)}\mathrm{d}\zeta}\]

where \(r_0\) is the initial core radius, \(\alpha_L\) is Lamb’s constant (\(\alpha_L = 1.25643\)), \(\nu\) is the kinematic viscosity, and \(\varepsilon\) is the filament strain:

(41)\[\varepsilon = \frac{\Delta l}{l}\]

The parameter \(\delta\) is the ratio of apparent viscosity to actual viscosity. This accounts for core growth due to turbulent eddies present in the vortex core. Squire [Squ65] assumed that the eddy viscosity is only function of the kinematic viscosity and vortex circulation yielding the following expression for \(\delta\)

(42)\[\delta = 1 + a_1\left(\frac{\Gamma_v}{\nu}\right) = 1 + a_1 Re_v\]

The constant \(a_1\) is another empirical parameter that depends on the way eddy viscosity varies radially in the vortex core. Ramasamy and Leishman [RL06] found it lies in a broad range from \(5\times 10^{-5}\) to \(4\times 10^{-4}\). A value of \(a_1 = 6.5\times 10^-5\) was used in this work.

To avoid integrating the strain integral over the wake age every timestep, the integral is transformed into an ODE that can be time-marched. The integral is assigned a variable \(\eta\) and the wake age derivative is taken:

(43)\[ \begin{align}\begin{aligned}\eta &= \int_{\zeta_0}^{\zeta} \frac{1}{\Omega(1 + \varepsilon)}\mathrm{d}\zeta\\\frac{\mathrm{d}\eta}{\mathrm{d}\zeta} &= \frac{1}{\Omega(1 + \varepsilon)}\end{aligned}\end{align} \]

To make it a time derivative instead of a wake age derivative the following is used:

(44)\[\frac{\mathrm{d}\eta}{\mathrm{d}t} = \frac{\mathrm{d}\eta}{\mathrm{d}\zeta}\frac{\mathrm{d}\zeta}{\mathrm{d}t}\]

The change of wake age in time is just the rotor speed \(\Omega\) and we can plug in equation (43) to yield:

(45)\[\frac{\mathrm{d}\eta}{\mathrm{d}t} = \frac{1}{(1 + \varepsilon)}\]

This ODE can now be time-marched, solving for \(\eta\), which reduces the cost of computing core growth.

Wake Model

Fig. 2 A component view of both the trailed tip wake and shed wake vortex filaments of a single blade

The rotor wake consists of two components: the shed wake and the trailed tip wake. The shed wake accounts for the circulation left behind a blade that undergoes a change in circulation. The trailed tip wake accounts for the tip vortex rolled up at the end of blade. Both components of the wake are comprised of straight vortex filament segments that are emitted from the blade at every timestep. Fig. 2 shows a breakdown of the wake components for a single blade. The trailed tip vortex elements are advected downstream by the induced velocities computed from each rotor’s dynamic inflow model and the freestream velocity. As the dynamic inflow model is derived from potential flow theory, the flow field of each rotors dynamic inflow model can be superimposed to get the aggregate flow field of the whole system. The shed wake is advected by the induced velocity of its own rotor’s dynamic inflow model and freestream velocity allowing it to deform as well. Thus the total induced velocity acting upon a trailed vortex element is:

(46)\[\boldsymbol v_i(x, y, z) = \sum_{j = 0}^{N_R} V_j(x, y, z)\]

and for a shed vortex element

(47)\[\boldsymbol v_i(x, y, z) = V_r(x, y, z)\]

where \(r\) in this case is index of the rotor that shed the vortex element. This gives the final expression for the advection velocity of each vortex element

(48)\[\boldsymbol v_{\mathrm{ve}} = \boldsymbol{\dot x}_{\mathrm{ve}} = \boldsymbol{v}_i + \boldsymbol{V}_\infty\]

Which can the be time-marched [GSS15] to yield to wake geometry. This can be integrated using the same time-marching integrator as the systems in equations (8) and (9). This is the same basic concept behind Beddoes’ wake model [Bed85], but using a dynamic inflow model to solve for the velocity induced by a nonuniformly loaded actuator disk in place of an empirical inflow model.

The induced velocity that is used to compute the blade loading is computed from the both the shed and trailed tip wake vortex filaments. This is computed via the Biot-Savart law with a finite vortex core radius to remove singularities:

(49)\[\boldsymbol v_v = \frac{\Gamma_v}{4\pi(h^2 + r_c^2)}\left[\hat l_v\cdot(\hat r_1 - \hat r_2)\right]\left(\hat r_1 \times \hat r_2\right)\]

where \(h\) is the distance from the vortex element to the point of interest, \(\hat l_v\) is the normalized element length vector, \(\hat r_1, \hat r_2\) are the normalized vectors between the element endpoints and the point of interest, and \(r_c\) is the vortex core radius. Thus the full Biot-Savart integration for any point on some rotor \(r\) is:

(50)\[\boldsymbol v_{i_r}(x, y, z) = \sum_{j = 0}^{N_R} \sum_{k = 0}^{N_B} \sum_{l=0}^{N_v} \boldsymbol v_v(\Gamma_{v_{j_{k_l}}}, x, y, z) + \sum_{k = 0}^{N_B} \sum_{s = 0}^{N_S} \sum_{l = 0}^{N_v} \boldsymbol v_v(\Gamma_{v_{r_{k_{s_l}}}}, x, y, z)\]

where \(N_R\) is the number of rotors in the system, \(N_B\) is the number of blades on the rotor, \(N_v\) is the number of elements in the vortex filament, and \(N_S\) is the number of shed vortex filaments.

Bed85(1,2)

TS Beddoes. A wake model for high resolution airloads. In International Conference on Rotorcraft Basic Research. 1985.

GSS15

Eric Greenwood, Fredric H. Schmitz, and Richard D. Sickenberger. A Semiempirical Noise Modeling Method for Helicopter Maneuvering Flight Operations. Journal of the American Helicopter Society, 60(2):1–13, April 2015.

GPP21(1,2)

Feyyaz Guner, JVR Prasad, and David A Peters. An Approximate Finite State Dynamic Wake Model for Predictions of Inflow below the Rotor. Journal of the American Helicopter Society, 2021.

HPP15

Jianzhe Huang, David A. Peters, and J. V. R. Prasad. Converged Velocity Field for Rotors by a Blended Potential Flow Method. GSTF Journal on Aviation Technology (JAT), April 2015.

MP02

Jorge A Morillo and David A Peters. Velocity field above a rotor disk by a new dynamic inflow model. Journal of aircraft, 39(5):731–738, 2002.

RL06

Manikandan Ramasamy and J Gordon Leishman. A generalized model for transitional blade tip vortices. Journal of the American Helicopter Society, 51(1):92–103, 2006.

RL07

Manikandan Ramasamy and J Gordon Leishman. A reynolds number-based blade tip vortex model. Journal of the American Helicopter Society, 52(3):214–223, 2007.

RG22

Robert F Rau and Eric Greenwood. A dynamic wake model for multirotor aircraft. In AIAA AVIATION 2022 Forum, 3959. 2022.

Squ65

HB Squire. The growth of a vortex in turbulent flow. Aeronautical Quarterly, 16(3):302–306, 1965.

vdW12

Berend G van der Wall. Extensions of prescribed wake modelling for helicopter rotor BVI noise investigations. CEAS Aeronautical Journal, 3(1):93–115, 2012.

Wei47

Jo Weissinger. The lift distribution of swept-back wings. Technical Report, NACA, 1947.

WG07

Adam M Wickenheiser and Ephrahim Garcia. Aerodynamic modeling of morphing wings using an extended lifting-line analysis. Journal of Aircraft, 44(1):10–16, 2007.