# A Component-Based Parametric Reduced-Order Modeling Method Combined with Substructural Matrix Interpolation and Automatic Sampling

Authors: Ying Liu; Hongguang Li; Yun Li; Huanyu Du
Journal: Shock and Vibration
 (2019)

---

## Abstract

An efficient parametric reduced-order modeling method combined with substructural matrix interpolation and automatic sampling procedure is proposed. This approach is based on the fixed-interface Craig-Bampton component mode synthesis method (CMS). The novel parametric reduced-order models (PROMs) are developed by interpolating substructural reduced-order matrices. To guarantee the compatibility of the coordinates, we develop a three-step adjustment procedure by reducing the local interface degrees of freedom (DOFs) and performing congruence transformation for the normal modes and interface reduced basis, respectively. In addition, an automatic sampling process is also introduced to dynamically fulfill the predefined error limits. It proceeds by first exploring the parameter space and identifying the sampling points with maximum error indicators for all the parameter-dependent substructures. The exact error of the assembled model at the optimal parameter point is subsequently calculated to determine whether the automatic sampling procedure reaches a desired error tolerance. The proposed framework is then applied to the moving coil of electrical-dynamic shaker to illustrate the advantage and validity. The results indicate that this new approach can significantly reduce both the offline database construction time and online calculation time. Besides, the automatic procedure can sample the parameter space efficiently and fulfill the stopping criterion dynamically with assurance of the resulting PROM accuracy.


---

## Body

## 1. Introduction

The finite element method (FEM) has been commonly used in a number of serious, complicated engineering problems. FE models with complex geometry and coupled physical effects often yield high-dimensional ordinary differential equations (ODE) to provide the needed accuracy. However, the analysis of full-order models (FOMs) often requires very extensive computational resources. To alleviate this problem, reduced-order models (ROMs) have been presented and widely used for many large-scale engineering applications. At this moment, research studies of reduced-order modeling methods of linear models have matured [1–3]. Also, these techniques have already been applied in some areas like optimization and control [4, 5], uncertainty quantification [6], and parameter identification [7].However, systems are usually parameter-dependent, and the existing reduced-order modeling methods are not applicable for the models with parametric variations, while in some areas like optimization, system models need repeated calculation over a range of parameter values, which restricts the application range of these methods. This has promoted the study of parametric reduced-order modeling technique. This approach is aimed to generate parametric reduced-order models (PROMs) that approximate the original full-order dynamical system with high fidelity over the range of parameters.Various parametric reduced-order modeling techniques have been developed over the last two decades [8]. One way to acquire the new reduced system is building the global basis matrices over the parameter space by the singular value decomposition (SVD) or a rank-revealing QR factorization and multiplying the full-order matrices by the global basis matrices [9]. By contrast, one can also interpolate the local reduced bases at discrete parameter points to obtain the basis matrix for a new parameter value. Direct interpolation may lead to nonorthogonality of the new basis matrix and proves invalid. A better method was proposed by Amsallem and Farhat [10] to interpolate the local bases on a tangent space of the Stiefel manifold. It is verified that the matrices achieved by this method can preserve the desired properties like orthogonality. Another kind of method is to interpolate the local reduced model matrices by using matrix manifold interpolation or direct interpolation [11], but first of all, these methods should perform a congruence transformation of the local reduced matrices in order that they are in the same generalized coordinate system. Then, Panzer et al. [12] used direct Lagrange interpolation and obtained the reduced matrices at the parameter. Different from that, Amsallem et al. [13] projected all the matrices onto a tangent space. After Lagrange interpolation, the new matrices were projected back. Both of the two ways prove effective in achieving the new reduced matrices. Yet another option for acquiring the new reduced system is to interpolate the transfer functions of local reduced models. Baur and Benner [14, 15] took advantage of local reduced transfer functions at the discrete parameter points and obtained the reduced transfer function at the given parameter value by Lagrange interpolation. However, the necessity of explicitly realizing the transfer function makes this method time consuming and not universal.All the parametric reduced-order modeling methods mentioned above are based on the offline/online decomposition strategy. In the offline stage, ROMs for different parameter values are generated and adjusted, and the database is thereby built. However, as the dimension of parameters increases, the corresponding sampling points increase exponentially and this approach becomes dramatically expensive, especially for large-scale FE models. Also, it is not necessary to repeatedly reduce the entire FE model when only partial structures change with parameters. Therefore, it is a computationally cheap way to partition the entire model into several substructures and separate the parts with parametric variations. Among the various substructuring techniques [16–19], the fixed-interface component mode synthesis method (CMS) was developed by Craig and Bampton [20] and improved by Rubin [21]. This method has been widely used due to its numerical stability and accuracy. By combining the parametric reduced-order method and CMS method, we can obtain the ROM database more efficiently in offline stage.Regarding the parametric reduced-order modeling methods combined with the CB-CMS, Park [22] and Zhang [23] developed a Taylor series-based method to approximate the component-level stiffness and mass matrices. As for the transformation matrices, they collected all the component modes at the sample parameter values and constructed a global basis vector. Hong [24] proposed a new transformation matrix and a novel parameterization technique to capture the element-level nonlinearities of brick-type finite elements. Lee and Cho [25] performed congruence transformation for the component modes and interpolated them in a common coordinate system. However, these methods are all aimed at the basis vectors and involve large-scale full-order multiplications between substructural stiffness or mass matrices and collected transformation vectors in the online phase, which is relatively time consuming.In order to further reduce the online calculation time, a novel component-based parametric reduced-order modeling method manipulating substructural reduced-order matrices is proposed in this research. A three-step adjustment procedure is developed and implemented for each substructure, and the transformed reduced-order component matrices are deduced, which is coordinate compatible for direct interpolation. Compared with previous approaches, this method is rather efficient since the calculation is conducted at the reduced-order level.Another drawback of previous methods is that these PROMs might lose their accuracy if the system matrices are polynomially expanded in terms of a parameter. To capture the polynomial dependence, more parameter samples are needed. The most fundamental and critical issue is the sampling strategy of parameter space. The distribution of the parameter sample points should enable the surrogate model to capture the underlying parametric dependence to a maximum extent. Some scholars [26–28] use the greedy approaches to iteratively solve the HDM and determine the configuration at which the current ROM leads to larger errors. However, computing the exact error requires the solution of the HDM, which is time consuming. Instead, error indicators evaluated with ROMs have been developed. To find the optimal sampling points associated with the largest error indicators, one major approach calculates the a posteriori error indicator of large discrete sets of candidate parameter values and locates the points with maximum errors [29, 30]. The other one solves the model-constrained optimization problem using gradient-based approaches and selects the solution as the sampling point [31]. Both of these two approaches can determine the optimal sampling location in parameter space.However, Most of the existing automatic sampling approaches are based on the entire structure. Few literatures applied the automatic sampling strategy to component-based parametric reduced-order modeling framework. The conventional uniform or random sampling procedure is unable to capture the parameter dependence well. Besides, the lack of error analysis provides no guaranty for the PROMs accuracy over the parameter space.Therefore, in this article, we have also developed an automatic parameter space sampling procedure for the novel component-based parametric reduced-order modeling method. In order to sample the parameter space automatically, an error indicator for rapid error evaluation is introduced. For each single substructure, we identify the optimal sample point by a direct exhaustive exploration of large discrete candidate parameter sets. All the components at the optimal parameter points are subsequently synthesized, and the exact errors between the directly computed CMS ROM and the interpolated CMS ROM are calculated to monitor the convergence rate and used as the termination criteria.This paper is organized as follows: Section2 gives an overview of the CB-CMS method. Section 3 describes the component-based parametric reduced-order modeling method, deduces the transformation matrices, presents the error analysis, and outlines the automatic sampling framework. In Section 4, the novel method is applied to the moving coil of electrical-dynamic shaker FE model to illustrate the advantage and validity of the proposed approach.


## 2. Reduced-Order Modeling with CB-CMS Method

