In wind tunnel tests, pressures are usually measured at each pressure tap on a prototype scale model of a solar panel array with wind azimuth direction and tilt angle of modules. At each pressure tap, the time history of the pressure coefficient C p ( t ) can be obtained as Equation 1 (Abiola-Ogedengbe et al., 2015; Letchford et al., 1993; Holmes, 2007; Simiu and Yeo, 2019): C p t = p t − p ∞ 0.5 ρ V ∞ 2 , (1) where p ( t ) is the instantaneous pressure at the pressure tap, p ∞ is the atmospheric pressure, ρ is the air density, and V ∞ is the mean wind velocity. Although wind tunnel tests provide wind-induced pressure measurements, these pressures are obtained from prototype scale solar panel array models with substantial experimental resources. While mean wind loads are generally not significantly affected by geometric scales, peak wind loads are sensitive to geometric scales (Aly and Bitsuamlak, 2013). This may result, in certain cases, in unsatisfactory estimation of wind load time histories and dynamic responses of the solar panel array. Further, wind tunnel tests typically do not account for the spatial correlation along the horizontal length of the solar panel array, which would affect the dynamic response of the solar modules, and in particular the excitation of certain modes.

In this regard, a novel simulation framework is proposed to overcome the limitation of wind loads generation in wind tunnel tests that is affected by the geometric do not account for the spatial correlation. Specifically, instead of employing constant mean wind velocities in Equation 1, stochastic wind velocities that vary in one spatial dimension are modeled as a stochastic wave that is continuous in both space and time (Benowitz and Deodatis, 2015; Chen et al., 2018) as follows. The frequency-wavenumber spectrum can be expressed in Equation 2 as S ω , κ = 1 2 π ∫ − ∞ ∞ S ω ⋅ γ ξ , ω e i κ ξ d ξ , (2) where S ( ω , κ ) is the power spectrum, γ ( ξ , ω ) is the coherence function, ω is the frequency, ξ is the spatial distance, and κ is the wavenumber. This paper employs the Kaimal two-sided spectrum (Benowitz and Deodatis, 2015; Chen et al., 2018; Zhao et al., 2023) for wind velocity simulation, defined as: S k ω = 1 2 ⋅ 200 2 π ⋅ u ∗ 2 ⋅ z U z ⋅ 1 1 + 50 ω z 2 π U z 5 3 , (3) where z is the height above the ground, U ( z ) is the average wind velocity at height z , and u ∗ is the velocity friction component defined as u ∗ = k U z ln z z 0 . (4)

In Equation 4, k represents von Karman’s constant and k = 0.4 in this paper, z 0 represents ground roughness, specified as 0.01 m since the solar panel site at MI RTC is open terrain. Similarly, the Davenport coherence function, which is popular in wind engineering, is selected and defined as γ ξ , ω = e − λ ω ξ / 2 π U z , (5)

where λ is a decay parameter, often chosen to be between 7 and 10 (Benowitz and Deodatis, 2015; Chen et al., 2018). In this paper, λ = 10 is chosen. Further, using Equations 3, 5, the frequency-wavenumber spectrum S ( ω , ξ ) can be expressed in Equation 6 as S ω , κ = 1 2 ⋅ 200 2 π ⋅ u ∗ 2 ⋅ z U z ⋅ 1 1 + 50 ω z 2 π U z 5 3 ⋅ λ ω π U z κ 2 + λ ω 2 π U z 2 . (6)

The Spectral Representation Method (SRM) can then be used to simulate sample realizations of a stochastic wave from its frequency-wavenumber spectrum (Benowitz and Deodatis, 2015; Peng et al., 2016; Peng et al., 2017; Chen et al., 2018; Zhao et al., 2023) as follows: u x , t = 2 ∑ i = 1 N κ ∑ j = 1 N ω 2 ⋅ S ω j , κ i Δ κ Δ ω ⋅ cos κ i x + ω j t + φ i j + cos κ i x − ω j t + φ ̃ i j , (7)

