# Parallel Calculation of Random Vibration for Reentry Vehicles in Fluctuating Pressure Environment

Authors: Xuanhua Fan; Jun Liang; Keying Wang; Huan Zhou
Journal: Shock and Vibration
 (2021)

---

## Abstract

The reentry vehicle is in an unconstrained free-flight state during the reentry process. The pressure of air fluctuations acts on the external surface of the vehicle and induces a random vibration environment with a wide-frequency range, which makes great influences on the structural dynamic properties of vehicles. This paper aims at the parallel calculation of the wide-frequency random vibration problems for the free body under multipoint fluctuating pressure loads. The calculation method is based on the self-developed JAUMIN framework and PANDA platform. The novel algorithm based on the modal superposition method reduces the computational complexity, realizes the parallel calculation, and achieves the maximum numerical simulation calculation capacity of 1.045 billion degrees of freedom (DOFs) while improving the calculation efficiency and the accuracy of the results. Finally, this article also uses typical examples to verify the correctness and parallel scalability of algorithms and programs.


---

## Body

## 1. Introduction

Reentry vehicles suffer from high levels of reentry fluctuating pressure caused by high-speed air currents in the dense atmosphere. The random vibration of the vehicles caused by these violent fluctuations may eventually lead to the vibration of the internal components and even the fatigue or failure of the structure [1, 2]. This fact forces vehicle producers to constantly look for improvements in the areas of quality and reliability. Life span, flawless operation, and safety are directly interconnected [3, 4].The mechanism of fluctuating pressure is very complicated. High Mach, high Reynolds, laminar-turbulent transition, turbulent vortex shedding, and some other typical airflows may generate high-level fluctuating pressure to the vehicles [5, 6]. A series of research works has been implemented to predict fluctuating pressure, such as theoretical analysis [7, 8], CFD (computational fluid dynamics) simulation [9, 10], wind tunnel test [11, 12], and flight tests [13, 14], and so on. In this paper, the assumed known fluctuating pressures are taken as the inputs for the random vibration analysis.Academically, the structural dynamics analysis of the reentry vehicle subjected to multi-point fluctuating pressure can be summarized as the random vibration problem of the free body including rigid modes. The solution to such problems is challenging, mainly due to the following reasons:(1)
The dynamic modes of the vehicles reentering the Earth’s atmosphere in free-flying status, including rigid-body motion and the elastic deformation. Because the stiffness matrix of finite element analysis for the free body is singular, modal analysis becomes difficult [15].(2)
The complex distribution and coupling characteristics of fluctuating pressure on the external surface of reentry vehicles make it difficult to describe the fluctuating pressure [2, 12]. The value of the fluctuating pressure varies with the position of the outer surface, exhibiting nonconforming distribution properties. Moreover, the fluctuating pressure in adjacent areas is strongly coupled. In the viewpoint of frequency spectra, the fluctuating pressure is a wide-frequency load with a range of more than several kilohertz [12, 16], posing another challenge for dynamic finite element analysis.(3)
Random vibration algorithms for the free body are different from those for constrained structures. Traditional commercial finite element codes mainly deal with random vibration analysis of constrained structures. Random vibration analysis of the free body should take the rigid-flexible coupling effect into consideration, including not only the effect of rigid body motion on elastic deformations but also the influence of elastic deformations on rigid body motion [17, 18]. Furthermore, both the auto-power spectral density (auto-PSD) and the cross-power spectral density (cross-PSD) of fluctuating pressure should be considered to describe the correlation of loads.In order to improve the analysis precision of random vibration in the mid- or high-frequency range, some local dynamic characteristics of the reentry vehicles should be described correctly. A refined model is usually required in finite element analysis. The corresponding calculation scale is as high as 100 million DOFs or above [19], which is beyond the capabilities of most commercial codes. Consequently, a large-scale parallel calculation has become particularly important. In other words, the large-scale refined finite element analysis is an effective way to increase the upper limit of random vibration analysis [20]. For large-scale random vibration analysis problems, the modal superposition method [21] is usually a feasible and effective way to obtain the dynamic vibration response of structures. After evaluating a set of orthogonal modal shape vectors, this method reduces a large number of global equilibrium equations to a relatively small number of uncoupled second-order differential equations through coordinate transformation and equation decoupling. The response under each modal coordinate is calculated, and then, these responses are combined to obtain an approximate solution of the structural dynamic response under the original physical coordinate [22].In the parallel realization of large-scale structural dynamics analysis, the design of the underlying framework is particularly important. At present, some well-known underlying frameworks, such as SIERRA [23], PETSc [24], PHG [25], and so on, have been developed or are under development. The basic design concept of these frameworks is to provide common data structures, parallel strategies, and interfaces with various application programs, that is, to develop massively parallel programs using the method of “integrating common features of some physical or mathematical problems to support frameworks and developing parallel application programs based on the framework.” The framework has been proven to be an effective way to develop massively parallel programs for various supercomputers [26, 27]. JAUMIN [26] (J adaptive unstructured mesh application infrastructure), which has been developed by the China Academy of Engineering Physics since 2011, is a basic parallel framework for unstructured mesh finite element applications. The JAUMIN framework supports the large-scale parallel simulations on adaptive unstructured meshes using massively parallel processing (MPP) machines. It is a complement to the previously developed JASMIN [28] framework, which aims at structured mesh applications in fluid mechanics or other related fields. Based on the JAUMIN framework, an application platform named PANDA (parallel adaptive nonlinear dynamic and static analysis) is being developed by our team. PANDA is a large-scale calculation platform for finite element analysis of engineering structures.This paper introduces the research details on the random vibration simulations of the reentry vehicles, including the algorithm design and parallel implementation of wide-frequency random vibration for the free body under multi-point fluctuating pressure. A corresponding large-scale parallel solution system is constructed. The rest of this paper is organized as follows: Section2 introduces the theory and algorithm for random vibration of the free body in detail, including the modal analysis of the free body, the description of fluctuating pressure, and the kernel theoretical derivation of free-body random vibration; Section 3 introduces the parallel implementation of random vibration in virtue of the JAUMIN framework and the PANDA platform; some representative examples are given to verify the correctness and parallel scalability of the algorithm and application in Section 4; and finally, this paper makes a brief conclusion in Section 5.


## 2. Theories and Algorithms

For vehicles, the mode superposition method is a preferred one for large-scale random vibration analysis in a wide-frequency range. The whole solution procedure includes modal analysis of the free body, description of fluctuating pressure, as well as the solution process and results output of structural responses. In this section, the above aspects will be illustrated.


### 2.1. Modal Analysis of Free Body

In virtue of finite element discretization and matrix assembly, the generalized eigenvalue equation corresponding to the modal analysis of free-body structures can be written as follows:(1)Kx=ω2Mx,where M and K are the mass matrix and the stiffness one, respectively; ω is the natural circular frequency; and x is the modal shape vector. Due to the totally free-flying status for vehicles in the entry atmosphere, the modal shapes of which include six rigid motion modes with ω = 0, that is, three translational ones and three rotational ones along the coordinate direction of X, Y, and Z. For free-body modal analysis, the key is to effectively solve the problem of eigenvalues caused by the singularity of the stiffness matrix and to effectively separate the rigid motion modes and elastic modes from the free-body modes.The existing mainstream algorithms for solving large-scale modal analysis are mainly subspace iterative methods [27, 29], such as the Krylov–Schur method [30] and the Jacobi–Davidson method [31]. However, these methods are mainly suitable for nonsingular matrix situations; for a free-body structure, the stiffness matrix K is a singular matrix. For this reason, the stiffness matrix needs to be processed to eliminate the singularity. The method we adopted in this paper is to introduce spectral transformation technology [32]. By introducing a positive shift value of σ into equation (1), we have(2)K+σMx=ω2+σMx.Taking (K + σM) as a matrix operator and substituting (K + σM) for K in modal analysis, the singularity of equation (1) can be eliminated, and equation (2) can be solved by subspace iterative methods. Comparing equation (1) with equation (2), it is found that whose corresponding modal shape vectors are unaltered, and the only eigenvalue difference between which is σ. Thus, the eigenpairs of equation (1) can be easily obtained by a transformation when equation (2) is computed by subspace iterative methods. In this way, the singularity problem of the free-body stiffness matrix is effectively solved.After above solution, we can getmt lowest-order eigenpairs of [(ωi, ϕi), i = 1,2, …,mt]. These eigenpairs will be used for the mode superposition of random vibration. Among these lowest eigenpairs, the first six ones are rigid body modes with ωi = 0, and the rest are elastic modes with ωi > 0. For the convenience of subsequent theoretical derivation, these eigenpairs are divided into two parts, and all output modal shape vectors are orthonormalized with mass matrix M, being satisfied with the following condition:(3)ϕiTMϕj=1i=j,0i≠j.Taking into account that the free-body structure may have a partially restricted degree of freedom, we usem1 and m2 to represent the number of rigid body modes and elastic modes, respectively, where m1 is equal to 6 for reentry vehicles, and for other noncompletely free cases, m1 is less than 6. Denoting the rigid modal shape as ϕRi and the elastic modal shape as ϕEi, we can get(4)ϕRi=ϕi,fori=1,2,…,m1,ϕEi=ϕi+m1,fori=1,2,…,m2.


### 2.2. Description of Fluctuating Pressures

Generally, the reentry process of a vehicle is treated as a stationary random process in the finite element analysis of random vibration [33]. The basic inputs for random vibration are the auto-PSDs and cross-PSDs of fluctuating pressures. According to these inputs, the auto-PSDs of displacement, the velocity, and acceleration for concerned points inside the structures and the equivalent root-mean-square (RMS) responses of full field are calculated via random vibration analysis.It is obvious that the fluctuating pressure distributions on the external surface of reentry vehicles are extremely complex. Consequently, it is impossible to describe the fluctuating pressure at each discrete finite element node of the outside surface, especially for the case of large-scale computations.In this paper, a treatment of geometry area decomposition on the outside surface is adopted to approximate the pressure distributions. The external surface of the vehicle is divided into numerous areas, and a uniform pressure value is specified in each area. In order to represent the real fluctuating pressure distributions more exactly, the geometry areas can be further subdivided to approximate the real loads.Without loss of generality, we assume that there arer different fluctuating pressures acting on the external surface, being denoted as p1,p2,…,pr. The acting areas of these pressures constitute the whole load region of the reentry vehicle. It should be noted that these fluctuating pressures cannot be used directly in the implementation of random vibration. Instead of that, the auto-PSDs Spipi (i = 1,2, …, r) and the cross-PSDs Spipj (i = 1,2, …, r; j = 1,2, …, r, i ≠ j) of these fluctuating pressures are the direct input data for random vibration of free body. For the stationary random process, the relationship between fluctuating pressures and their PSDs can be shown as follows [34]:(5)Splpmω=limT⟶∞12Tplωpm∗ω.Theoretically, the fluctuating pressures on two different geometry areas might be totally correlated, partly correlated, or uncorrelated. Taking each case into consideration, the input PSDs are designed as anr × r matrix as follows:(6)Sp1p1Sp1p2⋯Sp1prSp2p1Sp2p2⋯Sp2pr⋮⋮⋱⋮Sprp1Sprp2⋯Sprpr,where the diagonal elements represent the auto-PSDs and the nondiagonal elements represent the cross-PSDs. The matrix elements in equation (6) are nonzero in the condition that the fluctuating pressures between these two areas are correlated, otherwise zero.In this way, the fluctuating pressure loads on the outer surface of the reentry vehicle are treated as anr-dimensional spectral curve matrix; each element in the matrix corresponds to a spectral curve or 0. The advantage of this processing is that the matrix can describe various situations such as completely related, partially related, and completely unrelated.Although the fluctuating pressures are not the direct input for the random vibration, their distributing information are still necessary for the assembly of dynamic equations and the computation of modal participation factors. The discrete finite-element dynamic equation for the reentry vehicle with rigid motion can be described as follows:(7)Mu¨+Cu˙+Ku=f,where C is the damping matrix, u is the displacement vector for all degrees of freedom, and f is the fluctuating pressure vector. In the construction of equation (7), f is assembled via the integration of each pressure on its corresponding element surface and then combining them together. f can be described as follows:(8)f=∑l=1rflpl,where fl is the coefficient vector of discrete load with unit fluctuating pressure on the acting area of pl.In order to determinefl, a unit pressure load in this paper is supposed to apply to the acting area of pl, and finite element integration and discretization are carried out to form fl. It can be seen in the following theoretical derivations and implementations that neither the fluctuating pressure pl nor the total discretized load f needs to be explicitly given or calculated. The actual input for random vibration of the free body is the load factor vector fl and the PSDs of fluctuating pressures.