For a damped structure, the FE equations of motion can be written in a matrix form as(1)M⋅x¨+C⋅x˙+K⋅x=F,where M,C,K∈ℝn×n are the mass, damping, and stiffness matrices of the system, respectively, X is the system response in physical coordinate, F is the force vector, and n is the degree of freedom (DOF) of the system.There are usually more than hundreds of thousands of DOFs for large FE models, which result in large system matrices and require vast amounts of computer memory and process time for solving. In order to reduce computational cost, the fixed-interface CMS method is employed to partition and generate substructures for its simplicity and computational stability.To apply the CMS method, the complex structure of interest is partitioned into substructures. The generation pass can be simplified if we omit the damping matrix and define master DOFs where loads are applied. The corresponding damping, force, and contact conditions are only included in the assembling stage. Therefore, for a typical undamped componenti, the equation of motion can be written as(2)Mi⋅x¨i+Ki⋅xi=Fi,i=1,2,…,N,where N denotes the number of subdomains; the physical DOFs of the stiffness and mass matrices can be partitioned into interior and interface DOFs as(3)Mi=MiOOMiOAMiAOMiAA,Ki=KiOOKiOAKiAOKiAA,where superscript A denotes the interface “active” DOFs including the excitation and response points and O denotes “omitted” DOFs of the interior. Then, the equation of motion for the ith substructure can be expressed as(4)MiOOMiOAMiAOMiAAx¨iOx¨iA+KiOOKiOAKiAOKiAAxiOxiA=FiOFiA.If the firstm low-order modes of the interior DOFs are adopted as the retained component modes, the following coordinate transformation matrix from physical coordinate to modal coordinate can be obtained by(5)xiOxiA=ΘixiOmxiA.The transformation matrixΘi consists of the low-order normal modes ΦiOm, and the constraint modes ΨiOA can be expressed as(6)Θi=ΦiOmΨiOA0iAmIiAA.Here the firstm low-order normal modes matrix ΦiOm can be obtained by solving the following eigen problem:(7)KiOO−λiαMiOOϕiO,α=0,α=1,2,…,m.And by assuming no external force on the interior DOFs, i.e.,FiO=0, the following constraint modes can be derived through the upper part in equation (4):(8)ΨiOA=−KiOO−1KiOA.Then, we can obtain the reduced stiffness and mass matrices of theith substructure by multiplying the transformation matrix to the original matrices:(9)M˜i=ΘiTMiΘi=M˜immM˜imAM˜iAmM˜iAA,K˜i=ΘiTKiΘi=K˜immK˜imAK˜iAmK˜iAA.The low-dimensional matrices are then synthesized after all the components are reduced. In general, the final reduced-order matrices of the entire assembled model can be characterized by one common set for the global interface coordinate and one unique set for each component’s internal coordinate. For illustration, we introduce a simple structure that consists of two substructures and one interface. The two components are denoted by using subscript 1 and 2, and the interface is designated by 3. Based on the interface compatibility, the following equation is tenable:(10)x1A=x2A=x3A.The displacement vectors of the two components and assembled model are then expressed as(11)x=x1Omx1Ax2Omx2A=I0000I0I000Ix1Omx2Omx3A=I0000I0I000Ix˜.The dynamic equation of the entire FE model can thereby be condensed as(12)M˜⋅x˜¨+C˜⋅x˜˙+K˜⋅x˜=F˜,where the structural damping is assumed in the form C˜=αM˜+βK˜, the force vector F˜ is generated based on the interface DOFs, and the assembled mass and stiffness matrices M˜ and K˜ are written as(13)M˜=M˜1mm0M˜1mA0M˜2mmM˜2mAM˜1AmM˜2AmM˜1AA+M˜2AA,(14)K˜=K˜1mm0K˜1mA0K˜2mmK˜2mAK˜1AmK˜2AmK˜1AA+K˜2AA.Conclusively, the entire calculation flow of fixed-interface CB-CMS method can be summarized as: partition the entire structure, generate the reduced-order substructures, synthesize the components, and calculate the reduced responses.


## 3. Proposed Parametric Reduced-Order Modeling Method

Global parametric reduced-order models (PROMs) [32] have been developed for fast reanalysis of structures with parametric variability in their properties. However, it is quite computationally expensive to calculate and collect system-level eigenvectors at multiple parameter sets. Therefore, Park [22] and Hong [24] developed component-based PROMs. Although these methods can reduce the offline computational cost to some extent, it is still relatively time consuming in the online phase, especially when the interface DOFs are relatively large scale, since these methods involve large-scale full-order online multiplications between substructural stiffness or mass matrices and collected transformation vectors. Another drawback is that these PROMs might lose their accuracy if the system matrices are polynomially expanded in terms of a parameter. To capture the polynomial dependence, more parameter samples are needed. It is a fundamental problem how the relative error of PROMs evolves as the number of parameter samples grows. These issues are addressed herein by developing a novel PROM using matrix interpolation and automatic sampling procedure. This new method performs the interpolation for the reduced-order substructural matrices after adjusting their coordinates and samples the parameter space dynamically to meet the predefined errors limits. In general, the procedure for this method can be summarized as constructing the ROMs of the substructures at initial parameter values, generating and transforming new components matrices at the next optimal sampling point to establish the compatible offline components database, interpolating the components ROMs for the relevant matrices at a new parameter value in the online phase, and then continuing the rest of the substructuring scheme.


### 3.1. Adjustment of the Offline Reduced-Order Component Matrices