where N κ is the wave number; Δ κ = κ u / N κ , where κ u is the cut-off wave number; N ω is the number of discretized frequency; Δ ω = ω u / N ω , where ω u is the cut-off frequency; φ i j and φ ̃ i j are two independent sets of random phase angles uniformly distributed over [ 0,2 π ] . In this paper, the following values are selected: N κ = 8192 , Δ κ = 0.004 m − 1 , κ u = 32.99 m − 1 , N ω = 1048 , Δ ω = 0.024 rad/s, ω u = 8 π rad/s.

Further, based on the frequency and wave number values, the time duration is T = 2 π / Δ ω = 261.75 s, the number of time step is N t = 2 ∗ N ω = 2096 , and the time step is Δ t = T / N t = 0.1249 s.

Further, the horizontal domain is discretized for spatial correlation consideration. In this paper, the number of space N x = 2 ∗ N k = 16384 , and Δ x = L / N x = 0.002 m , where L is the horizontal length of the solar panel array. Note that only the grid points located on the module surface are retained for stochastic wind loads estimation.

Equation 7 requires substantial computational resources. However, the Fast Fourier Transform (FFT) can be employed to dramatically increase the computational efficiency of simulation (Benowitz and Deodatis, 2015; Zhao et al., 2023). To achieve this, sample realizations produced by Equation 7 can instead be expressed as u x , t = Re [   ∑ l = 0 N κ − 1 ∑ m = 0 N ω − 1 B l m 1 ⁡ exp i ω m t + i κ l x + B l m 2 ⁡ exp − i ω m t + i κ l x ] . (8)

In Equation 8, Re [ ⋅ ] is the real part of the realizations, and B l m n = 2 S ω m , κ l ⋅ Δ ω Δ κ e i ϕ l m n . (9)

Equation 8 can then be re-written as u x , t = Re FFT κ FFT ω B 1 + FFT κ IFFT ω B 2 . (10) where FFT ( ⋅ ) and IFFT ( ⋅ ) are the Fast Fourier Transform and Inverse Fast Fourier Transform, respectively. Further, the total wind velocities simulated on a solar panel array can be expressed as u ̃ x , t = U z + u x , t . (11)

In this section, the constant mean wind velocity at the reference height of the solar panel array is assumed to be 9 m/s at the reference height of horizontal center-line of solar panels, and the wind velocity profile follows the wind power law, as it was in Section 2.1. Sample realizations of wind velocities in both space and time are shown in Figure 1. Further, the accuracy of the simulated wind velocity can be validated by comparing the power spectrum of the simulated wind velocity to the Kaimal spectrum. The power spectrum of the simulated u ( x , t ) can be estimated from (Benowitz and Deodatis, 2015) S j j ω = | ∫ 0 T u x j , t e − i ω t d t | 2 2 π T , (12) where u ( x j , t ) denotes the simulated wind velocities at location x j = 2 m, and T denotes the duration of time. Figure 2 shows that the estimated power spectrum from Equation 12 matches the desired Kaimal spectrum from Equation 3. Additionally, the accuracy of the simulated wind can be validated by comparing the estimated and target cross spectra. For sample realizations u ( x j , t ) and u ( x k , t ) , cross spectra cross spectra can be obtained from Equation 13 shown as S j k ω = ∫ 0 T u x j , t e − i ω t d t ∫ 0 T u x k , t e − i ω t d t 2 π T , (13)

and the coherence functions are estimated as γ ξ j k , ω = S j k ω S j j ω S k k ω . (14)

In this section, three locations ( x 1 = 3 m, x 2 = 8 m, and x 3 = 17 m) are used to assess the estimated cross spectra. Therefore, the separation distances are ξ 12 = 5 m, ξ 13 = 14 m, ξ 23 = 9 m. Estimated coherence functions in Equation 14 and target coherence functions in Equation 5 are shown in Figure 3. As it can be seen, the agreement between estimated and target coherence functions can be satisfied. Note that the vertical spatial correlation is not considered because the module height is only 1.012 m. In this manner, variations in vertical coherence are expected to be limited compared with the along-row correlation in this study.

