An immersed lifting and dragging line model for the vortex particle-mesh method

Abstract

The numerical study of the wake of full-scale slender devices such as aircraft wings and wind turbine blades requires high-fidelity large eddy simulation tools. The broad spectrum of scales involved entails the use of coarse-grain models for the devices. Actuator disk or line methods have been developed for that purpose, and are to date the most employed in that context. These methods transfer the force from the device to the fluid using 3D mollification functions. A novel immersed lifting and dragging line method is here presented, together with its implementation in a hybrid vortex particle-mesh flow solver. The line models the effect of blades or wings on the flow through the generation of vorticity, in a Lagrangian manner. This has several advantages over the treatment inherent to actuator methods: The absence of mollification in the spanwise direction and the relaxation of the classical Courant–Friedrichs–Lewy condition contribute to keeping the accuracy and the efficiency at a high level. The method is thoroughly verified against theory using results on various airfoils and wings. Finally, the efficiency of the approach is illustrated with the simulation of the wake of an elliptical wing. The role of parasitic drag in the development of wake instabilities is pointed out, by comparing configurations with none, moderate and high profile drag. In the last case, the simulation captures the fine scales vortex dynamics and the establishment of a turbulent wake over a length of 30 wing spans.

Appendices

Appendix A. Reconstruction of the effective upstream velocity

In the presence of parasitic drag, the effective upstream velocity $\mathbf {U}_e$ differs from the measured $\mathbf {U}_{0}$. We here obtain the velocity induced by the dragging vortex sheets, which will be used as a correction velocity $U_\text {corr}$ to recover the effective upstream velocity in the ILDL model: $U_e = U_\text {0}-U_\text {corr}$.

Fig. 16

Dragging dipole model centered at O and of width d

First, we consider a 2D steady, inviscid wake model with singular shed vortex sheets (Fig. 16) and an upstream velocity $U_e$. With small drag forces such as those experienced on airfoils, we reasonably assume that these vortex sheets remain parallel. Using the Biot–Savart law, the velocity evaluated at the center O of the dipole is

$$\begin{aligned} U_0 = U_e - \frac{\gamma }{2} = U_e - \frac{\mu }{2d}. \end{aligned}$$

(17)

Using this expression to substitute $\mu$ in Eq. (10), we obtain

$$\begin{aligned} \left( \frac{U_0}{U_e}\right) ^2 - \frac{U_0}{U_e} + \frac{C_d}{4} \frac{c}{d} = 0, \end{aligned}$$

(18)

which provides the velocity correction factor

$$\begin{aligned} \frac{U_e}{U_0} = \frac{2}{1+\sqrt{1-C_d \frac{c}{d}}}. \end{aligned}$$

(19)

The factor $C_d \frac{c}{d}$ plays an important role in this expression, and it is clear that it should be maintained well below unity for the correction to be valid.

In the ILDL, the vortex sheets are regularized using a 1D Gaussian regularization (of parameter $\sigma _d$) along y. The drag-induced mollified vorticity field is

$$\begin{aligned} \omega _\sigma ^{d}(x,y)&= \frac{\partial V_\sigma }{\partial x} - \frac{\partial U_\sigma }{\partial y} \nonumber \\&= \frac{\gamma }{\sqrt{\pi }\,\sigma _d} H(x) \left( \exp \left( { - \frac{\left( {y}+\frac{d}{2}\right) ^2}{{\sigma _d}^2}}\right) - \exp \left( { - \frac{\left( y-\frac{d}{2}\right) ^2}{{\sigma _d}^2}}\right) \right) , \end{aligned}$$

(20)

where H(x) is the Heavyside step function. As $\frac{\partial V_\sigma }{\partial x}$ is small, we obtain the horizontal velocity at $x=0$ as

$$\begin{aligned} U_\sigma ^d(0,y) = - \int _{-\infty }^{y} \omega _\sigma (0,\zeta ) d\zeta = -\frac{ \gamma }{4} \, \left[ {{\,\mathrm{erf}\,}}\left( \frac{y+\frac{d}{2}}{\sigma _d} \right) - {{\,\mathrm{erf}\,}}\left( \frac{y-\frac{d}{2}}{\sigma _d} \right) \right] , \end{aligned}$$

(21)

In case of pointwise velocity sampling, the correction $U_\text {corr}$ is based on the evaluation of $U_\sigma ^d$ at $y=0$, and

$$\begin{aligned} U_e = U_0 - U_\sigma ^d(0,0) = U_0 + \frac{\mu }{2d} \, {{\,\mathrm{erf}\,}}\left( \frac{d}{2\sigma _d} \right) , \end{aligned}$$

(22)

and we obtain, again using Eq. (10),

$$\begin{aligned} \frac{U_e}{U_0} = \frac{2}{1+ \sqrt{1- C_d \frac{c}{d} {{\,\mathrm{erf}\,}}\left( \frac{d}{2\sigma _d} \right) }}. \end{aligned}$$