To further reduce the online computational cost, we conduct the interpolation for the reduced-order substructural matrices. However, direct interpolation may lead to inaccuracy as these matrices are expressed in different modal coordinates and are incompatible. Thus, adjustment of the reduced-order component matrices should be performed at first. In this section, a three-step adjusting method is proposed. We first transform the substructural matrices according to the minimum Frobenius norm criteria between different normal modes. Then, combined with the local interface reduction method, the secondary transformation relation of the constraint DOFs is deduced. The reduced matrices corresponding to the interface DOFs are subsequently adjusted to guarantee that the resulting substructural matrices are compatible to interpolate.Consider an arbitraryith parameter-dependent substructure. The structural dynamical equation is described in parametric form as follows:(15)Mip⋅x¨i+Kip⋅xi=Fip,i=1,2,…,N.The uniform or random sampling method can be applied to sample the continuous parameter fieldp and generate a rich set of data with full coverage of the parameter space. An automatic sampling process is proposed in the following section.Assume now that one sampling method is utilized to obtain a discrete set of parameter samplespj for j = 1,…, P. The large-scale FE model is then generated at each corresponding discrete parameter value:(16)Mi,j⋅x¨i,j+Ki,j⋅xi,j=Fi,j,i=1,2,…,N,j=1,2,…,P.In order to perform the congruence transformation, the previous normal mode matrixΦiOm should satisfy the orthogonality condition ΦiOmTΦiOm=I which is guaranteed by a Gram–Schmidt procedure. The reduced-order matrices change accordingly.Assume that the reduced stiffness and mass matricesK˜i,j,M˜i,j can be rearranged as(17)M˜i,j=M˜i,jmmM˜i,jmAM˜i,jAmM˜i,jAA,K˜i,j=K˜i,jmmK˜i,jmAK˜i,jAmK˜i,jAA.The elements of the mass matrix can be further expressed as(18)M˜i,jmm=Φi,jOmTMiOOΦi,jOm,M˜i,jmA=Φi,jOmTMi,jOOΨi,jOA+Φi,jOmTMiOA,M˜i,jAm=Ψi,jOATMi,jOOΦi,jOm+MiAOΦi,jOm,M˜i,jAA=Ψi,jOATMi,jOOΨi,jOA+MiAOΨi,jOA+Ψi,jOATMiOA+MiAA,where Φi,jOm are the normal modes at jth parameter value and Ψi,jOA are the corresponding constraint modes. The elements of the stiffness matrix can be expressed in the same way.These different constraint modesΨi,jOA are described in physical coordinates such that they are compatible for interpolation, while the normal modes Φi,jOm are incompatible since they are expressed in each modal coordinate and represent different spaces. Therefore, the first step of transformation is aimed at adjusting the matrices related to the normal modes Φi,jOm.Assume that one reduced model is chosen as standard, and the corresponding matrices areM˜i,0,K˜i,0,Φi,0Om. In order to calculate the congruence transformation matrix for the jth reduced model, the following minimization problems are introduced:(19)Ri,j=minRi,j∈ΟmΦi,jOmRi,j−Φi,0OmF2,where Οm denotes a set of orthogonal matrices in ℝm×m, the subscript F represents the Frobenius norm, and Ri,j is a m-dimensional matrix.By expanding equation (19), the equivalent maximization problems is obtained as(20)Ri,j=maxRi,j∈ΟmtrRi,jTΦi,jOmTΦi,0Om.According to [33], equation (20) has an analytical solution as(21)Ri,j=Ui,jVi,jT,where Ui,j and Vi,j are calculated by the singular value decomposition as follows:(22)Φi,jOmTΦi,0Om=Ui,j∑i,jVi,jT.By multiplyingRi,j, all the normal modes Φi,jOm are transformed to coordinate compatible Φi,jOm∗ such that(23)Φi,jOm∗=Φi,jOmRi,j.The corresponding projection of the displacement vector can be expressed as(24)xi,jOxi,jA=Θi,jxi,jOmxi,jA=Zi,jx¯i,jOmxi,jA.The primary transformation matrixZi,j is deduced as follows:(25)Zi,j=Θi,j⋅Ri,j00I=Φi,jOmRi,jΨi,jOA0i,jAmIi,jAA.The transformed reduced stiffness and mass matricesK˜i,j∗andM˜i,j∗ can be rewritten as(26)M˜i,j∗=M˜i,jmm∗M˜i,jmA∗M˜i,jAm∗M˜i,jAA∗,(27)K˜i,j∗=K˜i,jmm∗K˜i,jmA∗K˜i,jAm∗K˜i,jAA∗.The elements of the mass matrix in equation (26) can be described as(28)M˜i,jmm∗=Φi,jOm∗TMiOOΦi,jOm∗=Ri,jTΦi,jOmTMi,jOOΦi,jOmRi,j=Ri,jTM˜i,jmmRi,j,M˜i,jmA∗=Φi,jOm∗TMi,jOOΨi,jOA+Φi,jOm∗TMiOA=Ri,jTΦi,jOmTMi,jOOΨi,jOA+Ri,jTΦi,jOmTMiOA=Ri,jTM˜i,jmA,M˜i,jAm∗=Ψi,jOATMi,jOOΦi,jOm∗+MiAOΦi,jOm∗=Ψi,jOATMi,jOOΦi,jOmRi,j+MiAOΦi,jOmRi,j=M˜i,jAmRi,j,M˜i,jAA∗=Ψi,jOATMi,jOOΨi,jOA+MiAOΨi,jOA+Ψi,jOATMiOA+MiAA=M˜i,jAA.Therefore, the entire reduced mass matrix is written as(29)M˜i,j∗=Ri,jTM˜i,jmmRi,jRi,jTM˜i,jmAM˜i,jAmRi,jM˜i,jAA.And the reduced stiffness matrix can be expressed in the same way such that(30)K˜i,j∗=Ri,jTK˜i,jmmRi,jRi,jTK˜i,jmAK˜i,jAmRi,jK˜i,jAA.Although the substructural matrices after the first-step transformation are compatible for interpolation, the size of these matrices could be relatively large scale, which is dominated by the constraint DOFs corresponding to theA partition in equations (29) and (30). To address this issue, Castanier et al. [34] suggested the characteristic constraint (CC) modes. However, this technique is performed at system level and thereby not applicable to our case. In the present study, the local interface reduction is implemented for the secondary transformation.The eigen problem on constraint DOFs of the mass and stiffness matrices of componenti at the jth parameter value is constructed by(31)K˜i,jAA∗θi,jO,α=λiαM˜i,jAA∗θi,jO,α,α=1,2,…,r.By selecting ther eigenvectors from the lowest eigenvalue, the reduced basis of the interface DOF is obtained as follows:(32)Ξi,jAr=θi,jA,1⋯θi,jA,r.Similar to the previous normal mode matrix, these reduced basis vectorsΞi,jAr should satisfy the orthogonality condition Ξi,jArTΞi,jAr=I. This is guaranteed by performing a Gram–Schmidt procedure for these vectors.The displacement vectors are further expressed as(33)xi,jOxi,jA=Θi,jxi,jOmxi,jA=Zi,jx¯i,jOmxi,jA=Γi,jx¯i,jOmxi,jAr,where Γi,j is the secondary transformation matrix and can be described as(34)Γi,j=Ζi,j⋅I00Ξi,jAr=Φi,jOmRi,jΨi,jOAΞi,jAr0i,jAmΞi,jAr.The reduced stiffness and mass matricesK˜i,j∗andM˜i,j∗ deduced by equations (29) and (30) can be rearranged as(35)M˜i,j∗∗=M˜i,jmm∗∗M˜i,jmA∗∗M˜i,jAm∗∗M˜i,jAA∗∗,K˜i,j∗∗=K˜i,jmm∗∗K˜i,jmA∗∗K˜i,jAm∗∗K˜i,jAA∗∗.Considering the reduced basisΞi,jAr, the elements of the mass matrix can be transformed as follows:(36)M˜i,jmm∗∗=Ri,jTΦi,jOmTMiOOΦi,jOmRi,j=Ri,jTM˜i,jmmRi,j,M˜i,jmA∗∗=Ri,jTΦi,jOmTMi,jOOΨi,jOAΞi,jAr+Ri,jTΦi,jOmTMiOAΞi,jAr=Ri,jTM˜i,jmAΞi,jAr,M˜i,jAm∗∗=Ξi,jArTΨi,jOATMi,jOOΦi,jOmRi,j+Ξi,jArTMiAOΦi,jOmRi,j=Ξi,jArTM˜i,jAmRi,j,M˜i,jAA∗∗=Ξi,jArTΨi,jOATMi,jOOΨi,jOAΞi,jAr+Ξi,jArTMiAOΨi,jOAΞi,jAr+Ψi,jOATMiOAΞi,jAr+Ξi,jArTMiAAΞi,jAr=Ξi,jArTM˜i,jAAΞi,jAr.Thus, the entire reduced mass matrix can be rewritten as(37)M˜i,j∗∗=Ri,jTM˜i,jmmRi,jRi,jTM˜i,jmAΞi,jArΞi,jArTM˜i,jAmRi,jΞi,jArTM˜i,jAAΞi,jAr.It should be noted that the projection has little effect on the results which is proved in [33]. And the reduced stiffness matrix can be expressed as(38)K˜i,j∗∗=Ri,jTK˜i,jmmRi,jRi,jTK˜i,jmAΞi,jArΞi,jArTK˜i,jAmRi,jΞi,jArTK˜i,jAAΞi,jAr.Closely resembling the first step, the reduced basis vectorsΞi,jAr represent different local modal spaces and are still incompatible. Therefore, the goal of the third step is to adjust, in a similar way to the first step, the substructural matrices related to these reduced basis vectors.Assuming that one reduced modelM˜i,0∗∗,K˜i,0∗∗,Ξi,0Ar is chosen as standard, the following minimization problems is derived:(39)Ti,j=minTi,j∈ΟrΞi,jArTi,j−Ξi,0ArF2.According to equations (20)–(23), the solution to equation (39) is(40)Ti,j=Ui,jVi,jT,where Ui,j and Vi,j are calculated by the singular value decomposition as follows:(41)Ξi,jArTΞi,0Ar=Ui,j∑i,jVi,jT.By multiplyingTi,j, all the reduced basis vectors Ξi,jAr are transformed to coordinate compatible Ξi,jAr∗∗∗ such that(42)Ξi,jAr∗∗∗=Ξi,jArTi,j.The relation between the displacement vectors in each coordinate system can be expressed as(43)xi,jOxi,jA=Θi,jxi,jOmxi,jA=Zi,jx¯i,jOmxi,jA=Γi,jx¯i,jOmxi,jAr=ϒi,jx¯i,jOmx¯i,jAr.The tertiary transformation matrixϒi,j is(44)ϒi,j=Γi,j⋅I00Ti,j=Φi,jOmRi,jΨi,jOAΞi,jArTi,j0i,jAmΞi,jArTi,j.The final reduced stiffness and mass matricesK˜i,j∗∗∗andM˜i,j∗∗∗ can be rearranged as(45)M˜i,j∗∗∗=M˜i,jmm∗∗∗M˜i,jmA∗∗∗M˜i,jAm∗∗∗M˜i,jAA∗∗∗,K˜i,j∗∗∗=K˜i,jmm∗∗∗K˜i,jmA∗∗∗K˜i,jAm∗∗∗K˜i,jAA∗∗∗.And the elements of the mass matrix can be deduced as(46)M˜i,jmm∗∗∗=Ri,jTΦi,jOmTMiOOΦi,jOmRi,j=Ri,jTM˜i,jmmRi,j,M˜i,jmA∗∗∗=Ri,jTΦi,jOmTMi,jOOΨi,jOAΞi,jArTi,j+Ri,jTΦi,jOmTMiOAΞi,jArTi,j=Ri,jTM˜i,jmAΞi,jArTi,j,M˜i,jAm∗∗∗=Ξi,jArTi,jTΨi,jOATMi,jOOΦi,jOmRi,j+Ξi,jArTi,jTMiAOΦi,jOmRi,j=Ξi,jArTi,jTM˜i,jAmRi,j,M˜i,jAA∗∗∗=Ξi,jArTi,jTΨi,jOATMi,jOOΨi,jOAΞi,jArTi,j+Ξi,jArTi,jTMiAOΨi,jOAΞi,jArTi,j+Ξi,jArTi,jTΨi,jOATMiOAΞi,jArTi,j+Ξi,jArTi,jTMiAAΞi,jArTi,j=Ξi,jArTi,jTM˜i,jAAΞi,jArTi,j.In the end, the relation between the final reduced stiffness and mass matrices and the original ones is expressed as(47)K˜i,j∗∗∗=Ri,jTK˜i,jmmRi,jRi,jTK˜i,jmAΞi,jArTi,jΞi,jArTi,jTK˜i,jAmRi,jΞi,jArTi,jTK˜i,jAAΞi,jArTi,j,M˜i,j∗∗∗=Ri,jTM˜i,jmmRi,jRi,jTM˜i,jmAΞi,jArTi,jΞi,jArTi,jTM˜i,jAmRi,jΞi,jArTi,jTM˜i,jAAΞi,jArTi,j.After the three-step adjustment of the coordinates, the final reduced-order substructural matrices are compatible for further online interpolation. The online calculation efficiency is thereby significantly improved thanks to the dramatic reduction of computational complexity, compared to the traditional methods.It should be noted that all the computations mentioned above are conducted at the subdomain level. Moreover, the transformation calculations are performed in the reduced coordinate and thus need little computational cost.