### 2.3. Theories and Algorithm for Random Vibration of the Free Body

#### 2.3.1. Modal Decoupling and Separation of Rigid Displacements and Elastic Ones

The displacement vectoru in equation (7) can be divided into two parts: the rigid displacement vector uR and the elastic displacement vector uE [35]. When the modal superposition method is adopted, the two parts of displacement can be described by rigid modal shapes and elastic ones:(9)u=uR+uE=ΦRηR+ΦEηE,where ΦR and ΦE are the rigid modal matrix and the elastic one, respectively, and ηR and ηE are the rigid modal coordinate vector and the elastic one, respectively.It should be noted that equation (9) is approximately established due to the modal truncation in the modal superposition method [36]. Usually, the more the modal shapes we used, the better the precision we achieve. Since K is singular (KuR = 0), we achieve the following equation by substituting equation (9) into equation (7):(10)Mu¨E+Cu˙E+KuE=f−Mu¨R−Cu˙R.According to the orthogonality (equation (3)) and the small damping assumption, we can get two sets of decoupled equations by substituting uE=ΦEηE and uR=ΦRηR into equation (10) and premultiplying both sides of equation (10) by ΦΕΤ and ΦRΤ, as follows:(11)MRη¨R+CRη˙R=ΦRΤf,(12)MEη¨E+CEη˙E+KEηE=ΦΕΤf,where MR=ΦRTMΦR=IR, ME=ΦΕΤMΦE=IE, CR=ΦRTCΦR, CE=ΦETCΦE, and ΚE=ΦETΚΦE.Equations (11) and (12) give the rigid modal coordinate description [37] and the elastic modal coordinate description of equation (7), respectively. In other words, the dynamic equation of the free body is separated effectively with the rigid part and the elastic one. Comparing equation (11) with equation (12), we can find that both equations are nearly the same except for a missing singular stiffness term in equation (11). In the following, the relatively complex case of the elastic response of the structure is deduced first, and then, the rigid response and the rigid-elastic coupling response are deduced.It should be noted that the following derivation process of elastic displacement response, rigid body displacement response, and the total response of the coupling of the two is based on the free-body modal solution and the fluctuating pressure characterization. We combine the dynamics-related knowledge to carry out the theoretical derivation, and the entire calculation method has been granted a patent from the State Intellectual Property Office of China.


#### 2.3.2. Response of Elastic Displacement

For theith-order elastic modal coordinate in equation (12), we can get the corresponding motion equation by introducing the modal damping ratio as follows:(13)η¨Ei+2ωEiξEiη˙Ei+ωEi2ηEi=ϕEiTf,where ηEi is the ith-order elastic modal coordinate; ωEi and ϕEi are the ith-order elastic natural circular frequency and the elastic modal shape, respectively; and ξEi is the ith-order modal damping ratio. Substituting equation (8) into equation (13), we have(14)η¨Ei+2ωEiζEiη˙Ei+ωEi2ηEi=∑l=1rγEilpl,where γEil is the ith-order elastic modal participation factor of fluctuating pressure pl, γEil=ϕEiTfl.The modal participation factor is a crucial parameter in random vibration analysis of the free body, which transforms the fluctuating pressure information into the modal coordinate system and represents the influence degree of each fluctuating pressure on each modal coordinate. We can see that the discrete load vectorfl calculated in Section 2.2 is just devoted to the determination of modal participation factor [38] in equation (14).In the case of stationary random process, the steady-state solution of equation (14) in frequency domains is as follows:(15)ηEiω=HEiω∑l=1rγEilplω,where r is the number of fluctuating pressures and HEiω is the frequency response function (FRF) of the ith elastic mode. The expression of HEiω is as follows:(16)HEiω=1ωEi2−ω2+i2ξEiωEiω.Using the modal superposition method [39], we can get the steady-state solution of elastic displacement for the kth DOF as follows:(17)uEkω=∑i=1m2ϕEikηEiω,where m2 is the number of elastic modal shapes participated in the modal superposition and ϕEik is the kth element of ϕEi. According to the definition of auto-PSD in stationary ergodic random vibration, the auto-PSD of the kth DOF, denoted as SuEkuEkω, can be written as follows:(18)SuEkuEkω=limT⟶∞12TuEkωuEk∗ω.Substituting equation (15) into equation (17), and then substituting equation (17) into equation (18), we get(19)SuEkuEkω=∑l=1r∑m=1r∑i=1m2∑j=1m2ϕEikϕEjkγEilγEjmHEiωHEj∗ωSplpmω,where Splpmω is just the PSD of pl and pm and whose expression has been shown in equation (5). Equation (19) gives the relationship between the auto-PSD of the kth DOF and the input PSD of fluctuating pressures, which constitute the core theory formula of random vibration analysis. Both the correlated case and uncorrelated one of fluctuating pressures are taken into consideration in equation (19). For the uncorrelated case, Splpmω is equal to zero, and the corresponding computing term in equation (19) will be abandoned without participating in the superposition process in the program implementation.


#### 2.3.3. Response of Rigid Displacement

Generally, the response of rigid body motion [40] includes the nonstationary part and the stationary part, of which the latter one is focused on here. The PSD responses due to the rigid motion can be deduced according to the decoupled equation (11) and the theoretical derivation of elastic responses in Section 2.3.2. As mentioned above, the main difference between equations (11) and (12) lies in that no stiffness term appears in equation (11) and the natural circle frequencies for the rigid motion are zero. For the ith rigid body mode, the decoupled equation corresponding to equation (14) can be written as follows:(20)η¨Ri=∑l=1rγRilpl,where ηRi is the ith rigid modal coordinate and γRil is the ith rigid modal participation factor of fluctuating pressure pl, with γRil=ϕRiTfl.Being similar to equation (15), the steady-state solution of equation (20) for the stationary random process in frequency domains can be written as follows:(21)ηRiω=HRiω∑l=1rγRilplω,where HRiω is the FRF of the ith rigid mode. Due to the zero natural frequency characteristic of rigid mode, the expression of HRiω is degenerated as follows:(22)HRiω=−1ω2.It is found in equation (22) that the FRF of rigid mode is a real number with a constant variable of −1/ω2. Being similar with the theoretical derivation of elastic response, we can get the PSD expression of rigid motion by replacing the subscript E with R and replacing the summation superscript m2 with m1. The PSD of rigid response for the kth DOF is as follows:(23)SuRkuRkω=∑l=1r∑m=1r∑i=1m1∑j=1m1ϕRikϕRjkγRilγRjmHRiωHRj∗ωSplpmω.According to equations (22) and (23), we can find some interesting characteristics for random vibration of the free body induced by the rigid motion:(1)
Due to the zero natural frequency character of rigid body mode, the PSD of displacement around the zero frequency tends towards infinity, that is, the rigid PSD response presents divergent behavior at zero frequency. According to the PSD relationship between displacement and velocitySu˙u˙ω=ω2Suuω as well as that between displacement and acceleration Su˙u˙ω=ω4Suuω, we can find that the PSD of velocity on zero frequency is still divergent; however, the PSD of acceleration on zero frequency is convergent.(2)
Although the PSDs of displacement and velocity are divergent on zero frequency, they are convergent on other frequencies greater than zero. Moreover, rigid body modes have an obvious effect on the PSD response in nonzero frequency domains, which often play a leading role around the frequency of zero.


#### 2.3.4. Overall Response of Displacement

In Sections2.3.2 and 2.3.3, the theoretical calculating formulas of PSD responses for elastic deformation and rigid body motion were derived, respectively. The overall PSD response of displacement considering both the rigid body motion and elastic deformation can also be derived according to the deriving process of the elastic case. Getting rid of the subscript E in (19) and replacing the summation superscript m2 with mt (mt = m1 + m2), we can get the total PSD response of displacement for the kth DOF (denoted as Sukukω) as follows:(24)Sukukω=∑l=1r∑m=1r∑i=1mt∑j=1mtϕikϕjkγilγjmHiωHj∗ωSplpmω,where ϕik is the kth element of the ith-order modal shape vector, with following relationships with the rigid body mode and the elastic one:(25)ϕik=ϕRik,fori=1,2,…,m1,ϕi+m1k=ϕEik,fori=1,2,…,m2.Similarly, the modal participation factorγil and the FRF Hi(ω) in equation (24) are determined according to the following relationships:(26)γil=γRil,fori=1,2,…,m1,γi+m1l=γEil,fori=1,2,…,m2,Hiω=HRiω,fori=1,2,…,m1,Hi+m1ω=HEiω,fori=1,2,…,m2.It needs to be noted that the overall PSD response is not simply equal to the sum of the rigid PSD response and the elastic one; instead of that, a cross term is included in equation (24). The cross term, noted as SuERkuERkω, reflects the coupling effect between rigid body mode and elastic one. Comparing the expressions among equations (19), (23), and (24), we can get the expression of the cross term as follows:(27)SuERkuERkω=∑l=1r∑m=1r∑i=1m1∑j=m1+1mtϕikϕjkγilγjmHiωHj∗ωSplpmω+∑l=1r∑m=1r∑i=m1+1mt∑j=1m1ϕikϕjkγilγjmHiωHj∗ωSplpmω.Thus far, we have achieved the unified relational expressions between the displacement PSD responses of each DOF in three different cases and input PSDs of fluctuating pressure. The coupling relationship between the rigid body modes and the elastic ones is revealed in equation (27) for overall PSD response. As mentioned in Section 2.3.3, the PSD responses for velocity and acceleration can be obtained by multiplying a multiple of ω2 and ω4 on the PSD response of displacement.


#### 2.3.5. Root-Mean-Square (RMS) Response

The RMS response [41] is an important full-field output for random vibration, which describes the distributing information of equivalent displacement, velocity, or acceleration in the whole structure. Physically, the RMS response is obtained by making numerical integration of the PSD curve on each DOF. Taking the overall displacement as an example, the RMS response of the kth DOF, noted as σuk, can be expressed as follows:(28)σuk2=∫ωsωeSukukωdω=∫ωsωe∑l=1r∑m=1r∑i=1mt∑j=1mtϕikϕjkγilγjmHiωHj∗ωSplpmωdω,where ωs is the starting circle frequency and ωe is the ending circle frequency, respectively. In order to reduce the floating-point computations in equation (28), an optimization measure is adopted in this paper. Since the modal shapes do not relate to the integral frequency ω, they can be separated from the integral sign as follows:(29)σuk2=∑i=1mt∑j=1mtϕikϕjkQijwithQij=∫ωsωe∑l=1r∑m=1rγilγjmHiωHj∗ωSplpmωdω.Actually,Qij can be taken as an element of mt × mt covariance matrix in the modal coordinate space. The modal covariance matrix is first calculated, and the RMS response is then solved. In this way, the numbers of floating-point computations are greatly reduced at the cost of storing an mt × mt covariance matrix, whose memory space can be ignored for large-scale computations.The computational complexity between equation (29) and the optimized equation (29) can be roughly analyzed. It is assumed that there are n DOFs and p discrete frequency points in a structure model. Taking the quadruple summation term ϕikϕjkγilγjmHiωHj∗ωSplpmω in equation (28) as the basic calculation unit, it needs about n × p × r × r × mt × mt floating-point operations; instead of that, taking the double summation items ϕikϕjkQij and γilγjmHiωHj∗ωSplpmω in equation (29) as the basic calculation units, it just needs about (n × mt × mt + p × r × r) floating-point operations.For example, if the value of DOF for a structuren is 10,000, the value of frequency dispersion point p is 1,000, the value of multi-point excitation r is 10, and the value of modal participation mt is 100, the floating-point operations for equation (28) are about 1013, while those for equation (29) are only about 108. That is to say, the computational time for equation (29) is just 1/100,000 of that for equation (28).The optimization for the RMS solution can be understood as separating the integral variable with the modal spatial variable in the integral term. The integral result is obtained in the smaller modal space and superposed as a modal contribution to the physical space, which converts the original multiplicative floating-point operations into the sum of two smaller multiplication floating-point operations, thereby reducing the computational complexity.Through the above processing, we greatly reduce the amount of floating-point calculations for large-scale random vibration analysis. The following related examples in Section4 also show that the time to solve random vibration has been greatly reduced compared with modal analysis. This is important theoretical progress of our algorithms.Besides, the computing frequency point can be dispersed with a uniform discrete way or dispersed with an adaptive discretization based on the natural frequency and modal damping ratio, that is, we can increase the frequency points in the vicinity of the natural modal frequency and reduce the frequency points in the nonresonant area. This is because that the energy excited by external loads is mainly determined by the response near the natural frequency of the structure and the width of the resonance peak response curve is mainly determined by the structural damping. Generally, adaptive discretization achieves higher simulation accuracy with fewer frequency discrete points. For large-scale calculations, the amount of floating-point operations can also be reduced.