Next, a refinement of dynamic wind pressure estimation from Equation 1 is proposed to further increase computation efficiency. Specifically, the mean pressure coefficient distributions on the solar panel array can be obtained from pressure coefficient time histories. However, this still needs experimental resources to conduct wind tunnel tests. In this regard, an alternative computational fluid dynamics (CFD) simulation is used to generate mean pressure coefficient distributions on solar panel arrays. Specifically, the CFD is performed large-eddy simulation (LES) with standard Smagorinsky Sub-grid Scale (SGS) model (Smagorinsky, 1963). LES can adopt the inflow perfectly (Zhu et al., 2020), and the standard Smagorinsky SGS model is adopted due to its straightforward implementation and relative low computational cos. In addition, prior studies have reported that the mean quantities are not sensitive of the choice among Smagorinsky, dynamic Smagorinsky, or WALE SGS models (Cant and Lee, 2017).

Assuming the wind flow around the solar panel arrays as an incompressible fluid, the filtered Navier–Stokes equations that govern the motion of the large-scale eddies can be written in Equations 15, 16 as (Melaku et al., 2022; Melaku and Bitsuamlak, 2024) ∂ u ̃ i ∂ x i = 0 , (15)

and ∂ u ̃ i ∂ t + ∂ ∂ x j u ̃ i u ̃ j = − 1 ρ ∂ p ̃ ∂ x i + ν ∂ 2 u ̃ i ∂ x j ∂ x j − ∂ τ i j ∂ x j , (16)

where u ̃ i is the i th component of the filtered velocity, p ̃ is the filtered pressure, ν is the kinematic viscosity, and ρ is the air density of the wind. In this section, the properties of air ρ = 1.225   k g / m 3 and ν = 1.5 × 1 0 − 5   m 2 / s are applied. Further, τ i j is the shear stress on small-scale eddies and is expressed as (Melaku et al., 2022; Melaku and Bitsuamlak, 2024) τ i j = u i u j ̃ − u ̃ i u ̃ j . (17)

The shear stress in Equation 17 is expressed as subgrid-scale (SGS) stress and is often modeled using the eddy-viscosity hypothesis. In this study, the standard Smagorinsky SGS model can be adopted in Equations 18, 19 as τ i j − 1 3 δ i j τ k k = − 2 ν t S ̃ i j , (18)

with S ̃ i j = 1 2 ∂ u ̃ i ∂ x j + ∂ u ̃ j ∂ x i , (19)

where δ i j is the Kronecker delta, S ̃ i j is the strain rate tensor of the resolved velocities, and ν t is the turbulent eddy-viscosity that can be expressed in Equation 20 as (Melaku et al., 2022; Melaku and Bitsuamlak, 2024) ν t = C s Δ 2 | S ̃ | . (20)

In this equation, Δ is the filter width defined as the cube-root of the cell volume, and the coefficient C s represents Smagorinsky’s constant and can be expressed in Equation 21 as (Smagorinsky, 1963) C s = C k C k C ε 1 / 2 , (21)

where the coefficient C k and C ε are model constants. C k and C ε with typical values of 0.094 and 1.048 are applied in this study, respectively, yielding Smagorinsky’s constant C s = 0.1678 .

Note that the mean pressure coefficients from only the smooth flow atmospheric boundary layer (ABL) conditions need to be obtained. There are two reasons to propose this novel idea. First, it is observed that over a sufficiently long time period, the mean pressure coefficient distributions under different turbulence intensities exhibit similar patterns (Aly and Bitsuamlak, 2013). Second, substantial computational resources and time are needed to use CFD simulation to generate pressure coefficient time histories of turbulence inflow, and the duration of the pressure coefficient time histories is usually not sufficient to obtain reliable mean pressure coefficient values. However, CFD simulation based on smooth inflow does not require long time period to derive useful mean pressure coefficient values, which also reduces required computational resources.

To apply LES modeling with smooth inflow ABL conditions on a solar panel array model, a CFD simulation software called the Wind Engineering with Uncertainty Quantification (WE-UQ) is applied to estimate the mean pressure coefficient distributions.