### 3.2. Online Calculation Using the Reduced-Order Component Matrices

As for the interpolation method, one can use gridded data interpolation for regular parameter samples or scattered data interpolation for irregular parameter values. Since the interpolation is performed for reduced matrices, the process is quite computational efficient; some high-order methods like polynomial or radial basis function interpolation approach are therefore available. The detailed description can be found in [30]. For simplicity, we utilize the Lagrange interpolation method and regular parameter samples in this research.Consider a reduced-order datasetp1,K˜i,1∗∗∗,M˜i,1∗∗∗,Ξi,1Ar∗∗∗⋯pP,K˜i,P∗∗∗,M˜i,P∗∗∗,Ξi,PAr∗∗∗. The matrices of the reduced system at new parameter values are interpolated by use of Lagrange interpolation function such that(48)K˜i∗∗∗p¯=∑j=1PLjp¯K˜i,j∗∗∗,(49)M˜i∗∗∗p¯=∑j=1PLjp¯M˜i,j∗∗∗,(50)ΞiAr∗∗∗p¯=∑j=1PLjp¯Ξi,jAr∗∗∗,where Ljp¯ are the Lagrange basis functions satisfying(51)∑j=1PLjp¯=1.These new matrices are then assembled in the online stage. The compatibility conditions for the components can be expressed almost identically to those described in equation (10) such that(52)x¯1Arp¯=x¯2Arp¯=x¯3Arp¯.The relation between the displacement vectors of components and assembly, together with the assembled dynamic equation, can also be expressed in the same way as equations (11)–(14). As the interface DOFs contain all the nodes of interest, we focus on the DOFs denoted by A.After the reduced response obtained, we project it back to the physical coordinate in order to calculate the actual response of the interface DOFs by(53)xAp¯=ΞiAr∗∗∗p¯x¯3Arp¯.Similar to this, the reduced modal shape vectors of the assembled model can also be projected by(54)ψAp¯=ΞiAr∗∗∗p¯ψ¯3Arp¯,where ψ¯3Arp¯ is the reduced modal shape vectors of the assembled system and ψAp¯ is the corresponding full-order ones of the interface DOFs in the physical coordinate.Compared to the traditional methods, the online interpolation of the reduced-order components matrices could reduce the computation complexity significantly, thus providing a better performance with cheaper online computational cost.


### 3.3. Automatic Sampling Procedure

In the previous sections, it is assumed that the parameter domainΩ is sampled for P discrete parameter values. In order to capture the parameter dependence to the greatest extent, multiple parameter samples are needed. The traditional component-based parametric reduced-order modeling approach cannot provide the error estimation. Without the error evaluation, there is no guaranty that the parameter space can be sampled effectively and that the PROMs can satisfy the predefined error limits. Therefore, we also propose an automatic sampling procedure for the component-based parametric reduced-order modeling method.For the automatic sampling procedure, the definition of errors is the most critical issue. In this research, the exact relative error between the directly computed CMS ROM and the interpolated CMS ROM has to be calculated and applied to monitor the convergence of the algorithm and terminate the procedure automatically. The exact relative errors of eigenvalues, modal shapes, and responses are defined as follows:(55)ε1,assembly=λp−λp¯λp,(56)ε2,assembly=maxψpA−ψAp¯maxψpA,(57)ε3,assembly=maxxpA−xAp¯maxxpA,where λ,ψ,x are eigenvalues, modal shapes, and responses of the assembled system described in equations (53) and (54).It should be noted that the subscriptp and the p¯ in brackets represent the results calculated using the exact model and the PROM at parameter value p, respectively. The symbol  denotes the absolute value.To calculate the exact relative error, we have to generate the corresponding FE model and calculate the directly computed CMS solution. In practice, however, this process is relatively time consuming. Therefore, the exact relative error is only used to monitor whether the current PROM satisfies the predefined error limits at the parameter with maximum errors, but it is still an open question where the point locates. Identification of the parameter points with maximum errors relies on the inexpensive computation of a posteriori error indicator associated with a ROM solution. Therefore, we utilize a computational efficient equivalent error to indicate where to place the sampling point.The exact errors are defined in global assembled system to evaluate the overall error level, while the error indicators are defined in local component systems and the search process is performed at a subdomain level.As the frequency can reflect components’ relative error, we use the eigenvalues calculated by spline interpolation as substitute to the exact results and define the relative error between frequencies calculated by spline interpolation and PROMs as error indicator:(58)ε4,i=λp¯i,spline−λp¯i,PMORλp¯i,spline.Here,λp¯i,spline is the calculated by spline interpolation for eigenvalues at discrete parameter values of the ith substructure, while λp¯i,PMOR is obtained by solving the following order-reduced eigenvalue problem:(59)K˜iAA∗∗∗p¯−λiαM˜iAA∗∗∗p¯ϕiαp¯=0,α=1,2,…,m,where K˜iAA∗∗∗p¯ and M˜iAA∗∗∗p¯ are the reduced interface DOFs partition of the stiffness and mass matrices for the ith substructure described in equations (48) and (49).Considering spline interpolation is quite an accurate interpolation method, this simplified error indicatorε4,i can almost represent the exact relative error between interpolated components ROM and the accurate FOM.The goal of the automatic parameter sampling approach is to estimate, at each iterationNiter, the parameter μmax maximizing the error indicator over the set ℂ for each component. These components with maximum error indicators are subsequently assembled to evaluate the exact relative errors.This approach proceeds in two phases. In the first phase, the error estimator is computed for each component over the set:(60)ℂ=μ1,…,μNc⊂Ω.Then, we build the following error pool to construct a surrogate model for the errorε4,i as a function of the parameter values:(61)E=μ1,ε4,i1,μ2,ε4,i2,…,μNc,ε4,iNc.At each iterationNiter of the automatic sampling process, the parameter value μmax is chosen as the maximizer:(62)μmax=argmaxμ∈ΩEμ.The parameter values are then used in the next phase to generate corresponding components. After these components are assembled, the overall exact errors are evaluated by equations (55)–(57) for the synthesized system.To sum up, the automatic sampling procedure can be summarized as evaluating the error indicator over a candidate parameter set for each component, generating new components at the parameters maximizing the posteriori error indicator, assembling the components, and calculating the overall exact errors for the synthesized system.


### 3.4. Entire Calculation Procedure

