## 3 Multi-Physical Coupled Fields Simulation Analyses
### Equations for Rigid Body Motion and Solution Methodology
The Runge-Kutta algorithm is implemented as follows:  
\[y_{n+1}=y_{n}+\dot{y}\Delta t+\frac{\Delta t^{2}}{12}(S_{1}+2S_{2}+2S_{3}+S_{4}) \tag{6}\]  
\[S_{1}=f(x_{n},y_{n}) \tag{7}\]  
\[S_{2}=f(x_{n+1/2},y_{n}+\frac{\Delta t}{2}S_{1}) \tag{8}\]  
\[S_{3}=f(x_{n+1/2},y_{n}+\frac{\Delta t}{2}S_{2}) \tag{9}\]  
\[S_{4}=f(x_{n+1},y_{n}+\Delta tS_{3}) \tag{10}\]  
where \(x_{n}\) represents a series of discrete nodes; \(y_{n}\) represents the approximate solution corresponding to \(x_{n}\). Since the airfoil-based flutter energy harvester system is 0 at the initial position h, \(\dot{h}\), \(\dot{a}\), \(\alpha\) and \(V\), the initial value problem of the second-order differential Equations (1)-(3) has the following expression  
\[\ddot{h}= \left[(L+V\theta)I_{n}/S_{n}-(M-C_{n}\dot{\alpha}-(k_{x0}+k_{n2} \alpha^{2})\alpha)\right]/(m_{T}I_{n}/S_{n}-S_{n}) \tag{11}\] \[-I_{a}/(m_{T}I_{a}-S_{n}^{2})\cdot(C_{n}\dot{h}+(k_{x0}+k_{n2}h^ {2})h)\] \[\ddot{\alpha}= \left[Mm_{T}/S_{n}-(L-C_{n}\dot{h}-(k_{x0}+k_{x2}h^{2})h-\theta V )\right]/(m_{T}I_{a}/S_{n}-S_{n})\] (12) \[-m_{T}/(m_{T}I_{a}-S_{n}^{2})\cdot(C_{a}\dot{\alpha}+(k_{x0}+k_{x 2}\alpha^{2})\alpha)\]  
\[y_{n+1}=y_{n}+\dot{y}\Delta t+\frac{\Delta t^{2}}{12}(S_{1}+2S_{2}+2S_{3}+S_{4}) \tag{6}\]0  
As depicted in Figure 3, the solution process proceeds as follows: (1) initialize the flow field; (2) use UDFs to calculate amplitude, velocity, angular velocity, and angular acceleration; (3) update the dynamic mesh; (4) when the calculation is completed, return to step (2) until the calculation is completed. \[y_{n+1}=y_{n}+\dot{y}\Delta t+\frac{\Delta t^{2}}{12}(S_{1}+2S_{2}+2S_{3}+S_{4}) \tag{6}\]1 The process of the numerical simulation. The numerical simulation domain consists of a rectangular region measuring 20b \(\times\) 12b, as illustrated in Figure 4a. The boundary conditions are configured as follows: velocity inlet at the upstream boundary; pressure outlet at the downstream boundary; symmetric conditions at both upper and lower boundaries; and no-slip wall condition on the airfoil-based surface. To ensure mesh quality and capture flow characteristics accurately, a structured mesh strategy is employed, as shown in Figure 4b. Special attention is paid to the boundary layer region around the airfoil, where mesh refinement is implemented to resolve the near-wall flow physics. To handle the airfoil-based motion and prevent mesh distortion, a dynamic mesh technique is adopted, ensuring that the boundary layer mesh moves synchronously with the airfoil-based surface. The mesh update is performed at each time step according to the displacement calculated through the UDF framework. For spatial discretization, the Green-Gauss node-based gradient method available in FLUENT is selected. This approach demonstrates superior accuracy for unstructured meshes, particularly in cases involving triangular elements and significant mesh size variations. The method proves especially effective in handling mesh deformation during boundary motion simulations. In this study, ICAM CFD is utilized to generate a grid for computations. The grid's Element Quality surpasses 0.4, while the skewness remains below 0.6. To guarantee calculation accuracy and efficiency, we examine grid independence for the airfoil-based flutter energy harvesting system at a mounting angle of 0\({}^{\circ}\) under a wind speed of 9 m/s. Three varying grid resolutions are established sequentially, with an increase in the number of grids. The first boundary layer mesh height strictly adheres to the criteria of \(y^{+}\leq 30\), where the boundary layer mesh growth rate is set at 1.2.