(23)

This is consistent with Eq. (19) when $\sigma _d \rightarrow 0$. We also make sure, in our simulations, that $C_d \frac{c}{d}$ is well below unity.

In the present work, we use the “integral velocity sampling” [3], i.e., the principal value obtained as $U_0 = \int \xi _{\sigma _l}(\mathbf {x}) U(\mathbf {x}) d\mathbf {x}$, where $\xi _{\sigma _l}$ is the kernel used for the mollification of the bound vorticity, $\xi _{\sigma _l}(y)= \frac{1}{\sqrt{\pi }} \frac{1}{\sigma _l} \exp \left( {-\frac{y^2}{\sigma _l^2}}\right)$. Therefore, $U_\sigma ^d$ is used to derive an integrated correction,

$$\begin{aligned} U_e = \frac{1}{\sqrt{\pi }} \int _{-\infty }^{\infty } \exp \left( {-\frac{y^2}{\sigma _l^2}}\right) \, \left( U(0,y) - U_\sigma ^d(0,y)\right) \frac{dy}{\sigma _l}. \end{aligned}$$

(24)

In practice, this integral is evaluated numerically using the weights $w_j$ of Table 1, where U is the velocity measured at the location of the bound vortex particles.

In general, $\mathbf {U}_0$ is not aligned with $\mathbf {U}_e$; the correction velocity should be aligned with the local shedding direction as shown in Fig. 2, and we use Eq. (24) in its vector form.

Appendix B. Errors related to the ILDL model assumptions

1.1 B.1 Error in drag evaluation, due to wake expansion

The expansion of the wake is neglected in the derivation of Eq. (8). This assumption generates a small error which can be quantified using Eq. (13).

Again, we use our simplified drag model made of two parallel singular vortex sheets (Fig. 16). As the vortex sheets are horizontal, the vertical velocity V across them is nonzero. Consequently, there is a contribution to drag from the integral of $\mathbf {U}\wedge \varvec{\omega }$. This contribution can be assimilated to the error introduced by neglecting the expansion of the wake.

The vertical velocity induced by the two vortex sheets is

$$\begin{aligned} V^d(x,y)&= \frac{\gamma }{2\pi } \int _{0}^\infty \left[ \frac{(x-x')}{\left( (x-x')^2 + (y-\frac{d}{2})^2 \right) } - \frac{(x-x')}{\left( (x-x')^2 + (y+\frac{d}{2})^2 \right) } \right] \mathrm{d}x' \nonumber \\&= -\frac{\gamma }{4\pi } \log \left( \frac{x^2+ (y-\frac{d}{2})^2}{ x^2 + (y+\frac{d}{2})^2 } \right) . \end{aligned}$$

(25)

The model error in drag thus amounts to

$$\begin{aligned} \int \mathbf {U}\wedge \varvec{\omega }\, d\mathbf {x}\cdot \hat{\mathbf{e}}_x = -\int _0^\infty V^d\gamma \,\mathrm{d}x=\frac{\gamma ^2}{2\pi } \int _0^\infty \log \left( \frac{x^2}{x^2+d^2} \right) \mathrm{d}x = -\frac{\mu \gamma }{2}. \end{aligned}$$

(26)

This term is proportional to the square of $\mu$, and it is therefore indeed negligible for small drags.

1.2 Drag-induced lift

In the ILDL model, it is implicitly assumed that the effects of lift and drag can be added by linear superposition. In fact, even with wake expansion neglected, an interaction exists, which can be quantified again using the simple model of Fig. 16, to which we add a singular bound vortex of circulation $\varGamma$ at O (similar to what is presented in Fig. 2a).

As can be seen from Eq. (13), the integral of $\mathbf {U}\wedge \varvec{\omega }$ will have a nonzero contribution from the combination of the axial velocity induced by the bound vortex and the circulation of the vortex sheets. The horizontal velocity induced by the bound vortex is

$$\begin{aligned} U^l(x,y) = \frac{\varGamma }{2\pi } \frac{y}{(x^2+y^2)}. \end{aligned}$$

(27)

The interaction results in the appearance of some lift induced by drag,

$$\begin{aligned} \int \mathbf {U}\wedge \varvec{\omega }\, \mathrm{d}\mathbf {x}\cdot \hat{\mathbf{e}}_y = - 2 \int _0^\infty U^l\left( x,\frac{d}{2}\right) \, \gamma \, \mathrm{d}x = - \frac{\varGamma \gamma }{\pi } \int _0^\infty \frac{\frac{d}{2}}{\left( x^2+(\frac{d}{2})^2\right) } \mathrm{d}x = -\frac{\varGamma \gamma }{2}. \end{aligned}$$

(28)

When the mollification of the bound vortex and of the shed vortex sheets is taken into account, one can show that only a fraction of this term contributes to the drag-induced lift. The fraction depends on the values of $\frac{\sigma _l}{d}$ and $\frac{\sigma _d}{d}$, and it can be obtained numerically using the weights $w_i$.

