# Coupling Motion and Energy Harvesting of Two Side-by-Side Flexible Plates in a 3D Uniform Flow
## 2 Physical Problem and Numerical Methods
### Numerical Methods
#### 2.2.1 Mathematical Formulation  
Due to its efficiency and the relative simplicity of the computation, the lattice Boltzmann equation (LBE) with the BGK model has been widely employed to simulate a complex flow field. This method acts as an alternative to the conventional numerical methods in solving Navier-Stokes equations. The BGK-based LBM can be expressed as [35; 36; 41]:  
\[f_{i}\left(\mathbf{x}+c_{i}\Delta t,t+\Delta t\right)-f_{i}\left(\mathbf{x}, t\right)=-\frac{1}{\tau}\left(f_{i}\left(\mathbf{x},t\right)-f_{i}^{eq}\left( \mathbf{x},t\right)\right)+f_{i}^{F}\left(\mathbf{x},t\right)\Delta t \tag{7}\]  
where \(\tau\) is the dimensionless relaxation time, \(\tau=\mu/c_{s}^{2}+1/2\Delta t\); \(\Delta t\) is the time increment and \(t\) is the lattice time step. Additionally, \(f_{i}\left(\mathbf{x},t\right)\) is the distributed function for particles at position \(\mathbf{x}\) with velocity \(\mathbf{c}_{i}\). The lattice equilibrium distributed function \(f_{i}^{eq}\) and the forcing function \(f_{i}^{F}\) are given as [42]:  
\[f_{i}^{eq}=\omega_{l}\rho\left(1+\frac{\mathbf{c}_{i}\cdot\mathbf{u}}{c_{s}^{ 2}}+\frac{\left(\mathbf{c}_{i}\mathbf{c}_{i}-c_{s}^{2}\mathbf{I}\right): \mathbf{u}\mathbf{u}}{2c_{s}^{4}}\right) \tag{8}\]  
\[f_{i}^{F}=\left(1-\frac{1}{2\tau}\right)\omega_{i}\left(\frac{\mathbf{c}_{i}- \mathbf{u}}{c_{s}^{2}}+\frac{\mathbf{c}_{i}\cdot\mathbf{u}}{c_{s}^{4}}\mathbf{ c}_{i}\right)\cdot\mathbf{f}_{s} \tag{9}\]  
where \(\omega_{i}\) is the weighted factor, in the present \(D3Q19\) model [36; 43], \(\omega_{0}=1/3\), \(\omega_{1-6}=1/18\), \(\omega_{7-18}=1/36\). Additionally, \(c_{s}\) is the lattice speed of sound. The mass density \(\rho\) and fluid velocity \(\mathbf{u}\) can be obtained by:  
\[\rho=\sum_{i}f_{i} \tag{10}\]  
\[\rho\mathbf{u}=\sum_{i}f_{i}\mathbf{c}_{i}+\frac{\Delta t}{2}\mathbf{f}_{s} \tag{11}\]  
With the Dirac delta function in the immersed boundary model, the Lagrangian momentum forcing term is spread to the Eulerian grids, while the Eulerian velocities are interpolated to the Lagrangian meshes; the numerical schemes can be expressed as [36; 37]:  
\[\mathbf{f}_{s}(\mathbf{x},t)=\sum_{s_{1},s_{2}}\mathbf{F}\left(s_{1},s_{2},t \right)\delta_{\Delta}\left(\mathbf{x}-\mathbf{X}\left(s_{1},s_{2},t\right) \right)\Delta s_{1}\Delta s_{2} \tag{12}\]  
\[\mathbf{U}\left(s_{1},s_{2},t+\Delta t\right)=\sum_{x}\mathbf{u}\left( \mathbf{x},t+\Delta t\right)\delta_{\Delta}\left(\mathbf{x}-\mathbf{X}\left( s_{1},s_{2},t\right)\right)\Delta x^{3} \tag{13}\]  
Based on our careful validations and convergence examinations with different computational sizes and lattice spacings, the computational domain is modeled as a rectangular box, extending from \([0,0,0]\) to \(\left[L_{x},L_{y},L_{z}\right]\left(L_{x}=8,L_{y}=4,L_{z}=12\right)\) in the \(X\) (streamwise), \(Y\) (spanwise) and \(Z\) (transverse) directions, respectively. The grid spacing has a uniform value of \(\Delta x=\Delta y=\Delta z=0.02\). Here, both the domain size and grid spacing are scaled by the plate's length \(L\). To simulate the plates' flap in free stream, the Dirichlet boundary condition \(\left(u=U_{\infty},v=0,w=0\right)\) is adopted at the inflow (\(x=0\)), and a convective boundary condition \(\left(\partial\mathbf{u}/\partial t+U_{\infty}\partial\mathbf{u}/\partial x=0\right)\) is applied at the outflow (\(x=L_{x}\)). Meanwhile, the far-field boundary condition and periodic boundary condition are used in the transverse (\(z=0\) and \(z=L_{z}\)) and spanwise (\(y=0\) and \(y=L_{y}\)) direction, respectively. The fixed boundaries of the plates (\(x_{b}=0\)) are aligned with the \(Y\)-axis, with their mid-points \(C_{1}\) and \(C_{2}\) locatedin the mid-plane of the computational domain in the Eulerian coordinate system (\(x=L\),\(y=1/2L_{y}\)). Each of the plates is discretized using a three-node triangular element [38].