### 2.4. Algorithm Design of Random Vibration for Free Body

According to the relevant theoretical analysis in the above sections, as well as the characteristics of the JAUMIN framework and the PANDA platform, the algorithm of random vibration analysis for free bodies under the action of multi-point fluctuating pressure can be designed as follows:(1)
Finite Element Modeling and Matrix Assembly. According to the characteristics of the free-body structure and the distribution area of fluctuating pressures, the basic finite element model is constructed using the JAUMIN framework. Based on the PANDA platform, the parallel distributed mass matrix M and stiffness matrix K are constructed according to the discrete elements and physical properties. And then, the load coefficient vector fi is obtained according to the action areas of fluctuating pressure.(2)
Based on the mass matrixM and the stiffness matrix K, we construct the generalized eigenvalue problem (equation (1)) and then use the spectral transform (equation (2)) to carry out the parallel modal analysis to obtain the modal frequencies and the modal shapes, where the modal shapes are orthonormalized ones with the mass matrix M.(3)
According to the load coefficient vectorfi and the mode shape of each order, the modal participation factor of each fluctuating pressure is calculated.(4)
Taking the DOF of the nodes to be calculated as the first cycle and taking the discrete frequency of input spectrum curves as the second cycle, we calculate the PSD and RMS responses of displacement, velocity, or acceleration of each DOF and on each discrete frequency.(i)
According to the discrete frequency of the curve and the modal damping ratio corresponding to each mode, the decoupling equation in the modal coordinate system is constructed, and the real and imaginary parts of the FRF corresponding to each mode are calculated, respectively(ii)
According to the modal shape value of the respective DOF, the frequency response function of the discrete frequency point as well as the modal participation factor, the PSD responses of the elastic displacement, the rigid body displacement, and the overall displacement are calculated by combining equations (19), (23), and (24), respectively(iii)
The auto-PSDs of velocity and acceleration are calculated by the conversion relationship between displacement, velocity, and acceleration(iv)
The full-field RMS response is calculated based on equation (29)(5)
The outputs of the auto-PSD curves for the concerned nodes and the full-field RMS response are given.In the above algorithm design, we can solve the auto-PSD and RMS response of partial concerned points in structures and can also solve those of all DOFs in structures. The solution of each DOF can be achieved separately in the respective CPU processes according to the region decomposition of finite element models, that is, this solution process is a natural parallel one. It should be noted that if the auto-PSD curves of the overall DOFs of the nodes are solved, the demands of storage will be catastrophic for large-scale computation. In the actual calculation, we mainly consider calculating the RMS response of the whole field and the auto-PSD curve of a small number of points of interest.


## 2.1. Modal Analysis of Free Body

In virtue of finite element discretization and matrix assembly, the generalized eigenvalue equation corresponding to the modal analysis of free-body structures can be written as follows:(1)Kx=ω2Mx,where M and K are the mass matrix and the stiffness one, respectively; ω is the natural circular frequency; and x is the modal shape vector. Due to the totally free-flying status for vehicles in the entry atmosphere, the modal shapes of which include six rigid motion modes with ω = 0, that is, three translational ones and three rotational ones along the coordinate direction of X, Y, and Z. For free-body modal analysis, the key is to effectively solve the problem of eigenvalues caused by the singularity of the stiffness matrix and to effectively separate the rigid motion modes and elastic modes from the free-body modes.The existing mainstream algorithms for solving large-scale modal analysis are mainly subspace iterative methods [27, 29], such as the Krylov–Schur method [30] and the Jacobi–Davidson method [31]. However, these methods are mainly suitable for nonsingular matrix situations; for a free-body structure, the stiffness matrix K is a singular matrix. For this reason, the stiffness matrix needs to be processed to eliminate the singularity. The method we adopted in this paper is to introduce spectral transformation technology [32]. By introducing a positive shift value of σ into equation (1), we have(2)K+σMx=ω2+σMx.Taking (K + σM) as a matrix operator and substituting (K + σM) for K in modal analysis, the singularity of equation (1) can be eliminated, and equation (2) can be solved by subspace iterative methods. Comparing equation (1) with equation (2), it is found that whose corresponding modal shape vectors are unaltered, and the only eigenvalue difference between which is σ. Thus, the eigenpairs of equation (1) can be easily obtained by a transformation when equation (2) is computed by subspace iterative methods. In this way, the singularity problem of the free-body stiffness matrix is effectively solved.After above solution, we can getmt lowest-order eigenpairs of [(ωi, ϕi), i = 1,2, …,mt]. These eigenpairs will be used for the mode superposition of random vibration. Among these lowest eigenpairs, the first six ones are rigid body modes with ωi = 0, and the rest are elastic modes with ωi > 0. For the convenience of subsequent theoretical derivation, these eigenpairs are divided into two parts, and all output modal shape vectors are orthonormalized with mass matrix M, being satisfied with the following condition:(3)ϕiTMϕj=1i=j,0i≠j.Taking into account that the free-body structure may have a partially restricted degree of freedom, we usem1 and m2 to represent the number of rigid body modes and elastic modes, respectively, where m1 is equal to 6 for reentry vehicles, and for other noncompletely free cases, m1 is less than 6. Denoting the rigid modal shape as ϕRi and the elastic modal shape as ϕEi, we can get(4)ϕRi=ϕi,fori=1,2,…,m1,ϕEi=ϕi+m1,fori=1,2,…,m2.


## 2.2. Description of Fluctuating Pressures

Generally, the reentry process of a vehicle is treated as a stationary random process in the finite element analysis of random vibration [33]. The basic inputs for random vibration are the auto-PSDs and cross-PSDs of fluctuating pressures. According to these inputs, the auto-PSDs of displacement, the velocity, and acceleration for concerned points inside the structures and the equivalent root-mean-square (RMS) responses of full field are calculated via random vibration analysis.It is obvious that the fluctuating pressure distributions on the external surface of reentry vehicles are extremely complex. Consequently, it is impossible to describe the fluctuating pressure at each discrete finite element node of the outside surface, especially for the case of large-scale computations.In this paper, a treatment of geometry area decomposition on the outside surface is adopted to approximate the pressure distributions. The external surface of the vehicle is divided into numerous areas, and a uniform pressure value is specified in each area. In order to represent the real fluctuating pressure distributions more exactly, the geometry areas can be further subdivided to approximate the real loads.Without loss of generality, we assume that there arer different fluctuating pressures acting on the external surface, being denoted as p1,p2,…,pr. The acting areas of these pressures constitute the whole load region of the reentry vehicle. It should be noted that these fluctuating pressures cannot be used directly in the implementation of random vibration. Instead of that, the auto-PSDs Spipi (i = 1,2, …, r) and the cross-PSDs Spipj (i = 1,2, …, r; j = 1,2, …, r, i ≠ j) of these fluctuating pressures are the direct input data for random vibration of free body. For the stationary random process, the relationship between fluctuating pressures and their PSDs can be shown as follows [34]:(5)Splpmω=limT⟶∞12Tplωpm∗ω.Theoretically, the fluctuating pressures on two different geometry areas might be totally correlated, partly correlated, or uncorrelated. Taking each case into consideration, the input PSDs are designed as anr × r matrix as follows:(6)Sp1p1Sp1p2⋯Sp1prSp2p1Sp2p2⋯Sp2pr⋮⋮⋱⋮Sprp1Sprp2⋯Sprpr,where the diagonal elements represent the auto-PSDs and the nondiagonal elements represent the cross-PSDs. The matrix elements in equation (6) are nonzero in the condition that the fluctuating pressures between these two areas are correlated, otherwise zero.In this way, the fluctuating pressure loads on the outer surface of the reentry vehicle are treated as anr-dimensional spectral curve matrix; each element in the matrix corresponds to a spectral curve or 0. The advantage of this processing is that the matrix can describe various situations such as completely related, partially related, and completely unrelated.Although the fluctuating pressures are not the direct input for the random vibration, their distributing information are still necessary for the assembly of dynamic equations and the computation of modal participation factors. The discrete finite-element dynamic equation for the reentry vehicle with rigid motion can be described as follows:(7)Mu¨+Cu˙+Ku=f,where C is the damping matrix, u is the displacement vector for all degrees of freedom, and f is the fluctuating pressure vector. In the construction of equation (7), f is assembled via the integration of each pressure on its corresponding element surface and then combining them together. f can be described as follows:(8)f=∑l=1rflpl,where fl is the coefficient vector of discrete load with unit fluctuating pressure on the acting area of pl.In order to determinefl, a unit pressure load in this paper is supposed to apply to the acting area of pl, and finite element integration and discretization are carried out to form fl. It can be seen in the following theoretical derivations and implementations that neither the fluctuating pressure pl nor the total discretized load f needs to be explicitly given or calculated. The actual input for random vibration of the free body is the load factor vector fl and the PSDs of fluctuating pressures.


## 2.3. Theories and Algorithm for Random Vibration of the Free Body

### 2.3.1. Modal Decoupling and Separation of Rigid Displacements and Elastic Ones

The displacement vectoru in equation (7) can be divided into two parts: the rigid displacement vector uR and the elastic displacement vector uE [35]. When the modal superposition method is adopted, the two parts of displacement can be described by rigid modal shapes and elastic ones:(9)u=uR+uE=ΦRηR+ΦEηE,where ΦR and ΦE are the rigid modal matrix and the elastic one, respectively, and ηR and ηE are the rigid modal coordinate vector and the elastic one, respectively.It should be noted that equation (9) is approximately established due to the modal truncation in the modal superposition method [36]. Usually, the more the modal shapes we used, the better the precision we achieve. Since K is singular (KuR = 0), we achieve the following equation by substituting equation (9) into equation (7):(10)Mu¨E+Cu˙E+KuE=f−Mu¨R−Cu˙R.According to the orthogonality (equation (3)) and the small damping assumption, we can get two sets of decoupled equations by substituting uE=ΦEηE and uR=ΦRηR into equation (10) and premultiplying both sides of equation (10) by ΦΕΤ and ΦRΤ, as follows:(11)MRη¨R+CRη˙R=ΦRΤf,(12)MEη¨E+CEη˙E+KEηE=ΦΕΤf,where MR=ΦRTMΦR=IR, ME=ΦΕΤMΦE=IE, CR=ΦRTCΦR, CE=ΦETCΦE, and ΚE=ΦETΚΦE.Equations (11) and (12) give the rigid modal coordinate description [37] and the elastic modal coordinate description of equation (7), respectively. In other words, the dynamic equation of the free body is separated effectively with the rigid part and the elastic one. Comparing equation (11) with equation (12), we can find that both equations are nearly the same except for a missing singular stiffness term in equation (11). In the following, the relatively complex case of the elastic response of the structure is deduced first, and then, the rigid response and the rigid-elastic coupling response are deduced.It should be noted that the following derivation process of elastic displacement response, rigid body displacement response, and the total response of the coupling of the two is based on the free-body modal solution and the fluctuating pressure characterization. We combine the dynamics-related knowledge to carry out the theoretical derivation, and the entire calculation method has been granted a patent from the State Intellectual Property Office of China.


### 2.3.2. Response of Elastic Displacement