The drag-induced lift, as computed here, could thus be accounted for in the ILDL model by modifying Eq. (7). However, in this work, it was decided to keep the model in its simplest form. The error associated with this simplification is quantified in the analysis of the validation results.

Appendix C. Fourier-based Poisson solver for domains with unbounded and inflow/outflow directions

In this work, we rely on the Fourier-based solvers used by [5, 8] in order to solve Eq. (2) for the velocity. This solver handles a combination of periodic and unbounded directions through a Hockney–Eastwood approach [17] . We here extend that work to handle the unbounded directions together with the inflow and outflow conditions in the flow stream direction. Following [8], we re-write Eq. (2) as

$$\begin{aligned} \nabla ^2_{yz} \mathbf {u}+ \frac{\partial ^2 \mathbf {u}}{\partial x^2} = -\nabla \wedge \varvec{\omega }, \end{aligned}$$

(29)

where x is assumed to be the inflow/outflow direction and $\nabla ^2_{yz}$ is the Laplacian operator restricted to the $y-z$ plane. The Fourier-transform along x yields

$$\begin{aligned} \nabla ^2_{yz} \tilde{\mathbf {u}} -k_x^2 \tilde{\mathbf {u}} = -\widetilde{\nabla \wedge \varvec{\omega }}, \end{aligned}$$

(30)

where $\tilde{\mathbf {u}}(k_x,y,z)$ stands for the x-transformed field. For a given $k_x$ mode, this is a two-dimensional Helmholtz equation in an unbounded (y, z)-domain, which we can solve through the convolution with the Green’s function

$$\begin{aligned} G^{\text {2D}}_{\mathrm {i}|k_x|}(y,z) = {\left\{ \begin{array}{ll} \frac{1}{2\pi } \ln (\sqrt{y^2+z^2})&{} \text {if } k_x = 0\,,\\ \frac{1}{2\pi } K_0\left( |k_x| \, \sqrt{y^2+z^2}\right) &{} \text {otherwise}, \end{array}\right. } \end{aligned}$$

(31)

where $K_0$ denotes the modified Bessel function of the second kind. These expressions are singular at $(y,z)=(0,0)$. This is here circumvented by replacing them with the corresponding cell-averaged value [8], but we could alternatively use a mollified Helmholtz kernel [15, 16, 41]. The convolution itself can be carried out in Fourier space through the domain-doubling approach of [17].

In such a context, the imposition of inflow/outflow conditions simply corresponds to a restriction of the Fourier modes to be considered, depending on the symmetry conditions that need to be enforced on each side of the domain. We assume the domain to go from 0 to $L_x$ in x. Considering the velocity field, we wish to impose the following conditions on the streamwise velocity: $u_x(0,y,z) = u_x(L_x,y,z) = U_\infty$. This corresponds to odd symmetry conditions on $\mathbf {u}_{\varvec{\omega }}$ at the inflow and at the outflow. These are completed with even symmetry conditions on the transverse components $\frac{\partial u_y}{\partial x} (0,y,z) = \frac{\partial u_z}{\partial x} (0,y,z) = 0$ and $\frac{\partial u_y}{\partial x} (L_x,y,z) = \frac{\partial u_z}{\partial x} (L_x,y,z) = 0$. The allowed modes thus have to be harmonics of a base mode of period $2L_x$ with the proper phase

$$\begin{aligned} u_x&= U_\infty + \sum _i \sin \left( i\pi \frac{x}{ L_x} \right) \tilde{u}_x\left( \frac{i\pi }{L_x},y,z\right) \;, \end{aligned}$$

(32)

$$\begin{aligned} u_{(y,z)}&= \sum _i \cos \left( i\pi \frac{x}{ L_x} \right) \tilde{u}_{(y,z)}\left( \frac{i\pi }{L_x},y,z\right) \;. \end{aligned}$$

(33)

The modes of the corresponding vorticity field are then

$$\begin{aligned} \omega _x&= \sum _i \cos \left( i\pi \frac{x}{ L_x} \right) \tilde{\omega }_x\left( \frac{i\pi }{L_x},y,z\right) \;, \end{aligned}$$

(34)

$$\begin{aligned} \omega _{(y,z)}&= \sum _i \sin \left( i\pi \frac{x}{ L_x} \right) \tilde{\omega }_{(y,z)}\left( \frac{i\pi }{L_x},y,z\right) \;. \end{aligned}$$

(35)

Notice that, in general, any combination of symmetry conditions is accessible and it leads to domains with a periodicity of either 2 or 4 times $L_x$.

The practical implementation of this restriction is made straightforward through the use of discrete sine and cosine transforms (DCT and DST) in our FFT-based solver. Advantageously, and unlike the domain-doubling technique for the unbounded directions, this implementation does not lead to a twofold (or fourfold) increase in the unknowns as the DCT/DST transforms are performed in the domain with its original $L_x$ size.