The entire automatic calculation procedure is schematically outlined in Figure1. In what follows, we explain the general idea behind this procedure step by step.Figure 1
The entire calculation procedure.Step 1.
Set the maximum number of iterationsNmax, maximum number of candidates Nc, and maximum error tolerance Tol.Step 2.
Run the following steps for the parameter-dependent substructures to construct the initial offline database. First, choose the start and end points of the parameter domain as initial parameter samples. Then, generate the full-order substructure matricesKi and Mi at these two parameter values. The corresponding normal modes and constraint modes are subsequently calculated through the CMS method. In the end, perform the three-step transformation for reduced-order matrices and build the initial compatible ROM database.Step 3.
Use the subroutine to locate the parameter points with maximum error indicators. First, randomly select a set ofNc candidate parameters ℂ=μ1,…,μNc⊂Ω for each component. Compute K˜i∗∗∗μk and M˜i∗∗∗μk by the Lagrange interpolation method for each candidate parameter. The corresponding λi,PMORμk,λi,splineμk,ε4,iμk are calculated and the parameter points of interest μmax are then calculated through equation (62). We search the parameter candidate sets for the optimal sample points for each substructure. The reader should note that the optimum returned by such a procedure does not need to be very accurate as sampling the FOM in the neighborhood of the parameter of maximum error will likely decrease the error at that parameter. After that, we perform the following automatic sampling and database updating process on the basis of the initial database. At first, we interpolate the previous ROMs database for new substructure matrices K˜i∗∗∗μmax and M˜i∗∗∗μmax at the parameter value μmax, assemble all the new ROMs, and compute the interpolated λμmax,ψAμmax,xAμmax. Then, we generate the new model K˜i,μmax,M˜i,μmax and Φi,μmaxOm,Ψi,μmaxOA at μmax and perform the three-step transformation procedure to obtain the new ROM K˜i,μmax∗∗∗,M˜i,μmax∗∗∗. These new reduced matrices are assembled to calculate the exact λμmax,ψμmaxA,xμmaxA. The errors of eigenvalues ε1,assemblyμmax, modal shape vectors ε2,assemblyμmax, and responses ε3,assemblyμmax are determined between the exact and interpolated models. If the final stopping criterion is not fulfilled, update the offline ROMs database and iterate step 3 again until the overall PROMs are satisfactory. If the final stopping criterion is met, then exit.
Compared with the global methods, this component-based parametric reduced-order modeling approach can significantly reduce the parameter sample points number fromΟPN to ΟP⋅N. Besides, in comparison with the classical component-based PROMs, the proposed one can substantially reduce the online computational complexity. Furthermore, the automatic sampling method can sample the parameter space according to the error indicator and terminate the procedure automatically when the maximum actual error becomes smaller than a specified tolerance. This method can place the sampling points where the error reaches the maximum. Therefore the approach can provide a better performance with smaller relative errors and cheaper offline and online computational cost.


## 3.1. Adjustment of the Offline Reduced-Order Component Matrices