For theith-order elastic modal coordinate in equation (12), we can get the corresponding motion equation by introducing the modal damping ratio as follows:(13)η¨Ei+2ωEiξEiη˙Ei+ωEi2ηEi=ϕEiTf,where ηEi is the ith-order elastic modal coordinate; ωEi and ϕEi are the ith-order elastic natural circular frequency and the elastic modal shape, respectively; and ξEi is the ith-order modal damping ratio. Substituting equation (8) into equation (13), we have(14)η¨Ei+2ωEiζEiη˙Ei+ωEi2ηEi=∑l=1rγEilpl,where γEil is the ith-order elastic modal participation factor of fluctuating pressure pl, γEil=ϕEiTfl.The modal participation factor is a crucial parameter in random vibration analysis of the free body, which transforms the fluctuating pressure information into the modal coordinate system and represents the influence degree of each fluctuating pressure on each modal coordinate. We can see that the discrete load vectorfl calculated in Section 2.2 is just devoted to the determination of modal participation factor [38] in equation (14).In the case of stationary random process, the steady-state solution of equation (14) in frequency domains is as follows:(15)ηEiω=HEiω∑l=1rγEilplω,where r is the number of fluctuating pressures and HEiω is the frequency response function (FRF) of the ith elastic mode. The expression of HEiω is as follows:(16)HEiω=1ωEi2−ω2+i2ξEiωEiω.Using the modal superposition method [39], we can get the steady-state solution of elastic displacement for the kth DOF as follows:(17)uEkω=∑i=1m2ϕEikηEiω,where m2 is the number of elastic modal shapes participated in the modal superposition and ϕEik is the kth element of ϕEi. According to the definition of auto-PSD in stationary ergodic random vibration, the auto-PSD of the kth DOF, denoted as SuEkuEkω, can be written as follows:(18)SuEkuEkω=limT⟶∞12TuEkωuEk∗ω.Substituting equation (15) into equation (17), and then substituting equation (17) into equation (18), we get(19)SuEkuEkω=∑l=1r∑m=1r∑i=1m2∑j=1m2ϕEikϕEjkγEilγEjmHEiωHEj∗ωSplpmω,where Splpmω is just the PSD of pl and pm and whose expression has been shown in equation (5). Equation (19) gives the relationship between the auto-PSD of the kth DOF and the input PSD of fluctuating pressures, which constitute the core theory formula of random vibration analysis. Both the correlated case and uncorrelated one of fluctuating pressures are taken into consideration in equation (19). For the uncorrelated case, Splpmω is equal to zero, and the corresponding computing term in equation (19) will be abandoned without participating in the superposition process in the program implementation.


### 2.3.3. Response of Rigid Displacement

Generally, the response of rigid body motion [40] includes the nonstationary part and the stationary part, of which the latter one is focused on here. The PSD responses due to the rigid motion can be deduced according to the decoupled equation (11) and the theoretical derivation of elastic responses in Section 2.3.2. As mentioned above, the main difference between equations (11) and (12) lies in that no stiffness term appears in equation (11) and the natural circle frequencies for the rigid motion are zero. For the ith rigid body mode, the decoupled equation corresponding to equation (14) can be written as follows:(20)η¨Ri=∑l=1rγRilpl,where ηRi is the ith rigid modal coordinate and γRil is the ith rigid modal participation factor of fluctuating pressure pl, with γRil=ϕRiTfl.Being similar to equation (15), the steady-state solution of equation (20) for the stationary random process in frequency domains can be written as follows:(21)ηRiω=HRiω∑l=1rγRilplω,where HRiω is the FRF of the ith rigid mode. Due to the zero natural frequency characteristic of rigid mode, the expression of HRiω is degenerated as follows:(22)HRiω=−1ω2.It is found in equation (22) that the FRF of rigid mode is a real number with a constant variable of −1/ω2. Being similar with the theoretical derivation of elastic response, we can get the PSD expression of rigid motion by replacing the subscript E with R and replacing the summation superscript m2 with m1. The PSD of rigid response for the kth DOF is as follows:(23)SuRkuRkω=∑l=1r∑m=1r∑i=1m1∑j=1m1ϕRikϕRjkγRilγRjmHRiωHRj∗ωSplpmω.According to equations (22) and (23), we can find some interesting characteristics for random vibration of the free body induced by the rigid motion:(1)
Due to the zero natural frequency character of rigid body mode, the PSD of displacement around the zero frequency tends towards infinity, that is, the rigid PSD response presents divergent behavior at zero frequency. According to the PSD relationship between displacement and velocitySu˙u˙ω=ω2Suuω as well as that between displacement and acceleration Su˙u˙ω=ω4Suuω, we can find that the PSD of velocity on zero frequency is still divergent; however, the PSD of acceleration on zero frequency is convergent.(2)
Although the PSDs of displacement and velocity are divergent on zero frequency, they are convergent on other frequencies greater than zero. Moreover, rigid body modes have an obvious effect on the PSD response in nonzero frequency domains, which often play a leading role around the frequency of zero.


### 2.3.4. Overall Response of Displacement

In Sections2.3.2 and 2.3.3, the theoretical calculating formulas of PSD responses for elastic deformation and rigid body motion were derived, respectively. The overall PSD response of displacement considering both the rigid body motion and elastic deformation can also be derived according to the deriving process of the elastic case. Getting rid of the subscript E in (19) and replacing the summation superscript m2 with mt (mt = m1 + m2), we can get the total PSD response of displacement for the kth DOF (denoted as Sukukω) as follows:(24)Sukukω=∑l=1r∑m=1r∑i=1mt∑j=1mtϕikϕjkγilγjmHiωHj∗ωSplpmω,where ϕik is the kth element of the ith-order modal shape vector, with following relationships with the rigid body mode and the elastic one:(25)ϕik=ϕRik,fori=1,2,…,m1,ϕi+m1k=ϕEik,fori=1,2,…,m2.Similarly, the modal participation factorγil and the FRF Hi(ω) in equation (24) are determined according to the following relationships:(26)γil=γRil,fori=1,2,…,m1,γi+m1l=γEil,fori=1,2,…,m2,Hiω=HRiω,fori=1,2,…,m1,Hi+m1ω=HEiω,fori=1,2,…,m2.It needs to be noted that the overall PSD response is not simply equal to the sum of the rigid PSD response and the elastic one; instead of that, a cross term is included in equation (24). The cross term, noted as SuERkuERkω, reflects the coupling effect between rigid body mode and elastic one. Comparing the expressions among equations (19), (23), and (24), we can get the expression of the cross term as follows:(27)SuERkuERkω=∑l=1r∑m=1r∑i=1m1∑j=m1+1mtϕikϕjkγilγjmHiωHj∗ωSplpmω+∑l=1r∑m=1r∑i=m1+1mt∑j=1m1ϕikϕjkγilγjmHiωHj∗ωSplpmω.Thus far, we have achieved the unified relational expressions between the displacement PSD responses of each DOF in three different cases and input PSDs of fluctuating pressure. The coupling relationship between the rigid body modes and the elastic ones is revealed in equation (27) for overall PSD response. As mentioned in Section 2.3.3, the PSD responses for velocity and acceleration can be obtained by multiplying a multiple of ω2 and ω4 on the PSD response of displacement.


### 2.3.5. Root-Mean-Square (RMS) Response

The RMS response [41] is an important full-field output for random vibration, which describes the distributing information of equivalent displacement, velocity, or acceleration in the whole structure. Physically, the RMS response is obtained by making numerical integration of the PSD curve on each DOF. Taking the overall displacement as an example, the RMS response of the kth DOF, noted as σuk, can be expressed as follows:(28)σuk2=∫ωsωeSukukωdω=∫ωsωe∑l=1r∑m=1r∑i=1mt∑j=1mtϕikϕjkγilγjmHiωHj∗ωSplpmωdω,where ωs is the starting circle frequency and ωe is the ending circle frequency, respectively. In order to reduce the floating-point computations in equation (28), an optimization measure is adopted in this paper. Since the modal shapes do not relate to the integral frequency ω, they can be separated from the integral sign as follows:(29)σuk2=∑i=1mt∑j=1mtϕikϕjkQijwithQij=∫ωsωe∑l=1r∑m=1rγilγjmHiωHj∗ωSplpmωdω.Actually,Qij can be taken as an element of mt × mt covariance matrix in the modal coordinate space. The modal covariance matrix is first calculated, and the RMS response is then solved. In this way, the numbers of floating-point computations are greatly reduced at the cost of storing an mt × mt covariance matrix, whose memory space can be ignored for large-scale computations.The computational complexity between equation (29) and the optimized equation (29) can be roughly analyzed. It is assumed that there are n DOFs and p discrete frequency points in a structure model. Taking the quadruple summation term ϕikϕjkγilγjmHiωHj∗ωSplpmω in equation (28) as the basic calculation unit, it needs about n × p × r × r × mt × mt floating-point operations; instead of that, taking the double summation items ϕikϕjkQij and γilγjmHiωHj∗ωSplpmω in equation (29) as the basic calculation units, it just needs about (n × mt × mt + p × r × r) floating-point operations.For example, if the value of DOF for a structuren is 10,000, the value of frequency dispersion point p is 1,000, the value of multi-point excitation r is 10, and the value of modal participation mt is 100, the floating-point operations for equation (28) are about 1013, while those for equation (29) are only about 108. That is to say, the computational time for equation (29) is just 1/100,000 of that for equation (28).The optimization for the RMS solution can be understood as separating the integral variable with the modal spatial variable in the integral term. The integral result is obtained in the smaller modal space and superposed as a modal contribution to the physical space, which converts the original multiplicative floating-point operations into the sum of two smaller multiplication floating-point operations, thereby reducing the computational complexity.Through the above processing, we greatly reduce the amount of floating-point calculations for large-scale random vibration analysis. The following related examples in Section4 also show that the time to solve random vibration has been greatly reduced compared with modal analysis. This is important theoretical progress of our algorithms.Besides, the computing frequency point can be dispersed with a uniform discrete way or dispersed with an adaptive discretization based on the natural frequency and modal damping ratio, that is, we can increase the frequency points in the vicinity of the natural modal frequency and reduce the frequency points in the nonresonant area. This is because that the energy excited by external loads is mainly determined by the response near the natural frequency of the structure and the width of the resonance peak response curve is mainly determined by the structural damping. Generally, adaptive discretization achieves higher simulation accuracy with fewer frequency discrete points. For large-scale calculations, the amount of floating-point operations can also be reduced.


## 2.3.1. Modal Decoupling and Separation of Rigid Displacements and Elastic Ones

The displacement vectoru in equation (7) can be divided into two parts: the rigid displacement vector uR and the elastic displacement vector uE [35]. When the modal superposition method is adopted, the two parts of displacement can be described by rigid modal shapes and elastic ones:(9)u=uR+uE=ΦRηR+ΦEηE,where ΦR and ΦE are the rigid modal matrix and the elastic one, respectively, and ηR and ηE are the rigid modal coordinate vector and the elastic one, respectively.It should be noted that equation (9) is approximately established due to the modal truncation in the modal superposition method [36]. Usually, the more the modal shapes we used, the better the precision we achieve. Since K is singular (KuR = 0), we achieve the following equation by substituting equation (9) into equation (7):(10)Mu¨E+Cu˙E+KuE=f−Mu¨R−Cu˙R.According to the orthogonality (equation (3)) and the small damping assumption, we can get two sets of decoupled equations by substituting uE=ΦEηE and uR=ΦRηR into equation (10) and premultiplying both sides of equation (10) by ΦΕΤ and ΦRΤ, as follows:(11)MRη¨R+CRη˙R=ΦRΤf,(12)MEη¨E+CEη˙E+KEηE=ΦΕΤf,where MR=ΦRTMΦR=IR, ME=ΦΕΤMΦE=IE, CR=ΦRTCΦR, CE=ΦETCΦE, and ΚE=ΦETΚΦE.Equations (11) and (12) give the rigid modal coordinate description [37] and the elastic modal coordinate description of equation (7), respectively. In other words, the dynamic equation of the free body is separated effectively with the rigid part and the elastic one. Comparing equation (11) with equation (12), we can find that both equations are nearly the same except for a missing singular stiffness term in equation (11). In the following, the relatively complex case of the elastic response of the structure is deduced first, and then, the rigid response and the rigid-elastic coupling response are deduced.It should be noted that the following derivation process of elastic displacement response, rigid body displacement response, and the total response of the coupling of the two is based on the free-body modal solution and the fluctuating pressure characterization. We combine the dynamics-related knowledge to carry out the theoretical derivation, and the entire calculation method has been granted a patent from the State Intellectual Property Office of China.


## 2.3.2. Response of Elastic Displacement