WE-UQ has a user-friendly Graphical User Interface (GUI), which provides basic procedures and modeling the virtual wind tunnel and run CFD simulation (Melaku and Bitsuamlak, 2024). In this paper, the basic geometries of a 1:10 scale of the solar panel array model at Michigan Solar Regional Test Center (MI RTC), can be modeled and imported to WE-UQ for visualization, and the domain size of the wind tunnel can also be defined and visualized. Note that, although mean pressure coefficients are not significantly affected by geometry scale (Aly and Bitsuamlak, 2013), changing the scale may affect the Reynold number, which may still influence mean pressure coefficient values to some extent. Then, WE-UQ can define and visualize the domain size of the wind tunnel and establish the mesh, including boundary mesh on the virtual wind tunnel and mesh surrounding the imported model in the wind tunnel, as shown on the right side of Figure 4. Next, the boundary conditions of the wind tunnel can be configured and shown in Figure 5a to replicate the environmental conditions observed at the MI RTC solar panel array site, ensuring consistency between the CFD simulation and the actual testing environment. The constant mean wind velocity is 9 m/s at the reference height, which is the horizontal center-line of the solar panel array; the ground roughness is 0.01 m, which can make sure the power law constant mean wind profile imported to WE-UQ is identical with that from the stationary wind speed estimation from Section 2.1 and the identical with mean wind speed condition from wind tunnel tests in Section 2.1. Air density ρ = 1.225   k g / m 3 and kinematic viscosity of the air ν = 1.5 × 1 0 − 5   m 2 / s are applied in this simulation. Then, in numerical setup section from WE-UQ in Figure 5b, large-eddy simulation (LES) with Smagorinsky Sub-grid Scale (SGS) model with coefficients is selected for numerical simulation. Next, the CFD simulations are performed using the pimpleFoam solver. This solver is selected for its ability to maintain numerical stability and ensure accurate wind load predictions while allowing for relatively larger time steps in transient simulations. Since induced wind is specified as smooth inflow, a small time step of 0.0001 s is adopted to capture transient effects with high resolution. Next, several maximum Courant numbers C o are tested for computing mean pressure coefficient distributions, as shown in Figure 6. As it can be seen, the mean pressure coefficient distributions are similar for all selected C o values. Because C o = 1 is more computationally expensive than C o = 3 and C o = 5 , C o = 5 is chosen for efficient computation, further control of solver convergence, and numerical stability. Since smooth inflow has higher numerical stability in a transient simulation, pimpleFoam will exhibit fast convergence behavior. Then, the pressure coefficient time histories from smooth inflow is obtained, and one sample is shown in Figure 7. After an initial transient of 0.3 s, the time history becomes stationary. In this manner, the simulation duration is limited to 3 s, which provides an adequate averaging window to obtain mean pressure coefficient distributions.

Further, CFD simulations under a smooth flow ABL condition are conducted to examine whether changes in wind speed influence the mean pressure coefficients. Specifically, constant mean wind velocities of 9 m/s and 38 m/s at the reference height (the horizontal centerline of the solar panel array) are considered. Using the WE-UQ software for CFD simulation, the distributions of mean pressure coefficients from 38 m/s inflow wind speed on the solar panel model are obtained and presented in Figure 8a. Additionally, the ratio of mean pressure coefficient distributions between mean wind speeds 9 m/s and 38 m/s is calculated and presented in Figure 8b. This is done to assess whether the difference is significant. As it can be seen from Figures 8a,b, the mean pressure coefficient distributions along the solar panel arrays remain close under different values of inflow wind speed. This indicates that in CFD simulation for smooth flow ABL conditions, the change of inflow mean wind speed does not change the mean pressure coefficients significantly, and the slight difference of the pressure coefficients between two inflow speeds attribute the uncertainties generated by the mesh of the solar panel array in WE-UQ.

Additionally, as mentioned in Section 1, peak pressure coefficients are significantly influenced by the geometric scale of the solar panel array, while mean pressure coefficients are not significantly influenced by model size. Thus, the use of mean pressure coefficients rather than time histories reduces potential scaling errors. Using mean pressure coefficients, the dynamic wind pressure can be computed as p x , t = 0.5 ρ u x , t 2 C ̄ p + ρ ∞ , (22) where u ( x , t ) is again the simulated stationary stochastic wind velocity and C ̄ p is the mean pressure coefficient, both at a given location. In this context, a simple scheme can be built to estimate stationary wind loads based on spatially correlated stochastic wind speed and mean pressure coefficient distributions on the solar panel model:

