## 3 Multi-Physical Coupled Fields Simulation Analyses  
### Equations for Rigid Body Motion and Solution Methodology  
In this study, a dual-degree-of-freedom model is employed to characterize the system's motion. Given that the operational velocity of the system does not exceed 0.1 Mach, air is assumed to behave as an incompressible fluid. The numerical simulations are conducted using ANSYS FLUENT 2020R2. The computer's CPU model utilized for the simulations is the Intel Xeon Silver 4216 with 16 cores (Intel, Santa Clara, CA, USA). Since flow-induced vibration represents a typical transient problem, it is essential to account for the flow field phenomenon at each moment within a time period. Consequently, the unsteady Navier-Stokes equations are utilized to describe the flow. The control equation is expressed in dimensionless form as follows:  
\[\frac{\partial\overline{u_{i}}}{\partial\overline{x_{i}}}=0 \tag{4}\]  
\[\rho\frac{\partial\overline{u_{i}}}{\partial t}+\rho\frac{\partial\overline{u_ {i}u_{j}}}{\partial t}=-\frac{\partial p}{\partial x_{i}}+\mu\
abla^{2} \overline{u_{i}}-\rho\frac{\overline{u_{i}^{\prime}u_{j}^{\prime}}}{\partial x _{j}} \tag{5}\]  
where \(x_{i}\) and \(x_{j}\) are the Cartesian coordinates in the \(i\) and \(j\) directions, \(\overline{u_{i}}\) and \(\overline{u_{i}^{\prime}}\) are the time-averaged values of the corresponding instantaneous components of velocity and velocity wave term, \(\rho\) is the density of the fluid, \(\mu\) is the dynamic viscosity of the fluid, \(p\) is the pressure, and \(\rho\overline{u_{i}^{\prime}u_{j}^{\prime}}\) is the Reynolds stress tensor. To enhance the accuracy of the computational results, we employ the turbulence model proposed by [PERSON] [31; 32] to address the issue. This model incorporates a turbulent viscosity limit, which enhances the accuracy of calculations, particularly for backflow scenarios [34; 35]. The pressure-velocity coupling in this study is achieved through the Coupled algorithm in FLUENT. This solver employs a fully implicit approach where pressure and velocity fields are solved simultaneously. Compared to the segregated solver, the Coupled algorithm demonstrates superior convergence characteristics and enhanced numerical stability. Although it requires more computational memory per iteration, this approach significantly reduces the total number of iterations needed for convergence, thereby improving overall computational efficiency. In this paper, FLUENT User Defined Functions (UDFs) are employed to establish mesh motion, while the fourth-order Runge-Kutta method is utilized to solve the motion dynamics of the airfoil-based flutter energy harvester system. This approach enables the \begin{table}
\begin{tabular}{l c} \hline  & **Properties** & **Value** \\ \hline Span, \(s\) (m) & 0.1 \\ Total mass, \(\text{m}_{\text{I}}\) (kg) & 0.6 \\ Airfoil-based mass, \(\text{m}_{f}\) (kg) & 0.18 \\ Semi-chord, \(b\) (m) & 0.125 \\ Position of the airfoil-based pitch relative to the & 0.5 \\ semi-chord, \(a\) & 0.5 \\ Offset of airfoil-based gravity center, \(X_{a}\) & 0.335 \\ Stiffness coefficient in plunge, \(\text{k}_{\text{h}0}\) (N/m) & 240 \\ Stiffness coefficient in pitch, \(\text{k}_{\text{a}0}\) (0.28 Nm) & 0.28 \\ Damping coefficient in plunge, \(\text{C}_{\text{h}}\) (kg/s) & 0.056 \\ Damping coefficient in pitch, \(\text{C}_{\text{a}}\) (kg/s) & 0.056 \\ Resistance, \(R\) (K\(\Omega\)) & 250 \\ Capacitance, \(\text{C}_{p}\) (F) & 0.0000007 \\ electromechanical coefficient of the harvester, \(\theta\) (N/V) & \(2.7\times 10^{-4}\) \\ \hline \end{tabular}
\end{table} 
Table 2: Dimensional parameters of the harvester. computation of linear displacement, linear velocity, angular displacement, and angular velocity of the system at each time step. The solution from the preceding iterative calculation serves as the initial condition for the subsequent step, facilitating mesh updates.