For theith-order elastic modal coordinate in equation (12), we can get the corresponding motion equation by introducing the modal damping ratio as follows:(13)η¨Ei+2ωEiξEiη˙Ei+ωEi2ηEi=ϕEiTf,where ηEi is the ith-order elastic modal coordinate; ωEi and ϕEi are the ith-order elastic natural circular frequency and the elastic modal shape, respectively; and ξEi is the ith-order modal damping ratio. Substituting equation (8) into equation (13), we have(14)η¨Ei+2ωEiζEiη˙Ei+ωEi2ηEi=∑l=1rγEilpl,where γEil is the ith-order elastic modal participation factor of fluctuating pressure pl, γEil=ϕEiTfl.The modal participation factor is a crucial parameter in random vibration analysis of the free body, which transforms the fluctuating pressure information into the modal coordinate system and represents the influence degree of each fluctuating pressure on each modal coordinate. We can see that the discrete load vectorfl calculated in Section 2.2 is just devoted to the determination of modal participation factor [38] in equation (14).In the case of stationary random process, the steady-state solution of equation (14) in frequency domains is as follows:(15)ηEiω=HEiω∑l=1rγEilplω,where r is the number of fluctuating pressures and HEiω is the frequency response function (FRF) of the ith elastic mode. The expression of HEiω is as follows:(16)HEiω=1ωEi2−ω2+i2ξEiωEiω.Using the modal superposition method [39], we can get the steady-state solution of elastic displacement for the kth DOF as follows:(17)uEkω=∑i=1m2ϕEikηEiω,where m2 is the number of elastic modal shapes participated in the modal superposition and ϕEik is the kth element of ϕEi. According to the definition of auto-PSD in stationary ergodic random vibration, the auto-PSD of the kth DOF, denoted as SuEkuEkω, can be written as follows:(18)SuEkuEkω=limT⟶∞12TuEkωuEk∗ω.Substituting equation (15) into equation (17), and then substituting equation (17) into equation (18), we get(19)SuEkuEkω=∑l=1r∑m=1r∑i=1m2∑j=1m2ϕEikϕEjkγEilγEjmHEiωHEj∗ωSplpmω,where Splpmω is just the PSD of pl and pm and whose expression has been shown in equation (5). Equation (19) gives the relationship between the auto-PSD of the kth DOF and the input PSD of fluctuating pressures, which constitute the core theory formula of random vibration analysis. Both the correlated case and uncorrelated one of fluctuating pressures are taken into consideration in equation (19). For the uncorrelated case, Splpmω is equal to zero, and the corresponding computing term in equation (19) will be abandoned without participating in the superposition process in the program implementation.


## 2.3.3. Response of Rigid Displacement

Generally, the response of rigid body motion [40] includes the nonstationary part and the stationary part, of which the latter one is focused on here. The PSD responses due to the rigid motion can be deduced according to the decoupled equation (11) and the theoretical derivation of elastic responses in Section 2.3.2. As mentioned above, the main difference between equations (11) and (12) lies in that no stiffness term appears in equation (11) and the natural circle frequencies for the rigid motion are zero. For the ith rigid body mode, the decoupled equation corresponding to equation (14) can be written as follows:(20)η¨Ri=∑l=1rγRilpl,where ηRi is the ith rigid modal coordinate and γRil is the ith rigid modal participation factor of fluctuating pressure pl, with γRil=ϕRiTfl.Being similar to equation (15), the steady-state solution of equation (20) for the stationary random process in frequency domains can be written as follows:(21)ηRiω=HRiω∑l=1rγRilplω,where HRiω is the FRF of the ith rigid mode. Due to the zero natural frequency characteristic of rigid mode, the expression of HRiω is degenerated as follows:(22)HRiω=−1ω2.It is found in equation (22) that the FRF of rigid mode is a real number with a constant variable of −1/ω2. Being similar with the theoretical derivation of elastic response, we can get the PSD expression of rigid motion by replacing the subscript E with R and replacing the summation superscript m2 with m1. The PSD of rigid response for the kth DOF is as follows:(23)SuRkuRkω=∑l=1r∑m=1r∑i=1m1∑j=1m1ϕRikϕRjkγRilγRjmHRiωHRj∗ωSplpmω.According to equations (22) and (23), we can find some interesting characteristics for random vibration of the free body induced by the rigid motion:(1)
Due to the zero natural frequency character of rigid body mode, the PSD of displacement around the zero frequency tends towards infinity, that is, the rigid PSD response presents divergent behavior at zero frequency. According to the PSD relationship between displacement and velocitySu˙u˙ω=ω2Suuω as well as that between displacement and acceleration Su˙u˙ω=ω4Suuω, we can find that the PSD of velocity on zero frequency is still divergent; however, the PSD of acceleration on zero frequency is convergent.(2)
Although the PSDs of displacement and velocity are divergent on zero frequency, they are convergent on other frequencies greater than zero. Moreover, rigid body modes have an obvious effect on the PSD response in nonzero frequency domains, which often play a leading role around the frequency of zero.


## 2.3.4. Overall Response of Displacement

In Sections2.3.2 and 2.3.3, the theoretical calculating formulas of PSD responses for elastic deformation and rigid body motion were derived, respectively. The overall PSD response of displacement considering both the rigid body motion and elastic deformation can also be derived according to the deriving process of the elastic case. Getting rid of the subscript E in (19) and replacing the summation superscript m2 with mt (mt = m1 + m2), we can get the total PSD response of displacement for the kth DOF (denoted as Sukukω) as follows:(24)Sukukω=∑l=1r∑m=1r∑i=1mt∑j=1mtϕikϕjkγilγjmHiωHj∗ωSplpmω,where ϕik is the kth element of the ith-order modal shape vector, with following relationships with the rigid body mode and the elastic one:(25)ϕik=ϕRik,fori=1,2,…,m1,ϕi+m1k=ϕEik,fori=1,2,…,m2.Similarly, the modal participation factorγil and the FRF Hi(ω) in equation (24) are determined according to the following relationships:(26)γil=γRil,fori=1,2,…,m1,γi+m1l=γEil,fori=1,2,…,m2,Hiω=HRiω,fori=1,2,…,m1,Hi+m1ω=HEiω,fori=1,2,…,m2.It needs to be noted that the overall PSD response is not simply equal to the sum of the rigid PSD response and the elastic one; instead of that, a cross term is included in equation (24). The cross term, noted as SuERkuERkω, reflects the coupling effect between rigid body mode and elastic one. Comparing the expressions among equations (19), (23), and (24), we can get the expression of the cross term as follows:(27)SuERkuERkω=∑l=1r∑m=1r∑i=1m1∑j=m1+1mtϕikϕjkγilγjmHiωHj∗ωSplpmω+∑l=1r∑m=1r∑i=m1+1mt∑j=1m1ϕikϕjkγilγjmHiωHj∗ωSplpmω.Thus far, we have achieved the unified relational expressions between the displacement PSD responses of each DOF in three different cases and input PSDs of fluctuating pressure. The coupling relationship between the rigid body modes and the elastic ones is revealed in equation (27) for overall PSD response. As mentioned in Section 2.3.3, the PSD responses for velocity and acceleration can be obtained by multiplying a multiple of ω2 and ω4 on the PSD response of displacement.


## 2.3.5. Root-Mean-Square (RMS) Response

The RMS response [41] is an important full-field output for random vibration, which describes the distributing information of equivalent displacement, velocity, or acceleration in the whole structure. Physically, the RMS response is obtained by making numerical integration of the PSD curve on each DOF. Taking the overall displacement as an example, the RMS response of the kth DOF, noted as σuk, can be expressed as follows:(28)σuk2=∫ωsωeSukukωdω=∫ωsωe∑l=1r∑m=1r∑i=1mt∑j=1mtϕikϕjkγilγjmHiωHj∗ωSplpmωdω,where ωs is the starting circle frequency and ωe is the ending circle frequency, respectively. In order to reduce the floating-point computations in equation (28), an optimization measure is adopted in this paper. Since the modal shapes do not relate to the integral frequency ω, they can be separated from the integral sign as follows:(29)σuk2=∑i=1mt∑j=1mtϕikϕjkQijwithQij=∫ωsωe∑l=1r∑m=1rγilγjmHiωHj∗ωSplpmωdω.Actually,Qij can be taken as an element of mt × mt covariance matrix in the modal coordinate space. The modal covariance matrix is first calculated, and the RMS response is then solved. In this way, the numbers of floating-point computations are greatly reduced at the cost of storing an mt × mt covariance matrix, whose memory space can be ignored for large-scale computations.The computational complexity between equation (29) and the optimized equation (29) can be roughly analyzed. It is assumed that there are n DOFs and p discrete frequency points in a structure model. Taking the quadruple summation term ϕikϕjkγilγjmHiωHj∗ωSplpmω in equation (28) as the basic calculation unit, it needs about n × p × r × r × mt × mt floating-point operations; instead of that, taking the double summation items ϕikϕjkQij and γilγjmHiωHj∗ωSplpmω in equation (29) as the basic calculation units, it just needs about (n × mt × mt + p × r × r) floating-point operations.For example, if the value of DOF for a structuren is 10,000, the value of frequency dispersion point p is 1,000, the value of multi-point excitation r is 10, and the value of modal participation mt is 100, the floating-point operations for equation (28) are about 1013, while those for equation (29) are only about 108. That is to say, the computational time for equation (29) is just 1/100,000 of that for equation (28).The optimization for the RMS solution can be understood as separating the integral variable with the modal spatial variable in the integral term. The integral result is obtained in the smaller modal space and superposed as a modal contribution to the physical space, which converts the original multiplicative floating-point operations into the sum of two smaller multiplication floating-point operations, thereby reducing the computational complexity.Through the above processing, we greatly reduce the amount of floating-point calculations for large-scale random vibration analysis. The following related examples in Section4 also show that the time to solve random vibration has been greatly reduced compared with modal analysis. This is important theoretical progress of our algorithms.Besides, the computing frequency point can be dispersed with a uniform discrete way or dispersed with an adaptive discretization based on the natural frequency and modal damping ratio, that is, we can increase the frequency points in the vicinity of the natural modal frequency and reduce the frequency points in the nonresonant area. This is because that the energy excited by external loads is mainly determined by the response near the natural frequency of the structure and the width of the resonance peak response curve is mainly determined by the structural damping. Generally, adaptive discretization achieves higher simulation accuracy with fewer frequency discrete points. For large-scale calculations, the amount of floating-point operations can also be reduced.


## 2.4. Algorithm Design of Random Vibration for Free Body

According to the relevant theoretical analysis in the above sections, as well as the characteristics of the JAUMIN framework and the PANDA platform, the algorithm of random vibration analysis for free bodies under the action of multi-point fluctuating pressure can be designed as follows:(1)
Finite Element Modeling and Matrix Assembly. According to the characteristics of the free-body structure and the distribution area of fluctuating pressures, the basic finite element model is constructed using the JAUMIN framework. Based on the PANDA platform, the parallel distributed mass matrix M and stiffness matrix K are constructed according to the discrete elements and physical properties. And then, the load coefficient vector fi is obtained according to the action areas of fluctuating pressure.(2)
Based on the mass matrixM and the stiffness matrix K, we construct the generalized eigenvalue problem (equation (1)) and then use the spectral transform (equation (2)) to carry out the parallel modal analysis to obtain the modal frequencies and the modal shapes, where the modal shapes are orthonormalized ones with the mass matrix M.(3)
According to the load coefficient vectorfi and the mode shape of each order, the modal participation factor of each fluctuating pressure is calculated.(4)
Taking the DOF of the nodes to be calculated as the first cycle and taking the discrete frequency of input spectrum curves as the second cycle, we calculate the PSD and RMS responses of displacement, velocity, or acceleration of each DOF and on each discrete frequency.(i)
According to the discrete frequency of the curve and the modal damping ratio corresponding to each mode, the decoupling equation in the modal coordinate system is constructed, and the real and imaginary parts of the FRF corresponding to each mode are calculated, respectively(ii)
According to the modal shape value of the respective DOF, the frequency response function of the discrete frequency point as well as the modal participation factor, the PSD responses of the elastic displacement, the rigid body displacement, and the overall displacement are calculated by combining equations (19), (23), and (24), respectively(iii)
The auto-PSDs of velocity and acceleration are calculated by the conversion relationship between displacement, velocity, and acceleration(iv)
The full-field RMS response is calculated based on equation (29)(5)
The outputs of the auto-PSD curves for the concerned nodes and the full-field RMS response are given.In the above algorithm design, we can solve the auto-PSD and RMS response of partial concerned points in structures and can also solve those of all DOFs in structures. The solution of each DOF can be achieved separately in the respective CPU processes according to the region decomposition of finite element models, that is, this solution process is a natural parallel one. It should be noted that if the auto-PSD curves of the overall DOFs of the nodes are solved, the demands of storage will be catastrophic for large-scale computation. In the actual calculation, we mainly consider calculating the RMS response of the whole field and the auto-PSD curve of a small number of points of interest.