To further reduce the online computational cost, we conduct the interpolation for the reduced-order substructural matrices. However, direct interpolation may lead to inaccuracy as these matrices are expressed in different modal coordinates and are incompatible. Thus, adjustment of the reduced-order component matrices should be performed at first. In this section, a three-step adjusting method is proposed. We first transform the substructural matrices according to the minimum Frobenius norm criteria between different normal modes. Then, combined with the local interface reduction method, the secondary transformation relation of the constraint DOFs is deduced. The reduced matrices corresponding to the interface DOFs are subsequently adjusted to guarantee that the resulting substructural matrices are compatible to interpolate.Consider an arbitraryith parameter-dependent substructure. The structural dynamical equation is described in parametric form as follows:(15)Mip⋅x¨i+Kip⋅xi=Fip,i=1,2,…,N.The uniform or random sampling method can be applied to sample the continuous parameter fieldp and generate a rich set of data with full coverage of the parameter space. An automatic sampling process is proposed in the following section.Assume now that one sampling method is utilized to obtain a discrete set of parameter samplespj for j = 1,…, P. The large-scale FE model is then generated at each corresponding discrete parameter value:(16)Mi,j⋅x¨i,j+Ki,j⋅xi,j=Fi,j,i=1,2,…,N,j=1,2,…,P.In order to perform the congruence transformation, the previous normal mode matrixΦiOm should satisfy the orthogonality condition ΦiOmTΦiOm=I which is guaranteed by a Gram–Schmidt procedure. The reduced-order matrices change accordingly.Assume that the reduced stiffness and mass matricesK˜i,j,M˜i,j can be rearranged as(17)M˜i,j=M˜i,jmmM˜i,jmAM˜i,jAmM˜i,jAA,K˜i,j=K˜i,jmmK˜i,jmAK˜i,jAmK˜i,jAA.The elements of the mass matrix can be further expressed as(18)M˜i,jmm=Φi,jOmTMiOOΦi,jOm,M˜i,jmA=Φi,jOmTMi,jOOΨi,jOA+Φi,jOmTMiOA,M˜i,jAm=Ψi,jOATMi,jOOΦi,jOm+MiAOΦi,jOm,M˜i,jAA=Ψi,jOATMi,jOOΨi,jOA+MiAOΨi,jOA+Ψi,jOATMiOA+MiAA,where Φi,jOm are the normal modes at jth parameter value and Ψi,jOA are the corresponding constraint modes. The elements of the stiffness matrix can be expressed in the same way.These different constraint modesΨi,jOA are described in physical coordinates such that they are compatible for interpolation, while the normal modes Φi,jOm are incompatible since they are expressed in each modal coordinate and represent different spaces. Therefore, the first step of transformation is aimed at adjusting the matrices related to the normal modes Φi,jOm.Assume that one reduced model is chosen as standard, and the corresponding matrices areM˜i,0,K˜i,0,Φi,0Om. In order to calculate the congruence transformation matrix for the jth reduced model, the following minimization problems are introduced:(19)Ri,j=minRi,j∈ΟmΦi,jOmRi,j−Φi,0OmF2,where Οm denotes a set of orthogonal matrices in ℝm×m, the subscript F represents the Frobenius norm, and Ri,j is a m-dimensional matrix.By expanding equation (19), the equivalent maximization problems is obtained as(20)Ri,j=maxRi,j∈ΟmtrRi,jTΦi,jOmTΦi,0Om.According to [33], equation (20) has an analytical solution as(21)Ri,j=Ui,jVi,jT,where Ui,j and Vi,j are calculated by the singular value decomposition as follows:(22)Φi,jOmTΦi,0Om=Ui,j∑i,jVi,jT.By multiplyingRi,j, all the normal modes Φi,jOm are transformed to coordinate compatible Φi,jOm∗ such that(23)Φi,jOm∗=Φi,jOmRi,j.The corresponding projection of the displacement vector can be expressed as(24)xi,jOxi,jA=Θi,jxi,jOmxi,jA=Zi,jx¯i,jOmxi,jA.The primary transformation matrixZi,j is deduced as follows:(25)Zi,j=Θi,j⋅Ri,j00I=Φi,jOmRi,jΨi,jOA0i,jAmIi,jAA.The transformed reduced stiffness and mass matricesK˜i,j∗andM˜i,j∗ can be rewritten as(26)M˜i,j∗=M˜i,jmm∗M˜i,jmA∗M˜i,jAm∗M˜i,jAA∗,(27)K˜i,j∗=K˜i,jmm∗K˜i,jmA∗K˜i,jAm∗K˜i,jAA∗.The elements of the mass matrix in equation (26) can be described as(28)M˜i,jmm∗=Φi,jOm∗TMiOOΦi,jOm∗=Ri,jTΦi,jOmTMi,jOOΦi,jOmRi,j=Ri,jTM˜i,jmmRi,j,M˜i,jmA∗=Φi,jOm∗TMi,jOOΨi,jOA+Φi,jOm∗TMiOA=Ri,jTΦi,jOmTMi,jOOΨi,jOA+Ri,jTΦi,jOmTMiOA=Ri,jTM˜i,jmA,M˜i,jAm∗=Ψi,jOATMi,jOOΦi,jOm∗+MiAOΦi,jOm∗=Ψi,jOATMi,jOOΦi,jOmRi,j+MiAOΦi,jOmRi,j=M˜i,jAmRi,j,M˜i,jAA∗=Ψi,jOATMi,jOOΨi,jOA+MiAOΨi,jOA+Ψi,jOATMiOA+MiAA=M˜i,jAA.Therefore, the entire reduced mass matrix is written as(29)M˜i,j∗=Ri,jTM˜i,jmmRi,jRi,jTM˜i,jmAM˜i,jAmRi,jM˜i,jAA.And the reduced stiffness matrix can be expressed in the same way such that(30)K˜i,j∗=Ri,jTK˜i,jmmRi,jRi,jTK˜i,jmAK˜i,jAmRi,jK˜i,jAA.Although the substructural matrices after the first-step transformation are compatible for interpolation, the size of these matrices could be relatively large scale, which is dominated by the constraint DOFs corresponding to theA partition in equations (29) and (30). To address this issue, Castanier et al. [34] suggested the characteristic constraint (CC) modes. However, this technique is performed at system level and thereby not applicable to our case. In the present study, the local interface reduction is implemented for the secondary transformation.The eigen problem on constraint DOFs of the mass and stiffness matrices of componenti at the jth parameter value is constructed by(31)K˜i,jAA∗θi,jO,α=λiαM˜i,jAA∗θi,jO,α,α=1,2,…,r.By selecting ther eigenvectors from the lowest eigenvalue, the reduced basis of the interface DOF is obtained as follows:(32)Ξi,jAr=θi,jA,1⋯θi,jA,r.Similar to the previous normal mode matrix, these reduced basis vectorsΞi,jAr should satisfy the orthogonality condition Ξi,jArTΞi,jAr=I. This is guaranteed by performing a Gram–Schmidt procedure for these vectors.The displacement vectors are further expressed as(33)xi,jOxi,jA=Θi,jxi,jOmxi,jA=Zi,jx¯i,jOmxi,jA=Γi,jx¯i,jOmxi,jAr,where Γi,j is the secondary transformation matrix and can be described as(34)Γi,j=Ζi,j⋅I00Ξi,jAr=Φi,jOmRi,jΨi,jOAΞi,jAr0i,jAmΞi,jAr.The reduced stiffness and mass matricesK˜i,j∗andM˜i,j∗ deduced by equations (29) and (30) can be rearranged as(35)M˜i,j∗∗=M˜i,jmm∗∗M˜i,jmA∗∗M˜i,jAm∗∗M˜i,jAA∗∗,K˜i,j∗∗=K˜i,jmm∗∗K˜i,jmA∗∗K˜i,jAm∗∗K˜i,jAA∗∗.Considering the reduced basisΞi,jAr, the elements of the mass matrix can be transformed as follows:(36)M˜i,jmm∗∗=Ri,jTΦi,jOmTMiOOΦi,jOmRi,j=Ri,jTM˜i,jmmRi,j,M˜i,jmA∗∗=Ri,jTΦi,jOmTMi,jOOΨi,jOAΞi,jAr+Ri,jTΦi,jOmTMiOAΞi,jAr=Ri,jTM˜i,jmAΞi,jAr,M˜i,jAm∗∗=Ξi,jArTΨi,jOATMi,jOOΦi,jOmRi,j+Ξi,jArTMiAOΦi,jOmRi,j=Ξi,jArTM˜i,jAmRi,j,M˜i,jAA∗∗=Ξi,jArTΨi,jOATMi,jOOΨi,jOAΞi,jAr+Ξi,jArTMiAOΨi,jOAΞi,jAr+Ψi,jOATMiOAΞi,jAr+Ξi,jArTMiAAΞi,jAr=Ξi,jArTM˜i,jAAΞi,jAr.Thus, the entire reduced mass matrix can be rewritten as(37)M˜i,j∗∗=Ri,jTM˜i,jmmRi,jRi,jTM˜i,jmAΞi,jArΞi,jArTM˜i,jAmRi,jΞi,jArTM˜i,jAAΞi,jAr.It should be noted that the projection has little effect on the results which is proved in [33]. And the reduced stiffness matrix can be expressed as(38)K˜i,j∗∗=Ri,jTK˜i,jmmRi,jRi,jTK˜i,jmAΞi,jArΞi,jArTK˜i,jAmRi,jΞi,jArTK˜i,jAAΞi,jAr.Closely resembling the first step, the reduced basis vectorsΞi,jAr represent different local modal spaces and are still incompatible. Therefore, the goal of the third step is to adjust, in a similar way to the first step, the substructural matrices related to these reduced basis vectors.Assuming that one reduced modelM˜i,0∗∗,K˜i,0∗∗,Ξi,0Ar is chosen as standard, the following minimization problems is derived:(39)Ti,j=minTi,j∈ΟrΞi,jArTi,j−Ξi,0ArF2.According to equations (20)–(23), the solution to equation (39) is(40)Ti,j=Ui,jVi,jT,where Ui,j and Vi,j are calculated by the singular value decomposition as follows:(41)Ξi,jArTΞi,0Ar=Ui,j∑i,jVi,jT.By multiplyingTi,j, all the reduced basis vectors Ξi,jAr are transformed to coordinate compatible Ξi,jAr∗∗∗ such that(42)Ξi,jAr∗∗∗=Ξi,jArTi,j.The relation between the displacement vectors in each coordinate system can be expressed as(43)xi,jOxi,jA=Θi,jxi,jOmxi,jA=Zi,jx¯i,jOmxi,jA=Γi,jx¯i,jOmxi,jAr=ϒi,jx¯i,jOmx¯i,jAr.The tertiary transformation matrixϒi,j is(44)ϒi,j=Γi,j⋅I00Ti,j=Φi,jOmRi,jΨi,jOAΞi,jArTi,j0i,jAmΞi,jArTi,j.The final reduced stiffness and mass matricesK˜i,j∗∗∗andM˜i,j∗∗∗ can be rearranged as(45)M˜i,j∗∗∗=M˜i,jmm∗∗∗M˜i,jmA∗∗∗M˜i,jAm∗∗∗M˜i,jAA∗∗∗,K˜i,j∗∗∗=K˜i,jmm∗∗∗K˜i,jmA∗∗∗K˜i,jAm∗∗∗K˜i,jAA∗∗∗.And the elements of the mass matrix can be deduced as(46)M˜i,jmm∗∗∗=Ri,jTΦi,jOmTMiOOΦi,jOmRi,j=Ri,jTM˜i,jmmRi,j,M˜i,jmA∗∗∗=Ri,jTΦi,jOmTMi,jOOΨi,jOAΞi,jArTi,j+Ri,jTΦi,jOmTMiOAΞi,jArTi,j=Ri,jTM˜i,jmAΞi,jArTi,j,M˜i,jAm∗∗∗=Ξi,jArTi,jTΨi,jOATMi,jOOΦi,jOmRi,j+Ξi,jArTi,jTMiAOΦi,jOmRi,j=Ξi,jArTi,jTM˜i,jAmRi,j,M˜i,jAA∗∗∗=Ξi,jArTi,jTΨi,jOATMi,jOOΨi,jOAΞi,jArTi,j+Ξi,jArTi,jTMiAOΨi,jOAΞi,jArTi,j+Ξi,jArTi,jTΨi,jOATMiOAΞi,jArTi,j+Ξi,jArTi,jTMiAAΞi,jArTi,j=Ξi,jArTi,jTM˜i,jAAΞi,jArTi,j.In the end, the relation between the final reduced stiffness and mass matrices and the original ones is expressed as(47)K˜i,j∗∗∗=Ri,jTK˜i,jmmRi,jRi,jTK˜i,jmAΞi,jArTi,jΞi,jArTi,jTK˜i,jAmRi,jΞi,jArTi,jTK˜i,jAAΞi,jArTi,j,M˜i,j∗∗∗=Ri,jTM˜i,jmmRi,jRi,jTM˜i,jmAΞi,jArTi,jΞi,jArTi,jTM˜i,jAmRi,jΞi,jArTi,jTM˜i,jAAΞi,jArTi,j.After the three-step adjustment of the coordinates, the final reduced-order substructural matrices are compatible for further online interpolation. The online calculation efficiency is thereby significantly improved thanks to the dramatic reduction of computational complexity, compared to the traditional methods.It should be noted that all the computations mentioned above are conducted at the subdomain level. Moreover, the transformation calculations are performed in the reduced coordinate and thus need little computational cost.


## 3.2. Online Calculation Using the Reduced-Order Component Matrices

As for the interpolation method, one can use gridded data interpolation for regular parameter samples or scattered data interpolation for irregular parameter values. Since the interpolation is performed for reduced matrices, the process is quite computational efficient; some high-order methods like polynomial or radial basis function interpolation approach are therefore available. The detailed description can be found in [30]. For simplicity, we utilize the Lagrange interpolation method and regular parameter samples in this research.Consider a reduced-order datasetp1,K˜i,1∗∗∗,M˜i,1∗∗∗,Ξi,1Ar∗∗∗⋯pP,K˜i,P∗∗∗,M˜i,P∗∗∗,Ξi,PAr∗∗∗. The matrices of the reduced system at new parameter values are interpolated by use of Lagrange interpolation function such that(48)K˜i∗∗∗p¯=∑j=1PLjp¯K˜i,j∗∗∗,(49)M˜i∗∗∗p¯=∑j=1PLjp¯M˜i,j∗∗∗,(50)ΞiAr∗∗∗p¯=∑j=1PLjp¯Ξi,jAr∗∗∗,where Ljp¯ are the Lagrange basis functions satisfying(51)∑j=1PLjp¯=1.These new matrices are then assembled in the online stage. The compatibility conditions for the components can be expressed almost identically to those described in equation (10) such that(52)x¯1Arp¯=x¯2Arp¯=x¯3Arp¯.The relation between the displacement vectors of components and assembly, together with the assembled dynamic equation, can also be expressed in the same way as equations (11)–(14). As the interface DOFs contain all the nodes of interest, we focus on the DOFs denoted by A.After the reduced response obtained, we project it back to the physical coordinate in order to calculate the actual response of the interface DOFs by(53)xAp¯=ΞiAr∗∗∗p¯x¯3Arp¯.Similar to this, the reduced modal shape vectors of the assembled model can also be projected by(54)ψAp¯=ΞiAr∗∗∗p¯ψ¯3Arp¯,where ψ¯3Arp¯ is the reduced modal shape vectors of the assembled system and ψAp¯ is the corresponding full-order ones of the interface DOFs in the physical coordinate.Compared to the traditional methods, the online interpolation of the reduced-order components matrices could reduce the computation complexity significantly, thus providing a better performance with cheaper online computational cost.


