*3.1. Numerical Methods*

The fluid–structure interaction solution method in conjunction with the weak coupling solution technique were used to numerically calculate the evolutions of wind field over the computational domain and to obtain the aerodynamic force acting on the harvester. The two-dimensional Reynolds number average Navier–Stokes (RANS) equation and the k-ω SST (shear stress transport) turbulent model [26] were used to numerically simulate the flow around a square cross-section. This turbulence model was chosen because, compared to other turbulence models, the k-ω model of SST is a hybrid of the k-ω and k-ε models. It activates the k-ω model near the wall and the k-ε model in free flow. It has better performance when predicting the boundary layer flow of the reverse pressure gradient. Therefore, this paper adopts the SST k-ω turbulence model in CFD simulation.

The interaction effects between the wind flows and the harvester structure were embedded into the CFD to implement the link with the main program. At the same time, the structure dynamic equation was solved and the structure domain calculated, and then the motion lift coefficient and displacement of the structure were obtained. At the same time, the overset grid technology was used to update the grid of the computational domain in real time.

User-defined functions (UDFs) were compiled (in C Language) and embedded into the CFD software to realize the simulation of fluid–structure interaction. In this paper, DEFINE\_CG\_MOTION was used as the motion function to define the motion mode. The overset grid model was used to define computational grid movement and data exchange. The Newmark-β method was used in UDF (refer to the Appendix A for details). This method is a numerical integration method for solving differential equations. It has been widely used to solve problems with oscillating systems [27].

Based on the basic assumptions of the central difference method, the explicit difference method of dynamic equations was adopted. The basic assumption of the Newmark-β method is that the acceleration changes linearly in each time increment, and the characteristics of the system remain constant during this interval. The equilibrium requirement of the force system acting on the mass at any time is as follows: the Newmark-β method approximates the speed and displacement of the system at time (*t* + Δ*t*) through the following two assumptions:

$$
\dot{y}(t + \Delta t) = \dot{y}(t) + (1 - \gamma)\Delta t \ddot{y}(t) + \gamma \Delta t \ddot{y}(t + \Delta t) \tag{4}
$$

$$y(t + \Delta t) = y(t) + \Delta t \dot{y}(t) + (\frac{1}{2} - \beta)\Delta t^2 \ddot{y}(t) + \beta \Delta t^2 \ddot{y}(t + \Delta t) \tag{5}$$

The two parameters *γ* and *β* in the formula can be selected as required, and different combinations correspond to different processing effects. When *γ* = <sup>1</sup> <sup>2</sup> and *<sup>β</sup>* <sup>=</sup> <sup>1</sup> <sup>2</sup> , it is the central difference method [28]. Therefore, the acceleration and velocity can be obtained.

$$\ddot{y}(t+\Delta t) = \frac{1}{\beta \Delta t^2} (y(t+\Delta t) - y(t)) - \frac{1}{\beta \Delta t} \dot{y}(t) - (\frac{1}{2\beta} - 1)\ddot{y}(t) \tag{6}$$

$$\dot{y}(t+\Delta t) = \frac{\gamma}{\beta \Delta t}(y(t+\Delta t) - y(t)) + (1 - \frac{\gamma}{\beta})\dot{y}(t) + (1 - \frac{\gamma}{2\beta})\Delta t \ddot{y}(t) \tag{7}$$

The motion of the system at time *t* + Δ*t* can be expressed as:

$$M\ddot{y}(t+\Delta t) + C\dot{y}(t+\Delta t) + Ky(t+\Delta t) = F\_{\mathcal{Y}}(t+\Delta t) \tag{8}$$

Integrating the above formulas leads to:

$$
\tilde{K}y(t+\Delta t) = \tilde{F}\_{\mathcal{Y}}(t+\Delta t) \tag{9}
$$

where

$$\overline{K} = \frac{1}{\beta \Delta t^2} M + \frac{\gamma}{\beta \Delta t} C + K \tag{10}$$

$$\begin{split} \overline{F}(t + \Delta t) &= F(t + \Delta t) + M \Big[ \frac{1}{\overline{\beta}\Delta t^2} y(t) + \frac{1}{\overline{\beta}\Delta t} \dot{y}(t) + (\frac{1}{2\overline{\beta}} - 1)\ddot{y}(t) \Big] \\ &+ \mathbb{C} \Big[ \frac{\gamma}{\overline{\beta}\Delta t} y(t) + (\frac{\gamma}{\overline{\beta}} - 1)\dot{y}(t) + (\frac{\gamma}{2\overline{\beta}} - 1)\Delta t \ddot{y}(t) \Big] \end{split} \tag{11}$$

The numerical simulation flowchart is shown in Figure 3.

**Figure 3.** Flowchart of transient numerical simulation.