Select the appropriate power spectrum of stationary wind turbulence. Kaimal spectrum is selected in this paper.

Further, the torque tube in the solar panel array is used to rotate the modules with respect to the position of the sun. So it must be designed to withstand the torque induced by twisting during heavy winds. Additionally, the distribution of twist along the torque tube gives an indication of the global response of the structure. Monitoring the torsional movement of the torque tube under wind loads is useful to ensure that the torque tube can prevent system failures and the dynamic response can be further understood.

In this regard, by applying the proposed simulation framework in section 2, 9 m/s and 38 m/s mean wind speed at the same reference height are used to estimate the torsional strain of the torque tube under different wind conditions, and the results are shown in Figure 17. As it can be seen, compared to the proposed framework based on 9 m/s mean wind speed, 38 m/s mean wind speed generates much larger torsional strain time series. Then, the maximum strains in Figure 17c from 9 m/s and 38 m/s mean wind speed are used respectively to calculate the twisting angle of the torque tube at location T3 shown in Figure 10. At 9 m/s, the twist angle is 0.04 ° ; at 38 m/s, the twist angle is 0.84 ° . Then, the peak differential displacements between the two ends of the upper frame of the right end solar panel module which is attached to T3 are derived respectively. At 9 m/s, the peak displacement is 1.3 mm; at 38 m/s, the peak displacement is 29 mm. This indicates that the deformation of the solar panel module under strong wind conditions may induce permanent damage to the module and detach the module from the torque tube.

Further, to account for the effect of spatial variation on dynamic response of the torque tube, the peak differential settlement of the module under wind velocity without spatial variation is derived. Specifically, the fundamental spectral representation method (SRM) is used to generate stationary wind velocity time histories that are applied uniformly across the horizontal domain (Abrous and Wamkeue, 2015). In this case, the peak differential displacement is only 0.4 mm, indicating that neglecting spatial variation can substantially underestimate the dynamic response of the solar panel array.

Next, the dominant frequency components contributing to the highest strain in the time series data are evaluated. This is done by applying the Inverse Fast Fourier Transform (IFFT) to reconstruct the signals corresponding to certain frequencies, and applying frequency domain decomposition of the strain time series. Results for an example strain time history simulated at location T3 are shown in Figure 18. Figure 18 shows that more excitation of higher modes from 38 m/s mean wind speed is again demonstrated.

Further, the dynamic characteristics of the torsional strains are also needed for design purposes. In this manner, it is appropriate to observe the time-varying frequency content of the response of the torsional strain under a representative stochastic wind load sample from each case described in Section 2. Since the torsional strains at location T3 are greater than at other positions on the torque tube, this makes the characteristics of torsional strains more distinguishable for analysis. Therefore, the spectrograms of torsional strain at T3 in Figure 10, under the two cases of induced wind loading, are selected and shown in Figure 19. In Figure 19, the dashed-dot black lines represent standard deviation of frequency values. As it can be observed from Figures 19a,b, the time-frequency content from 38 m/s mean wins speed is generally much more powerful than the content from 9 m/s mean wind speed, which causes the modes to be more excited. Further, Figure 19c shows that with the high mean wind speed, the torsional strain is much more powerful than the torsional strain in Figure 19b, and modes can be much more excited and visualized in spectrogram.

Next, from Section 2, 200 samples of stationary wind velocities based on 9 m/s mean wind speed and 38 m/s mean wind speed are used to simulate corresponding samples of torsional strain time series, respectively. Then, the ensemble average of the spectrogram values of the strains is shown in Figure 20. As it can be seen, the dominant frequencies at which the structural modes are excited still remain obvious in both cases, which can validate the accuracy of the proposed framework to identify the mode excitation. Figure 20 also shows that the time-frequency content from 38 m/s mean wins speed is still generally much more powerful than the content from 9 m/s mean wind speed.