## 3.3. Automatic Sampling Procedure

In the previous sections, it is assumed that the parameter domainΩ is sampled for P discrete parameter values. In order to capture the parameter dependence to the greatest extent, multiple parameter samples are needed. The traditional component-based parametric reduced-order modeling approach cannot provide the error estimation. Without the error evaluation, there is no guaranty that the parameter space can be sampled effectively and that the PROMs can satisfy the predefined error limits. Therefore, we also propose an automatic sampling procedure for the component-based parametric reduced-order modeling method.For the automatic sampling procedure, the definition of errors is the most critical issue. In this research, the exact relative error between the directly computed CMS ROM and the interpolated CMS ROM has to be calculated and applied to monitor the convergence of the algorithm and terminate the procedure automatically. The exact relative errors of eigenvalues, modal shapes, and responses are defined as follows:(55)ε1,assembly=λp−λp¯λp,(56)ε2,assembly=maxψpA−ψAp¯maxψpA,(57)ε3,assembly=maxxpA−xAp¯maxxpA,where λ,ψ,x are eigenvalues, modal shapes, and responses of the assembled system described in equations (53) and (54).It should be noted that the subscriptp and the p¯ in brackets represent the results calculated using the exact model and the PROM at parameter value p, respectively. The symbol  denotes the absolute value.To calculate the exact relative error, we have to generate the corresponding FE model and calculate the directly computed CMS solution. In practice, however, this process is relatively time consuming. Therefore, the exact relative error is only used to monitor whether the current PROM satisfies the predefined error limits at the parameter with maximum errors, but it is still an open question where the point locates. Identification of the parameter points with maximum errors relies on the inexpensive computation of a posteriori error indicator associated with a ROM solution. Therefore, we utilize a computational efficient equivalent error to indicate where to place the sampling point.The exact errors are defined in global assembled system to evaluate the overall error level, while the error indicators are defined in local component systems and the search process is performed at a subdomain level.As the frequency can reflect components’ relative error, we use the eigenvalues calculated by spline interpolation as substitute to the exact results and define the relative error between frequencies calculated by spline interpolation and PROMs as error indicator:(58)ε4,i=λp¯i,spline−λp¯i,PMORλp¯i,spline.Here,λp¯i,spline is the calculated by spline interpolation for eigenvalues at discrete parameter values of the ith substructure, while λp¯i,PMOR is obtained by solving the following order-reduced eigenvalue problem:(59)K˜iAA∗∗∗p¯−λiαM˜iAA∗∗∗p¯ϕiαp¯=0,α=1,2,…,m,where K˜iAA∗∗∗p¯ and M˜iAA∗∗∗p¯ are the reduced interface DOFs partition of the stiffness and mass matrices for the ith substructure described in equations (48) and (49).Considering spline interpolation is quite an accurate interpolation method, this simplified error indicatorε4,i can almost represent the exact relative error between interpolated components ROM and the accurate FOM.The goal of the automatic parameter sampling approach is to estimate, at each iterationNiter, the parameter μmax maximizing the error indicator over the set ℂ for each component. These components with maximum error indicators are subsequently assembled to evaluate the exact relative errors.This approach proceeds in two phases. In the first phase, the error estimator is computed for each component over the set:(60)ℂ=μ1,…,μNc⊂Ω.Then, we build the following error pool to construct a surrogate model for the errorε4,i as a function of the parameter values:(61)E=μ1,ε4,i1,μ2,ε4,i2,…,μNc,ε4,iNc.At each iterationNiter of the automatic sampling process, the parameter value μmax is chosen as the maximizer:(62)μmax=argmaxμ∈ΩEμ.The parameter values are then used in the next phase to generate corresponding components. After these components are assembled, the overall exact errors are evaluated by equations (55)–(57) for the synthesized system.To sum up, the automatic sampling procedure can be summarized as evaluating the error indicator over a candidate parameter set for each component, generating new components at the parameters maximizing the posteriori error indicator, assembling the components, and calculating the overall exact errors for the synthesized system.


## 3.4. Entire Calculation Procedure

The entire automatic calculation procedure is schematically outlined in Figure1. In what follows, we explain the general idea behind this procedure step by step.Figure 1
The entire calculation procedure.Step 1.
Set the maximum number of iterationsNmax, maximum number of candidates Nc, and maximum error tolerance Tol.Step 2.
Run the following steps for the parameter-dependent substructures to construct the initial offline database. First, choose the start and end points of the parameter domain as initial parameter samples. Then, generate the full-order substructure matricesKi and Mi at these two parameter values. The corresponding normal modes and constraint modes are subsequently calculated through the CMS method. In the end, perform the three-step transformation for reduced-order matrices and build the initial compatible ROM database.Step 3.
Use the subroutine to locate the parameter points with maximum error indicators. First, randomly select a set ofNc candidate parameters ℂ=μ1,…,μNc⊂Ω for each component. Compute K˜i∗∗∗μk and M˜i∗∗∗μk by the Lagrange interpolation method for each candidate parameter. The corresponding λi,PMORμk,λi,splineμk,ε4,iμk are calculated and the parameter points of interest μmax are then calculated through equation (62). We search the parameter candidate sets for the optimal sample points for each substructure. The reader should note that the optimum returned by such a procedure does not need to be very accurate as sampling the FOM in the neighborhood of the parameter of maximum error will likely decrease the error at that parameter. After that, we perform the following automatic sampling and database updating process on the basis of the initial database. At first, we interpolate the previous ROMs database for new substructure matrices K˜i∗∗∗μmax and M˜i∗∗∗μmax at the parameter value μmax, assemble all the new ROMs, and compute the interpolated λμmax,ψAμmax,xAμmax. Then, we generate the new model K˜i,μmax,M˜i,μmax and Φi,μmaxOm,Ψi,μmaxOA at μmax and perform the three-step transformation procedure to obtain the new ROM K˜i,μmax∗∗∗,M˜i,μmax∗∗∗. These new reduced matrices are assembled to calculate the exact λμmax,ψμmaxA,xμmaxA. The errors of eigenvalues ε1,assemblyμmax, modal shape vectors ε2,assemblyμmax, and responses ε3,assemblyμmax are determined between the exact and interpolated models. If the final stopping criterion is not fulfilled, update the offline ROMs database and iterate step 3 again until the overall PROMs are satisfactory. If the final stopping criterion is met, then exit.
Compared with the global methods, this component-based parametric reduced-order modeling approach can significantly reduce the parameter sample points number fromΟPN to ΟP⋅N. Besides, in comparison with the classical component-based PROMs, the proposed one can substantially reduce the online computational complexity. Furthermore, the automatic sampling method can sample the parameter space according to the error indicator and terminate the procedure automatically when the maximum actual error becomes smaller than a specified tolerance. This method can place the sampling points where the error reaches the maximum. Therefore the approach can provide a better performance with smaller relative errors and cheaper offline and online computational cost.


## 4. Numerical Example