## 3. Parallel Implementation

The parallel implementation for random vibration analysis of free body is mainly based on our independently developed PANDA platform, which is a parallel adaptive finite element computing service platform based on the data structure of the JAUMIN framework. The JAUMIN framework has been developed by the China Academy of Engineering Physics since 2011, which is mainly aimed at the large-scale numerical calculation of unstructured grid. Following, a brief introduction about the JAUMIN framework and the PANDA platform is firstly given, and then the parallel implementations of random vibration analysis for free bodies are described.


### 3.1. The JAUMIN Framework

The JAUMIN framework is mainly based on the hierarchical and modular design idea, using C++, MPI [42], and OpenMP mixed programming, which contains about 120,000 lines of code. The basic architecture of the JAUMIN framework is shown in Figure 1, which is divided into three levels.Figure 1
Architecture of the JAUMIN framework.The bottom layer is the high-performance computing support layer, composed of multiple modules for high-performance computing of an unstructured grid. In this layer, the relevant modules are divided into three sublayers. The first sublayer is a toolbox module, mainly to provide some basic tools, such as operating environment management, input database, memory management, and so on. The second layer includes three modules, in which the patch hierarchy structure and patch data are mainly used for the storage of finite element grid and field data and the data communication module is mainly used for data communication between different grid regions. The third sublayer includes grid adaptive module, load balancing module, and math operation module, which are used for adaptive refinement of finite element mesh, data balance processing of each process and related mathematical operation, and so on.The middle layer is the numerical general layer, including some modules for general numerical algorithms and application toolboxes, mathematical solvers for matrix-vector operations, time integrators, mesh generator, multi-physics coupling, and so on.The top layer is the applied interface layer, mainly providing the basic program interface for different applications. JAUMIN provides C++ interfaces for parallel programming of various application software. Based on this layer, users can write their own application subroutines based on the parallel data structure of the JAUMIN framework. These subroutines together constitute the application program of JAUMIN.In order to adapt to the large-scale parallel computing environment, a two-level partitioning strategy is adopted in the JAUMIN framework, shown as in Figure2. First, the finite element model is partitioned into submeshes (such as P0,P1 in Figure 2) considering load balance; each CPU processor is assigned with a submesh. Second, the submesh is again partitioned into a group of smaller submeshes, named “patch.” The patch is composed of some mesh elements and is the basic object for memory management and numerical computation in JAUMIN. Each patch includes some basic information such as patch geometry, patch topology, and patch data, which constitutes the core data structures of the framework.Figure 2
A two-level partitioning strategy in the JAUMIN framework.JAUMIN integrates a universal unstructured-mesh data model, a pattern-based parallel computing model, and a component-based programming model. These features enable JAUMIN to effectively shield parallel implementation details and significantly reduce programming complexity, thus serving as a bridge between supercomputers and applications. JAUMIN facilitates the development of large-scale parallel unstructured mesh applications and supports applications to scale to tens of thousands of CPU cores. Recently, JAUMIN has been widely used in a series of applications such as large deformation elastoplastic hydrodynamics, material damage and fracture, large complex structural impact dynamics, structural static and modal vibration analysis, and so on. For a more detailed introduction to the JAUMIN framework, we recommend literature [26].


### 3.2. The PANDA Platform

The PANDA platform is a large-scale finite element parallel computing platform based on the data structure of the JUMIN framework, which is mainly for large-scale finite element calculation of complex engineering equipment. At present, the PANDA platform includes static, modal and vibration, impact dynamics, and other types of finite element analysis, mainly coded with C++, including codes more than 100,000 lines. For the complex engineering structure, the pre- and postprocessing software are combined with the JAUMIN framework to complete the finite element modeling and grid domain decomposition, and the PANDA platform is responsible for the parallel matrix assembly and solution according to the structural material properties and external load conditions. The basic structure of the PANDA platform and its relationship with the JAUMIN framework are shown in Figure3.Figure 3
The PANDA platform and its relationship with the JAUMIN framework.The free-body random vibration analysis mentioned in this paper is taken as an analysis module in the PANDA-StaVib software, as shown in Figure4. In addition to the static analysis [43], the analysis software PANDA-StaVib includes some structural dynamic modules such as modal analysis, single- or multi-point harmonic response analysis [44], seismic response spectrum analysis, single- or multi-point basic excitation random vibration analysis [45], free-body random vibration analysis with multi-point fluctuating pressures, and so on, covering most of dynamic frequency domain analysis types. Although the above dynamic analysis modules are relatively independent, some common features such as the calculation of modal participation factor, description of input spectral curves, and modal analysis interfaces for modal superposition can be shared with each other, which is also compatible with the large-scale program design principle of “extracting common features, stripping personality,” reducing unnecessary code duplications.Figure 4
Analysis modules in PANDA-StaVib software.The PANDA platform mainly uses the finite method for numerical calculation and also includes the main hierarchical structure such as the finite element support layer, the finite element common layer, and the application personality layer. Among them, the finite element support layer provides data support services for the software, including the registration of JAUMIN variables, model data management, degree of freedom mapping, point data output, and toolboxes commonly used in matrix-vector operations; the finite element common layer is useful for structural mechanics finite element calculations where the common part is extracted, mainly including elements, materials, field data, constraints, loads and the corresponding solution process of each analysis type; and the finite element personality layer establishes a personality analysis module for different structural analysis types, as mentioned in this paper, random vibration analysis is one of the modules.


### 3.3. Parallel Implementation of Free-Body Random Vibration Analysis

The basic flow of the free-body random vibration analysis can be described in Figure5. In virtue of domain decomposition and load balancing management of the JAUMIN framework, the data information for each domain of finite element models is evenly distributed to each processor. In each processor, the PANDA platform combines the model parameter information to generate the mass matrix M, the stiffness matrix K, and the load coefficient vector fl. The information of matrices and vectors stored in each processor is parallel distributed, which provides the basis of data parallelism for the subsequent random vibration analysis of the free body.Figure 5
The basic calculation process of free-body random vibration analysis.The natural frequency and modal shapes, including all rigid body modes and partial low order elastic modes, are obtained with the free-body modal analysis in the PANDA platform. Using these data, the frequency response function for different discrete frequencies and the modal participation factors for different modes can be calculated and provided for the random vibration response. The spectral curve manager is used for the auto-PSD and cross-PSD curve management among different fluctuating pressures to provide the load input parameters for the calculation process. Based on these inputs, the random vibration response analysis is performed on each degree of freedom using the algorithm designed in Section2.4. Finally, the results are directed to the output stream, and the equivalent RMS response for all fields as well as the auto-PSD response curves for the concerned points are outputted.In the above solution process, the main parallel stages include the model domain decomposition and the matrix assembly, modal analysis, modal superposition, vibration response calculation phase, and so on. The modal analysis is mainly carried out by the eigenvalue solution library of the JAUMIN framework. Currently, both the Krylov–Schur [46] method and Jacobi–Davidson [47] method are provided.In the implementation of the procedure, the entire response solution process consists mainly of the following four C ++ classes:(1)
Spectral Curve Management Class. This class is responsible for dealing with the spectral curves excited by a number of different fluctuating pressures, including the curve description, the curve dispersion, and so on. The curve discrete methods support the uniform discretization and the adaptive one, which is implemented according to the modal natural frequency and the modal damping ratio. These discrete data are called when the random vibration response is solving.(2)
Modal Participation Factor Class. The participation factors in different modal coordinates corresponding to each fluctuating pressure are calculated according to the modal shape and the load coefficient vector. These modal participation factors include the rigid-body parts and the elastic ones, which are stored in the form of an array and called when the random vibration response is solving.(3)
Random Vibration Response Class. This class is used to calculate the PSD response curve of each focused DOF and the equivalent RMS response in the whole field. The input of this class includes the modal participation coefficient, the frequency response function, each modal shape, and the excitation PSD curves.(4)
Output Class. This class is responsible for the output of the PSD response curves and the full-field equivalent RMS response. The corresponding calculation data are given according to the data format specified by the postprocessing software of the JAUMIN framework, which is convenient for curve drawing and full-field deformation image display.From the above introduction, the parallel implementation of the free-body random vibration algorithm in this paper, on the one hand, benefits from the parallel framework JUAMIN, which provides software and hardware bridges for parallel computing and support for parallel domain decomposition and load balancing processing, and provides commonly used mathematical solution library. On the other hand, the design of the PANDA platform is based on the established framework. It is mainly responsible for the design problems of the finite element calculation process and builds a common support structure for the solution of various types of problems. Such a design will be very conducive to the combination of software and supercomputers and can give full play to the parallel solution advantages of supercomputers. Developers can focus on the serial flow of the algorithm, and related parallel operations are handed over to the JAUMIN framework and the PANDA platform to complete. This is also the important difference and progress between the parallel implementation of this paper and the traditional commercial finite element software.


## 3.1. The JAUMIN Framework

The JAUMIN framework is mainly based on the hierarchical and modular design idea, using C++, MPI [42], and OpenMP mixed programming, which contains about 120,000 lines of code. The basic architecture of the JAUMIN framework is shown in Figure 1, which is divided into three levels.Figure 1
Architecture of the JAUMIN framework.The bottom layer is the high-performance computing support layer, composed of multiple modules for high-performance computing of an unstructured grid. In this layer, the relevant modules are divided into three sublayers. The first sublayer is a toolbox module, mainly to provide some basic tools, such as operating environment management, input database, memory management, and so on. The second layer includes three modules, in which the patch hierarchy structure and patch data are mainly used for the storage of finite element grid and field data and the data communication module is mainly used for data communication between different grid regions. The third sublayer includes grid adaptive module, load balancing module, and math operation module, which are used for adaptive refinement of finite element mesh, data balance processing of each process and related mathematical operation, and so on.The middle layer is the numerical general layer, including some modules for general numerical algorithms and application toolboxes, mathematical solvers for matrix-vector operations, time integrators, mesh generator, multi-physics coupling, and so on.The top layer is the applied interface layer, mainly providing the basic program interface for different applications. JAUMIN provides C++ interfaces for parallel programming of various application software. Based on this layer, users can write their own application subroutines based on the parallel data structure of the JAUMIN framework. These subroutines together constitute the application program of JAUMIN.In order to adapt to the large-scale parallel computing environment, a two-level partitioning strategy is adopted in the JAUMIN framework, shown as in Figure2. First, the finite element model is partitioned into submeshes (such as P0,P1 in Figure 2) considering load balance; each CPU processor is assigned with a submesh. Second, the submesh is again partitioned into a group of smaller submeshes, named “patch.” The patch is composed of some mesh elements and is the basic object for memory management and numerical computation in JAUMIN. Each patch includes some basic information such as patch geometry, patch topology, and patch data, which constitutes the core data structures of the framework.Figure 2
A two-level partitioning strategy in the JAUMIN framework.JAUMIN integrates a universal unstructured-mesh data model, a pattern-based parallel computing model, and a component-based programming model. These features enable JAUMIN to effectively shield parallel implementation details and significantly reduce programming complexity, thus serving as a bridge between supercomputers and applications. JAUMIN facilitates the development of large-scale parallel unstructured mesh applications and supports applications to scale to tens of thousands of CPU cores. Recently, JAUMIN has been widely used in a series of applications such as large deformation elastoplastic hydrodynamics, material damage and fracture, large complex structural impact dynamics, structural static and modal vibration analysis, and so on. For a more detailed introduction to the JAUMIN framework, we recommend literature [26].


## 3.2. The PANDA Platform