In this section, the novel PROM is illustrated with its application to the moving coil of electrical-dynamic shakers. Since the geometric model has perfect cyclic symmetry property, it can be simplified to one over sixteen of the entire model. The simplified moving coil model consists of a drive coil, a stiffener plate, and a table. At the position shown in Figure2(b), only axial displacement is permitted and the normal displacement at the symmetric boundary is zero. The model is excited by equivalent 1000 N surface force imposed on the interface of the drive coil and the stiffener plate. The measuring points are also marked in Figure 2.Figure 2
The moving coil of electrical-dynamic shakers. (a) The entire geometric model. (b) The simplified FE model.
(a)
(b)For simplicity, we partition the entire model into three substructures, and the subdomains along with the corresponding shape parameters are marked in Figure3. The material is linear elastic with E = 69 GPa, ν = 0.33, and ρ = 2.7e3 Kg/cm3, roughly corresponding to aluminum. The finite element model consists of 15608 SOLID elements and 21246 nodes. The size of the full-order model is 59487, and the number of interface DOF is 1140, which is 1.9% of the full DOFs.Figure 3
The interfaces and shape parameters of the substructures. (a) Substructure 1. (b) Substructure 2. (c) Substructure 3.
(a)
(b)
(c)Regarding the parameter space, each substructure contains one shape parameter. In this application, the thickness of stiffener platesT, changing from 10 mm to 32 mm, the length of drive coil L, changing from 242 mm to 342 mm, and the diameter of table D, changing from 920 mm to 1020 mm, span the three-dimensional parameter space. The original parameter configuration is T = 18 mm, L = 292 mm, and D = 970 mm.It should be noted that even though the shape parameters of the mesh vary, the underlying mesh topology remains the same. The preprocess was conducted in HYPERMESH to parameterize the mesh. To be specific, we first import the simplified geometry model into HYPERMESH and generate the original regular mesh. The HYPERMORPH module is subsequently utilized to parameterize the mesh without changing the underlying topology. Two FE models corresponding to the initial and final parameter sets are generated, respectively. Their node coordinates are then exported for use of generation of new FE model at arbitrary parameter values in ANSYS. The interface and measuring points are set as master nodes. The boundary condition and excitation condition shown in Figure2(b) are also initialized in HYPERMESH.The sampling scheme is automatically implemented in MATLAB including calling ANSYS to generate FE model and corresponding substructure matrices, importing the component matrices files, synthesizing the substructures, and conducting the automatic methods outlined in Figure1. All numerical calculations are performed on a Windows Machine with Intel Core (TM) i5-4210U CPU at 2.60 GHz and 8 GB of RAM.In this research, Latin hypercube sampling (LHS) is used to select the random candidates. The number of random candidates is set toNc = 20 for all substructures. The reduced-order basis applied in this approach is built by 10-order orthogonal eigenvectors after a Gram–Schmidt process. We set the exact errors of the assembled systems as the stopping criterion. Only components with the maximum error indicator are assembled for evaluation. The performance of this approach is assessed by monitoring the maximum errors of frequency, modal shapes, and responses between the directly computed CB-CMS ROM and interpolated CB-CMS ROM. The parameter sampling points for each component change from 2 to 7, and this number varies simultaneously for all components.Figure4 presents the evolution of exact maximum errors of the assembled system defined in equations (55)–(57). The results indicate that the maximum errors decrease significantly as the parameter sample points increase. It can also be observed that the responses has larger errors, compared with the other two variables. Therefore, it can be used as the stopping criteria of the algorithm.Figure 4
Error evolutions for the automatic sampling methods.The convergence performance of the proposed approach is also assessed. There is no guaranty that adding a parameter sample point will necessary lead to a monotonic decrease of maximum errors, but in practice, one can observe that the maximum error will decrease for most additional sample points due to the convergence property of the Lagrange interpolation method exploited in our approach.In order to illustrate the advantage of the proposed CMS-PROMs over the traditional global PROMs developed in [32], we compare the offline PROM database generation time and sampling points’ number.In the global parametric reduced-order modeling procedure, 7 uniformly distributed sample points are set for each dimension of the three-dimensional parameter space as same as the proposed automatic procedure. In ANSYS, it is relatively time consuming to change the nodes coordinate by use of NMODIF command. Therefore, the total time includes FE model generation and model reduction time.One can observe from Figure5 that the proposed procedure outperforms the global method thanks to its ability to partition the entire parameter space into separate ones with the help of CMS method. This new method can reduce significantly the sample points’ number from 73  =  343 to 7 ∗ 3 = 21. The offline model generation and reduction time decrease a lot accordingly. The computational benefit of the proposed approach can be apparently observed with a speedup ratio of about 19 times. In general, the proposed approach can improve the calculation efficiency and reduce the offline ROMs database construction time to a great extent.Figure 5
Comparison of offline database construction time between the global PROM and the novel CMS-PROM.In addition, this novel CMS-PROM surpasses not only the global PROM but also the classical CMS-PROM [24] in the aspects of online calculation time and accuracy. To validate the performance, we first compare the online modal analysis time for these two CMS-PROMs.It can be observed from Figure6 that the online calculation time of classical CMS-PROM is almost 120 times larger than its counterpart calculated using the proposed CMS-PROM. This is due to the fact that the classical method still involves time-consuming online multiplications between stiffness or mass matrices and collected transformation vectors. As the number of interface DOFs and parameter samples increases, the calculation time quickly grows. The proposed approach, by contrast, performs the online calculation only for the reduced-order component matrices, thus leading to cheaper online calculation cost.Figure 6
Comparison of online calculation time between the classical CMS-PROM and the proposed CMS-PROM.Furthermore, to assess the accuracy of these CMS-PROMs, we subsequently perform accuracy comparisons for first 5-order frequencies. The query parameter point isT = 30.31 mm, L = 335.22 mm, and D = 968.19 mm at which the error indicator reaches the maximum.As shown in Table1, the proposed approach leads to more accurate PROM than its classical counterpart. The results suggest that the proposed CMS-PROM can achieve a higher level of accuracy at about Ο10−6 than the classical CMS-PROM at about Ο10−3. Moreover, the expected values of the frequency errors confirm the results observed in Figure 4.Table 1
Comparison of nature frequencies between different CMS ROMs.
Algorithm
Mode
1
2
3
4
5
Direct (Hz)
4.7920
1440.0916
1552.6679
1747.8217
1806.7117
Classical (Hz)
4.7920
1441.5584
1559.7724
1748.2735
1806.9418
Relative error (%)
0
0.102
0.458
0.0258
0.0127
Proposed (Hz)
4.7920
1440.0927
1552.6695
1747.8352
1806.7189
Relative error (%)
0
7.638 × 10−5
1.031 × 10−4
7.724 × 10−4
3.985 × 10−4Then, we compare the harmonic response for the direct CMS ROM and interpolated CMS ROM. The frequency range is from 1 Hz to 1900 Hz with the frequency interval 0.1 Hz. Rayleigh damping coefficient is applied as alpha = 2.5, beta = 9 ✕ 10−5, and the detecting point locates at the center of top surface as shown in Figure 2(b).To display the distinct comparison, different scales are adopted for thex-axis. For example, the results between 0 and 1500 Hz are plotted using a base 10 logarithmic scale while for the range of 1500∼1900 Hz, linear scales are used.From Figure7, one can observe that the harmonic response predicted by the proposed CMS-PROM is almost indistinguishable from its counterpart calculated by direct CMS ROM. In other words, the novel CMS-PROM can reproduce the model with almost negligible errors, while the relative errors of the harmonic response obtained using the classical CMS-PROM are pretty large within some frequency ranges.Figure 7
Comparison of harmonic responses between different CMS ROMs. (a) Frequency range: 0–1500 Hz. (b) Frequency range: 1500–1900 Hz.
(a)
(b)To investigate exact errors and assess the accuracy of these CMS-PROMs, the first 5-order resonance responses are compared in Table2.Table 2
Comparison of resonance responses between different CMS ROMs.
Algorithm
Mode
1
2
3
4
5
Direct (dB)
62.4761
−10.1305
−14.4657
−38.0917
−40.6774
Classical (dB)
62.4761
−10.7511
−15.5202
−37.9199
−40.8626
Relative error (%)
0
6.126
7.290
0.451
0.455
Proposed (dB)
62.4761
−10.1395
−14.4689
−38.0975
−40.6728
Relative error (%)
0
0.0888
0.0221
0.0152
0.0113As shown in Table2, the results associated with the proposed approach show a higher level of accuracy than the other classical approach. The results also reveal that the relative errors of the proposed PROM are below 0.1%, and very good predictions can be obtained using this novel approach.To sum up, all of the above results illustrate the fact that the novel component-based parametric reduced-order modeling method can not only sample the parameter space automatically and efficiently but also lead to more accurate ROMs for further calculation.


## 5. Conclusions

An automatic component-based parametric reduced-order modeling method has been developed in this paper. On the basis of introducing the CB-CMS, we combined it with the congruence transformation method and local interface reduction method to achieve the compatibility between discrete CMS ROMs. The automatic sampling procedure is then described in detail. During this process, we defined the error indicator to locate next optimal sampling points and the exact error to monitor the convergence rate and termination criteria. In addition, the entire analysis procedure was summarized in the schematics. For verification of the presented method, we analyzed the moving coil of electrical-dynamic shaker. Numerical experiments illustrate that this novel PROM can provide perfect predictions of frequencies and harmonic responses. The proposed approach can automatically determine the minimum parameter sample points to satisfy the predefined error and combine the benefits of component-based parametric reduced-order modeling method and those of automatic sampling procedure.



---

*Source: 6407437-2019-07-08.xml*