The PANDA platform is a large-scale finite element parallel computing platform based on the data structure of the JUMIN framework, which is mainly for large-scale finite element calculation of complex engineering equipment. At present, the PANDA platform includes static, modal and vibration, impact dynamics, and other types of finite element analysis, mainly coded with C++, including codes more than 100,000 lines. For the complex engineering structure, the pre- and postprocessing software are combined with the JAUMIN framework to complete the finite element modeling and grid domain decomposition, and the PANDA platform is responsible for the parallel matrix assembly and solution according to the structural material properties and external load conditions. The basic structure of the PANDA platform and its relationship with the JAUMIN framework are shown in Figure3.Figure 3
The PANDA platform and its relationship with the JAUMIN framework.The free-body random vibration analysis mentioned in this paper is taken as an analysis module in the PANDA-StaVib software, as shown in Figure4. In addition to the static analysis [43], the analysis software PANDA-StaVib includes some structural dynamic modules such as modal analysis, single- or multi-point harmonic response analysis [44], seismic response spectrum analysis, single- or multi-point basic excitation random vibration analysis [45], free-body random vibration analysis with multi-point fluctuating pressures, and so on, covering most of dynamic frequency domain analysis types. Although the above dynamic analysis modules are relatively independent, some common features such as the calculation of modal participation factor, description of input spectral curves, and modal analysis interfaces for modal superposition can be shared with each other, which is also compatible with the large-scale program design principle of “extracting common features, stripping personality,” reducing unnecessary code duplications.Figure 4
Analysis modules in PANDA-StaVib software.The PANDA platform mainly uses the finite method for numerical calculation and also includes the main hierarchical structure such as the finite element support layer, the finite element common layer, and the application personality layer. Among them, the finite element support layer provides data support services for the software, including the registration of JAUMIN variables, model data management, degree of freedom mapping, point data output, and toolboxes commonly used in matrix-vector operations; the finite element common layer is useful for structural mechanics finite element calculations where the common part is extracted, mainly including elements, materials, field data, constraints, loads and the corresponding solution process of each analysis type; and the finite element personality layer establishes a personality analysis module for different structural analysis types, as mentioned in this paper, random vibration analysis is one of the modules.


## 3.3. Parallel Implementation of Free-Body Random Vibration Analysis

The basic flow of the free-body random vibration analysis can be described in Figure5. In virtue of domain decomposition and load balancing management of the JAUMIN framework, the data information for each domain of finite element models is evenly distributed to each processor. In each processor, the PANDA platform combines the model parameter information to generate the mass matrix M, the stiffness matrix K, and the load coefficient vector fl. The information of matrices and vectors stored in each processor is parallel distributed, which provides the basis of data parallelism for the subsequent random vibration analysis of the free body.Figure 5
The basic calculation process of free-body random vibration analysis.The natural frequency and modal shapes, including all rigid body modes and partial low order elastic modes, are obtained with the free-body modal analysis in the PANDA platform. Using these data, the frequency response function for different discrete frequencies and the modal participation factors for different modes can be calculated and provided for the random vibration response. The spectral curve manager is used for the auto-PSD and cross-PSD curve management among different fluctuating pressures to provide the load input parameters for the calculation process. Based on these inputs, the random vibration response analysis is performed on each degree of freedom using the algorithm designed in Section2.4. Finally, the results are directed to the output stream, and the equivalent RMS response for all fields as well as the auto-PSD response curves for the concerned points are outputted.In the above solution process, the main parallel stages include the model domain decomposition and the matrix assembly, modal analysis, modal superposition, vibration response calculation phase, and so on. The modal analysis is mainly carried out by the eigenvalue solution library of the JAUMIN framework. Currently, both the Krylov–Schur [46] method and Jacobi–Davidson [47] method are provided.In the implementation of the procedure, the entire response solution process consists mainly of the following four C ++ classes:(1)
Spectral Curve Management Class. This class is responsible for dealing with the spectral curves excited by a number of different fluctuating pressures, including the curve description, the curve dispersion, and so on. The curve discrete methods support the uniform discretization and the adaptive one, which is implemented according to the modal natural frequency and the modal damping ratio. These discrete data are called when the random vibration response is solving.(2)
Modal Participation Factor Class. The participation factors in different modal coordinates corresponding to each fluctuating pressure are calculated according to the modal shape and the load coefficient vector. These modal participation factors include the rigid-body parts and the elastic ones, which are stored in the form of an array and called when the random vibration response is solving.(3)
Random Vibration Response Class. This class is used to calculate the PSD response curve of each focused DOF and the equivalent RMS response in the whole field. The input of this class includes the modal participation coefficient, the frequency response function, each modal shape, and the excitation PSD curves.(4)
Output Class. This class is responsible for the output of the PSD response curves and the full-field equivalent RMS response. The corresponding calculation data are given according to the data format specified by the postprocessing software of the JAUMIN framework, which is convenient for curve drawing and full-field deformation image display.From the above introduction, the parallel implementation of the free-body random vibration algorithm in this paper, on the one hand, benefits from the parallel framework JUAMIN, which provides software and hardware bridges for parallel computing and support for parallel domain decomposition and load balancing processing, and provides commonly used mathematical solution library. On the other hand, the design of the PANDA platform is based on the established framework. It is mainly responsible for the design problems of the finite element calculation process and builds a common support structure for the solution of various types of problems. Such a design will be very conducive to the combination of software and supercomputers and can give full play to the parallel solution advantages of supercomputers. Developers can focus on the serial flow of the algorithm, and related parallel operations are handed over to the JAUMIN framework and the PANDA platform to complete. This is also the important difference and progress between the parallel implementation of this paper and the traditional commercial finite element software.


## 4. Numerical Examples

In the previous sections, the algorithm theories and parallel implementation for the random vibration of free bodies are presented, thereby forming a relatively complete multi-point random vibration analysis module. Then, two typical examples will be employed to verify the correctness and parallel scalability of the designed algorithm and parallel module.


### 4.1. Example of the Sector Body

An irregular sector body is taken as a benchmark example of our software, whose finite element model is shown in Figure6. The length in the Y direction at the left end of the model is 1.75 m; that in X direction at the right end is 1.25 m; and the thickness in Z direction is 0.6 m. The structure is in an unconstrained free state, with two correlated pressure loads acting at both ends. The material of the sector body is steel, with an elastic modulus of 210 GPa, the density of 7,800 kg/m3, and Poisson’s ratio of 0.3. The model is meshed with hexahedral elements, containing 150 elements and 291 nodes.Figure 6
FEM of the sector body.The random vibration analysis frequency is set from 1 Hz to 1,000 Hz. The auto-PSDs at the two ends and the cross-PSD between the two loads are shown in Figure7. All PSDs are assumed as the white noise spectra. The auto spectrum amplitude at the left end is 4,000 Pa2/Hz; the auto spectrum amplitude at the right end is 1,000 Pa2/Hz; and the real part of the cross spectrum between the left end and the right end is 700 Pa2/Hz, while the imaginary part is 1,200 Pa2/Hz. The solution of the modal analysis is focused on the first 20 orders, and the ☆ marked node in the model is used as the calculation point to compare with the commercial finite element software. The displacement PSD curves of the ☆ marked node, as well as the RMS responses of the whole field, are calculated.Figure 7
The PSD curves for excitation.The first 20-order modal frequencies calculated by PANDA-StaVib and the commercial finite element software are shown in Table1. We can see that the results are almost the same, with a maximum relative error of 0.1%. Especially, the six rigid-body modal frequencies calculated by PANDA-StaVib can be accurate to zero in virtue of the spectral transform technology, while those calculated by the commercial finite element software just approximate to zero. The latter one contains certain computer floating-point error information.Table 1
Comparison of modal frequencies between PANDA-StaVib and the commercial software.
Modal ordersCommercial software (Hz)PADNA-StaVib (Hz)Relative error (%)12.64E − 050—22.42E − 050—37.82E − 060—42.45E − 060—51.62E − 050—62.46E − 050—7220.69220.790.058309.28309.400.049434.69434.310.0910486.94487.240.0611618.83619.100.0412732.75732.580.0213748.80749.580.1014799.62800.060.0615816.32815.610.0916952.50952.690.02171,031.901,032.700.08181,099.401,099.800.04191,139.101,138.700.04201,161.101,160.500.05The PSD curves in theY direction of node marked with ☆ calculated by PANDA-StaVib and the commercial finite element software are shown in Figures 8–10. The calculation results indicate that the rigid-body response, the elastic response, and the overall response calculated by PANDA-StaVib are perfectly agreed with the results of the commercial finite element software, which verifies the correctness of response calculation results of PANDA-StaVib. The PSD curves for the three different types of response shown in Figure 11 indicate that the rigid-body modes have an obvious influence on the response of the structure in the low-frequency areas, especially near the zero frequency. The above effect is gradually weakened with the increase of the frequency. Thus, the coupling effect of rigid-body mode should be considered in the random vibration analysis of free body. Figure 12 gives the comparison of equivalent RMS displacement fringes between PANDA-StaVib and the commercial finite element software, which are also consistent with each other.Figure 8
The PSD curves for rigid-body response.Figure 9
The PSD curves for elastic response.Figure 10
The PSD curves for total response.Figure 11
The comparison of PSD responses.Figure 12
Comparison of RMS displacement fringes.
(a)(b)


### 4.2. The Example of Reentry Vehicle

As a typical application of reentry vehicle, this example focuses on the PSD response and the full-field RMS response of a certain vehicle under the action of multi-point fluctuating pressure loads. Figure13 shows the finite element model of the simplified vehicle and the distributions of the fluctuating pressure, where the red arrows indicates the direction of the fluctuating pressure and the different colors on the outer surface indicate different action areas of the fluctuating pressures. Starting from the front section of the aircraft, the pressure surface is divided into five annular regions, which are identified as 1, 2, 3, 4, and 5. Each region suffers different pressure loads and the fluctuating pressures between the adjacent regions are cross-correlated. Two materials are used for modeling the vehicle. The external shell is steel with the elastic modulus of 210 GPa, the density of 7,800 kg/m3, and Poisson’s ratio of 0.3, while the internal counterweight is aluminum with the elastic modulus of 70 GPa, the density of 2,700 kg/m3, and Poisson’s ratio of 0.3.Figure 13
The model and load distribution of a reentry vehicle.
(a)(b)The description of fluctuating pressures includes the auto-PSD curves of five annular regions and the cross-PSD curves of the adjacent part. The analysis frequencies vary from 1 Hz to 5,000 Hz. The load spectrum curves are shown in Figure14, where AutoPSDx stands for the auto PSD of region x, CrossPSD12-Re stands for the real part of cross-PSD between region 1 and region 2, CrossPSD12-Im stands for the imaginary part of cross-PSD between region 1 and region 2, CrossPSDother-Re stands for the real part of cross-PSD between other adjacent regions, and CrossPSDother-Im stands for the imaginary part of cross-PSD between other adjacent regions. The amplitudes of AutoPSD1 till AutoPSD5 from the front to the end are 1.2e7 pa2/Hz, 8e6 pa2/Hz, 6e6 pa2/Hz, 4e6 pa2/Hz, and 2e6 pa2/Hz, respectively. The real and imaginary parts of the cross-PSD for the first two adjacent regions are 1e6 pa2/Hz and 2e6 pa2/Hz, respectively. The real and imaginary parts for the other adjacent regions are all 1e6 pa2/Hz. The entire vehicle is in an unconstrained free state and the first 100 modes (including six rigid-body modes) are calculated in the free-body modal analysis.Figure 14
The PSD curves of fluctuating pressures.The initial model is meshed with hexahedral elements, including 13,604 nodes and 40,812 DOFs. TheX-direction PSD response curves of some chosen point are shown in Figure 15, where rigid-body response, elastic response, and the total coupled response of the chosen point are compared in the same coordinate system. It can be seen that the overall displacement response of the structure is mainly affected by the rigid-body part with the low frequencies beneath 1,000 Hz, while that is mainly affected by the elastic part with the frequencies above 1,000 Hz.Figure 15
The PSD curves of fluctuating pressures.The equivalent RMS displacement fringes inX-direction are shown in Figure 16. The calculation time statistics for different stages with a single processor are shown in Table 2. The results indicate that the computational time for the random vibration analysis is only about 10% of the whole computational time. The most time-consuming is the modal analysis phase, accounting for more than half of the overall calculation time. Generally, the larger the computational scale is, the longer the computational time is. In the case of the random vibration calculations over hundreds of millions of DOFs, the modal analysis often takes more than 95% of total computational time.Figure 16
The RMS displacement fringes of reentry vehicle.
(a)(b)(c)Table 2
Computing time statistics for different solution stages.
Response typesTotal time (s)Matrix assembly (s)Modal analysis (s)Vibration response (s)Others (s)Rigid-body response46.973.7134.713.674.88Elastic response51.433.9435.86.55.19Total response51.823.9535.866.945.07In order to verify the parallel scalability of our developed module in PANDA-StaVib, a large-scale parallel test is carried out on the Tianhe supercomputer. The initial test model of the vehicle is 2.4 million DOFs. With the adaptive mesh refinement function of JAUMIN framework, the initial model is refined three times. We get the 18 million DOFs, 143 million DOFs, and 1.045 billion DOFs by the first refinement, the second refinement, and the third refinement, respectively. Figure17 gives the speedup ratio curves for four different computational scales from 256 CPU processors to 2048 CPU processors (for the 1.045 billion DOFs, the calculation starts from 512 processors as the memory limitation for 256 ones). Table 3 shows the parallel computing time and the maximum memory cost for different scales with 2,048 CPU processors.Figure 17
Speedup curves for different DOFs of reentry vehicle.Table 3
Parallel computing time and maximum memory cost for different DOFs of reentry vehicle with 2,048 CPU processors.
Computing DOFsComputing time (s)Maximum memory cost (GB)2.4 million3997518 million1,808650143 million14,7604,2001,045 million98,64021,000It can be seen in Figure17 and Table 3 that the larger the calculation scale is, the better the parallel speedup performance of the PANDA-StaVib software behaves. For the 1.045 billion DOFs, the parallel speedup ratio curve is close to the ideal speedup ratio because the amount of parallel data communication among processors relative to the amount of computation in each processor is significantly reduced with the increase of the computing scale. It takes about 27.4 hours for the calculation of 1.045 billion DOFs in 2,048 processors, which pose a great challenge in the random vibration analysis of free body so far.


## 4.1. Example of the Sector Body

An irregular sector body is taken as a benchmark example of our software, whose finite element model is shown in Figure6. The length in the Y direction at the left end of the model is 1.75 m; that in X direction at the right end is 1.25 m; and the thickness in Z direction is 0.6 m. The structure is in an unconstrained free state, with two correlated pressure loads acting at both ends. The material of the sector body is steel, with an elastic modulus of 210 GPa, the density of 7,800 kg/m3, and Poisson’s ratio of 0.3. The model is meshed with hexahedral elements, containing 150 elements and 291 nodes.Figure 6
FEM of the sector body.The random vibration analysis frequency is set from 1 Hz to 1,000 Hz. The auto-PSDs at the two ends and the cross-PSD between the two loads are shown in Figure7. All PSDs are assumed as the white noise spectra. The auto spectrum amplitude at the left end is 4,000 Pa2/Hz; the auto spectrum amplitude at the right end is 1,000 Pa2/Hz; and the real part of the cross spectrum between the left end and the right end is 700 Pa2/Hz, while the imaginary part is 1,200 Pa2/Hz. The solution of the modal analysis is focused on the first 20 orders, and the ☆ marked node in the model is used as the calculation point to compare with the commercial finite element software. The displacement PSD curves of the ☆ marked node, as well as the RMS responses of the whole field, are calculated.Figure 7
The PSD curves for excitation.The first 20-order modal frequencies calculated by PANDA-StaVib and the commercial finite element software are shown in Table1. We can see that the results are almost the same, with a maximum relative error of 0.1%. Especially, the six rigid-body modal frequencies calculated by PANDA-StaVib can be accurate to zero in virtue of the spectral transform technology, while those calculated by the commercial finite element software just approximate to zero. The latter one contains certain computer floating-point error information.Table 1
Comparison of modal frequencies between PANDA-StaVib and the commercial software.
Modal ordersCommercial software (Hz)PADNA-StaVib (Hz)Relative error (%)12.64E − 050—22.42E − 050—37.82E − 060—42.45E − 060—51.62E − 050—62.46E − 050—7220.69220.790.058309.28309.400.049434.69434.310.0910486.94487.240.0611618.83619.100.0412732.75732.580.0213748.80749.580.1014799.62800.060.0615816.32815.610.0916952.50952.690.02171,031.901,032.700.08181,099.401,099.800.04191,139.101,138.700.04201,161.101,160.500.05The PSD curves in theY direction of node marked with ☆ calculated by PANDA-StaVib and the commercial finite element software are shown in Figures 8–10. The calculation results indicate that the rigid-body response, the elastic response, and the overall response calculated by PANDA-StaVib are perfectly agreed with the results of the commercial finite element software, which verifies the correctness of response calculation results of PANDA-StaVib. The PSD curves for the three different types of response shown in Figure 11 indicate that the rigid-body modes have an obvious influence on the response of the structure in the low-frequency areas, especially near the zero frequency. The above effect is gradually weakened with the increase of the frequency. Thus, the coupling effect of rigid-body mode should be considered in the random vibration analysis of free body. Figure 12 gives the comparison of equivalent RMS displacement fringes between PANDA-StaVib and the commercial finite element software, which are also consistent with each other.Figure 8
The PSD curves for rigid-body response.Figure 9
The PSD curves for elastic response.Figure 10
The PSD curves for total response.Figure 11
The comparison of PSD responses.Figure 12
Comparison of RMS displacement fringes.
(a)(b)


## 4.2. The Example of Reentry Vehicle

As a typical application of reentry vehicle, this example focuses on the PSD response and the full-field RMS response of a certain vehicle under the action of multi-point fluctuating pressure loads. Figure13 shows the finite element model of the simplified vehicle and the distributions of the fluctuating pressure, where the red arrows indicates the direction of the fluctuating pressure and the different colors on the outer surface indicate different action areas of the fluctuating pressures. Starting from the front section of the aircraft, the pressure surface is divided into five annular regions, which are identified as 1, 2, 3, 4, and 5. Each region suffers different pressure loads and the fluctuating pressures between the adjacent regions are cross-correlated. Two materials are used for modeling the vehicle. The external shell is steel with the elastic modulus of 210 GPa, the density of 7,800 kg/m3, and Poisson’s ratio of 0.3, while the internal counterweight is aluminum with the elastic modulus of 70 GPa, the density of 2,700 kg/m3, and Poisson’s ratio of 0.3.Figure 13
The model and load distribution of a reentry vehicle.
(a)(b)The description of fluctuating pressures includes the auto-PSD curves of five annular regions and the cross-PSD curves of the adjacent part. The analysis frequencies vary from 1 Hz to 5,000 Hz. The load spectrum curves are shown in Figure14, where AutoPSDx stands for the auto PSD of region x, CrossPSD12-Re stands for the real part of cross-PSD between region 1 and region 2, CrossPSD12-Im stands for the imaginary part of cross-PSD between region 1 and region 2, CrossPSDother-Re stands for the real part of cross-PSD between other adjacent regions, and CrossPSDother-Im stands for the imaginary part of cross-PSD between other adjacent regions. The amplitudes of AutoPSD1 till AutoPSD5 from the front to the end are 1.2e7 pa2/Hz, 8e6 pa2/Hz, 6e6 pa2/Hz, 4e6 pa2/Hz, and 2e6 pa2/Hz, respectively. The real and imaginary parts of the cross-PSD for the first two adjacent regions are 1e6 pa2/Hz and 2e6 pa2/Hz, respectively. The real and imaginary parts for the other adjacent regions are all 1e6 pa2/Hz. The entire vehicle is in an unconstrained free state and the first 100 modes (including six rigid-body modes) are calculated in the free-body modal analysis.Figure 14
The PSD curves of fluctuating pressures.The initial model is meshed with hexahedral elements, including 13,604 nodes and 40,812 DOFs. TheX-direction PSD response curves of some chosen point are shown in Figure 15, where rigid-body response, elastic response, and the total coupled response of the chosen point are compared in the same coordinate system. It can be seen that the overall displacement response of the structure is mainly affected by the rigid-body part with the low frequencies beneath 1,000 Hz, while that is mainly affected by the elastic part with the frequencies above 1,000 Hz.Figure 15
The PSD curves of fluctuating pressures.The equivalent RMS displacement fringes inX-direction are shown in Figure 16. The calculation time statistics for different stages with a single processor are shown in Table 2. The results indicate that the computational time for the random vibration analysis is only about 10% of the whole computational time. The most time-consuming is the modal analysis phase, accounting for more than half of the overall calculation time. Generally, the larger the computational scale is, the longer the computational time is. In the case of the random vibration calculations over hundreds of millions of DOFs, the modal analysis often takes more than 95% of total computational time.Figure 16
The RMS displacement fringes of reentry vehicle.
(a)(b)(c)Table 2
Computing time statistics for different solution stages.
Response typesTotal time (s)Matrix assembly (s)Modal analysis (s)Vibration response (s)Others (s)Rigid-body response46.973.7134.713.674.88Elastic response51.433.9435.86.55.19Total response51.823.9535.866.945.07In order to verify the parallel scalability of our developed module in PANDA-StaVib, a large-scale parallel test is carried out on the Tianhe supercomputer. The initial test model of the vehicle is 2.4 million DOFs. With the adaptive mesh refinement function of JAUMIN framework, the initial model is refined three times. We get the 18 million DOFs, 143 million DOFs, and 1.045 billion DOFs by the first refinement, the second refinement, and the third refinement, respectively. Figure17 gives the speedup ratio curves for four different computational scales from 256 CPU processors to 2048 CPU processors (for the 1.045 billion DOFs, the calculation starts from 512 processors as the memory limitation for 256 ones). Table 3 shows the parallel computing time and the maximum memory cost for different scales with 2,048 CPU processors.Figure 17
Speedup curves for different DOFs of reentry vehicle.Table 3
Parallel computing time and maximum memory cost for different DOFs of reentry vehicle with 2,048 CPU processors.
Computing DOFsComputing time (s)Maximum memory cost (GB)2.4 million3997518 million1,808650143 million14,7604,2001,045 million98,64021,000It can be seen in Figure17 and Table 3 that the larger the calculation scale is, the better the parallel speedup performance of the PANDA-StaVib software behaves. For the 1.045 billion DOFs, the parallel speedup ratio curve is close to the ideal speedup ratio because the amount of parallel data communication among processors relative to the amount of computation in each processor is significantly reduced with the increase of the computing scale. It takes about 27.4 hours for the calculation of 1.045 billion DOFs in 2,048 processors, which pose a great challenge in the random vibration analysis of free body so far.


## 5. Summary and Outlook

In this paper, the random vibration problem caused by the air fluctuating pressure during the reentry of vehicles was studied. The theoretical analysis and algorithm design were carried out systematically, including the free-body modal analysis, the fluctuating pressure description, and the random vibration solution of the free body. The rigid-body response, the elastic response, and the overall coupling response were discussed. Combined with the JAUMIN framework and the PANDA platform, the parallel implementation was presented, and a large-scale parallel computing ability was formed, achieving the scale more than one billion DOFs. The research shows the following:(1)
The rigid-body motion of reentry vehicle has a significant effect on the overall response of the structure and even plays a leading role in the low-frequency domains. The displacement and velocity response near the zero frequency are divergent and tend to infinity. In the middle- and high-frequency domains, the elastic response gradually plays an important role. In the whole frequency band, the coupling effect between the rigid-body motion and elastic deformation is obvious. The rigid body motion cannot be ignored for the overall response of free body.(2)
The optimization of the RMS response solution in this paper greatly reduces the computational complexity of the parallel implementation. The calculation time of the random vibration is much less than that of the modal analysis; thus, the large-scale free-body random vibration analysis becomes possible. The maximum computing scale goes beyond the general commercial finite element software about 1 to 2 orders of magnitude.(3)
“Developing the application software based on a common parallel calculation framework” has become a popular idea for large-scale parallel software development, which has a great significance to enhance the parallel performance and the vitality of application software. The parallel computing ability obtained in this paper is partly due to the parallel data structure provided by the JAUMIN framework and the PANDA platform.At present, our self-developed PANDA software includes statics analysis, modal and vibration analyses, impact dynamics analysis, and other analysis types. Among them, modal and vibration analyses have basically covered the main frequency domain analysis types of structural dynamics. The relevant software has been used in complex reentry equipment such as weapons and aircraft. Next, we will gradually extend to the parallel calculation study of time-domain analysis and multi-field dynamics analysis to form fine numerical prediction ability for complex structures, and further applications in other structural fields will be carried out.



---

*Source: 6864973-2021-09-20.xml*