Abstract
This study presents a technique that uses a model reduction method for the dynamic response analysis of a beam structure to a moving load, which can be modeled either as a moving point force or as a moving body. The nature of the dedicated condensation method tailored to address the moving load case is that the master degrees of freedom are reselected, and the coefficient matrices of the condensed model are recalculated as the load travels from one element to another. Although this process increases computational burden, the overall computational time is still greatly reduced because of the small scale of motion equations. To illustrate and validate the methodology, the technique is initially applied to a simply supported beam subjected to a single-point load moving along the beam. Subsequently, the technique is applied to a practical model for wheel-rail interaction dynamic analysis in railway engineering. Numerical examples show that the condensation model can solve the moving load problem faster than an analytical model or its full finite element model. The proposed model also exhibits high computational accuracy.
1. Introduction
A structure subjected to moving loads is a common situation in mechanical and civil engineering. Some examples include a shaft workpiece to a lathe tool, a bridge to cars and vehicles, and a rail to cranes or rail vehicles. Various models have been proposed for this dynamic analysis. Although they may be applied in different engineering domains, these models share similar characteristics. Two approaches are generally employed to build a model, that is, an analytical method or a finite element (FE) method.
Various analytical solutions to a uniform simply supported beam with a single span exist in the excellent monographs of Frba [1]. Mamandi and Kargarnovin [2] modeled a continuous beam with Timoshenko beam theory which is applicable for beams with a large height to span ratio. For slender beams, Bernoulli-Euler beam theory is promising. Hayashikawa and Watanabe [3] used a continuum method of dynamic analysis for multispan continuous beams. Zheng et al. [4] modified beam vibration functions, which are used as assumed modes to analyze the vibration of multispan nonuniform beams. Lalthlamuana and Talukdar [5] investigated dynamic interaction between a single-span bridge and moving flexible vehicles using a semianalytical approach. Analytical methods can provide accurate solutions to simple moving load problems, but they are sometimes cumbersome for modeling complex practical structures.
Alternatively, the FE method has been used as a versatile numerical approach to various complicated structures. Wu and Dai [6] used a transfer matrix method to evaluate the natural frequencies and mode shapes of a beam and applied a mode superposition method for the dynamic response analysis of the beam under a moving load. Wu et al. [7] developed a general analytical procedure to allow standard FE analyses to deal with the dynamic behavior of 3D structures with bodies moving in two dimensions. Rieker et al. [8] investigated the relationship between model accuracy and the number of elements used to discretize a structure for a moving load analysis.
However, the number of unknowns involved is frequently large in practical FE models. Thus, the computational burden may be great. A condensation method or a substructure approach is employed to reduce the operation scale. Yang and Lin [9] developed a vehicle-bridge coupled model, in which the bridge element and a part of the vehicle body resting on it were regarded as a substructure. The degrees of freedom (DOFs) associated with the vehicle body within the substructure were eliminated through a dynamic condensation method (DCM) based on iteration. Biondi et al. [10] proposed a variant component-mode synthesis method to model the train-track-bridge system. de Salvo et al. [11] tailored a variant of the component mode synthesis method to deal with multispan continuous beams by decomposing the entire structure in primary and secondary spans. Exact eigenfunctions were used as assumed local modes and the beam was reassembled to consider kinematical compatibility between spans.
To illustrate the methodology of employing a model reduction method to solve the moving load problem, three typical model reduction methods, namely, DCM [12], the improved reduced system (IRS) method [13], and the system equivalent reduction-expansion process (SEREP) [14], are introduced in this study. They serve as examples for investigating whether a reduction method can be employed to accelerate the dynamic response analysis of a beam structure under a moving load. The Guyan reduction method [15] performs well in static cases but is unsuitable for this dynamic moving load case. DCM is a representative physical model reduction technique. IRS combines the Guyan method and DCM. SEREP is a typical model reduction method. Other reduction methods are derived from these methods, but they are not discussed in this paper. Given that the load is moving along the beam, the master coordinates of the condensed model should be reselected; these coordinates are accompanied by a recalculation of the coefficient matrices. First, the beam is meshed into a series of equal beam elements. Then, the global system matrices and load vector of the full FE model are obtained. Second, based on the load location, the master DOFs of the condensed model are selected, and the elements in the global system matrices and the load vector are rearranged. Third, the transformation matrix that establishes the relationship between the master DOFs and the original coordinates is obtained. The system matrices and the load vector of the condensed model are calculated. The dynamic response is solved using a direct time integration method. Fourth, when the load moves on to the next element, the initial conditions of the new master coordinates are determined. Finally, by repeating the second to fourth steps, the operation is implemented continuously as the load traverses the entire beam. To expound and validate the methodology, the method is applied initially to a moving point load case and subsequently to a moving body case in practical railway engineering. Accuracy and high computational efficiency are discussed through numerical examples.
2. Moving Load FE Models
2.1. Full FE Model
A simply supported beam with a uniform cross-section is built to illustrate the methodology. The beam is subjected to a constant load . The load location is indicated by that traverses the beam at a constant speed , which starts at the left end and travels to the right end, as shown in Figure 1.
An FE method is employed to discretize the beam, as displayed in Figure 2.
The simply supported beam is discretized into nodes connected by 1D Bernoulli-Euler beam elements (the small circles on the beam denote the nodes of FEs rather than their hinges). The length of a single-beam element is , which is given by . The Arabic numerals beneath the beam denote the serial numbers of the nodes, whereas those above the beam represent the serial numbers of the elements between two adjacent nodes. Each node has two DOFs, that is, translation and rotation , . Therefore, each beam element has four DOFs, and the entire beam has DOFs (the vertical displacements of two beam endpoints are also counted, although they are always equal to zero). The motion of the beam is ruled by the following equation: where of order is the structural dynamic or effective stiffness matrix. and are the structural mass and stiffness matrices, respectively. of order is the displacement response vector, which is given by and of order is the external load vector. Considering that the external load can only be applied in nodes in the FE model, the point moving load must be converted into equivalent nodal forces and moments based on the generalized force equilibrium condition. The process is illustrated in Figure 3.
The equivalent forces or moments , , , and are evaluated by the external load and , which is the distance along the element to the application point of . Therefore, the overall external force vector is where The entries of matrix , excluding , , , and , are zero elements, which are represented by the symbol “⋯.” The variables on the right side of the column vector denote the position of the nonzero entries in the vector. is the serial number of the beam element where the external force is located. The concentrated force moves along the beam. Accordingly, the values of the four nonzero entries change with the periodical variation of . When the force moves to the next beam element, the position of these entries within the column vector shifts correspondingly.
2.2. Condensed FE Models
To make the model computationally efficient for solving the aforementioned problems, the condensed model is formulated using the dedicated variants of the two model reduction techniques tailored to the moving load case.
Given that the external force exerted on the beam is equivalent to two forces and two moments imposed on four DOFs, four coordinates are maintained in the condensed model. These coordinates are the DOFs of the beam element that is posited by the exerting force. The other DOFs are regarded as slaves. When the force moves on to the next beam element, the master DOFs have to be reselected accordingly; that is, the four DOFs for the next beam element are chosen as the master DOFs, which are shown in Figure 4.
The bold nodes of the beam model in Figure 4 are the ones maintained. Accordingly, the bold entries in the column vector are chosen as the four master DOFs. The four DOFs of the next beam element are selected as the master DOFs when the load moves on to the next element.
2.2.1. DCM
Assuming that the force is located on the th beam element, (1) can be partitioned into master DOFs () and slave DOFs to form two equations written in matrix form as follows: where is the frequency that forms the matrices of the condensed model. Accordingly, where superscript refers to the serial number of the beam element where the force is located. is the equivalent force vector and can be written alternatively as follows: where
Given that equivalent loads are only exerted on the master DOFs, no load acts on the slave DOFs. Thus, in (6). The second equation can be written as follows: which can be solved for the displacement at the slave DOFs as follows: This equation can be manipulated to obtain Introducing (13) into (6), premultiplying both sides by , and considering , the motion equation of the condensed model when the load acts on the th beam element is derived as follows: or where , , and are the dynamic stiffness, mass, and stiffness matrices of the condensed model, respectively. These variables are determined as follows:
2.2.2. IRS
The coefficient matrices of the condensed model are derived once more via IRS, which combines the Guyan method [15] and DCM by adding a term to the static transformation to consider inertial forces.
The Guyan reduction scheme is derived. Omitting the inertial item in (6), that is, considering the static form of motion equations, yields
Assuming that the forces on the slave DOFs are zero, the second equation can be solved for the displacement at the slave DOFs as follows: which can be manipulated to obtain Substituting (19) into (6), the eigenvalue problem of the condensed model based on the Guyan method is presented as follows: or where and are the natural frequency of the reduced model and its corresponding eigenvector at the master DOFs, respectively. and are the condensed mass and stiffness matrices obtained using the Guyan method, respectively. These variables are determined as follows: Using the binomial theorem, (12) can be rewritten as follows: where denotes an error of order . Introducing (21) into (23) and ignoring the terms in yields Thus, the relationship becomes available by relating the master DOFs to the full set of DOFs as follows: where This transformation can be regarded as a perturbation on the transformation from the static reduction method. Therefore, using the transformation matrix, the condensed mass and stiffness matrices obtained by IRS are respectively written as follows:
2.2.3. SEREP
When the coefficient matrices of the condensed model are derived using SEREP, the physical coordinates of the original model can be expressed in modal space as follows: or where is the modal matrix with columns that consist of modal vectors. is the modal coordinate vector. Considering only the master coordinates in (29) yields the following: In this case, the number of retained modes is assumed to be larger than the number of master DOFs to investigate the high-frequency vibration of the beam. Thus, the matrix is rank deficient. , which is nonsingular, is employed as an alternative solution to overcome this problem. Equation (30) is solved for the modal coordinate vector as follows: Substituting (31) into (28) yields an expression for the full coordinate vector in terms of the master coordinate vector, as follows: Therefore, the condensed mass and stiffness matrices obtained by SEREP can be respectively written as follows:
The dynamic time response of the condensed model, for example, the displacements, can be obtained by solving motion equations using a direct step-by-step time integration method, such as the central difference method.
Thus, the displacement of the beam at a contact point can be calculated with the four nodal displacements by the Hermitian interpolation function, which is given by where is the shape function determined as follows:
The master DOFs must be reselected as aforementioned when the force moves across the node and on to a new element. The stiffness and mass matrices of the condensed model should also be recalculated with the rearrangement of the coordinates in the original model.
The central difference integration method is a two-step algorithm. Thus, the displacements at moment can only be calculated with known displacements at two previous moments, that is, and , which are given by where , , and denote the displacements of the master DOFs at moments , , and , respectively, when force impacts on the th beam element. , , and are the constant coefficient matrices that are the functions of stiffness, the mass matrices of the condensed model , and the force vector , respectively.
Equation (36) is applicable in most cases, except in two cases in which the force has just crossed the node and moves on to the next element after only one or two time steps. These two cases are theoretically the same. Accordingly, only the first case is discussed in this paper. The situation is as follows. The force is one time step away from the left node of the th element, as illustrated in Figure 5.
(a)
(b)
(c)
On this occasion, and in (36) cannot be calculated because the master DOFs at these two time steps are the four nodal DOFs of the left (or prior) beam element , , . Consequently, cannot be directly calculated.
is used as an example to determine and . The transfer matrix is employed to establish the numerical relationship between and . The displacement vector of the full FE model (with entries in incorrect order) can be achieved by premultiplying by ; that is, . The four entries in vector are the third, fourth, th, and th entries of vector , respectively. Vector can be expressed in matrix form as follows: where is a 4 by matrix that is determined as follows: The serial numbers above the matrix denote the position of the nonzero entries in the matrix. The entries in the matrix, excluding the four unit entries, are zero entries.
Similarly, the transformation at moment is as follows: Substituting (37) and (39) into (36), can be addressed.
The computation process of is formally the same as that of .
Thus, the response of the condensed model can be solved continuously as the force moves from one element to the next. Time integration methods apart from the central difference method, such as Newmark and Wilson- methods, can also be used to solve the problem. The difference is the coordinate transformation of the master DOFs when the force moves from one element to the adjacent element.
3. Wheel-Rail Contact FE Model (Moving Body Model)
The model reduction method used to model the moving load problem can be extended to the case of actual railway engineering. In this case, the concentrated moving force is replaced by a double mass spring damping system that represents the railway vehicle wheel and sprung mass, which is assumed to run in the horizontal direction at a constant speed and is referred to as the “moving body.” The model is depicted in Figure 6. In this figure, , , , and denote the rail mass per meter, rail density, Young’s modulus of the rail, and the moment of inertia of the rail, respectively. refers to the rail bed spacing.
This model can be used to investigate the wheel-rail contact problem, which is a fundamental subject in railway engineering. To provide prominence to the study of contact behavior, a single-wheel model is adopted. The primary sprung mass is reduced to a lumped mass . The wheelset is modeled as a rigid body , which is linked with through a primary suspension .
The rail is modeled as a uniform Bernoulli-Euler beam meshed into a 1D beam element (the rail section between two rail beds is equally meshed into four beam elements, as shown in Figure 6). The rail is supported by discrete viscoelastic rail beds characterized as a series of springs and dampers . The reduction technique used to condense the rail model is the same as that of the moving load case discussed in the previous section.
The interaction between the moving body and the rail is achieved at the interface through wheel-rail contact force compatibility when considering the instantaneous lost contact of the wheel from the rail surface. The contact force can be described by nonlinear Hertzian theory that is commonly used in the wheel-rail interaction problem, which is given by where is the wheel displacement, is the displacement vector of the four nodal DOFs of the beam element on which the force is exerted (the master DOFs), and is the shape function defined previously. Thus, is the displacement of the rail at the contact point. represents rail irregularity. is the Hertzian contact coefficient determined as [16].
The motion equations of the condensed system are as follows: where and are the rigid body displacements of the sprung mass and wheelset, respectively; and is the vector of the master DOFs of the rail model. in (42) is multiplied by 2, which means that two parallel rails exist beneath the wheelset. The derivation of the condensed coefficient matrices , , and is nearly the same as that of the moving load model. However, the coefficient matrices of the two original models differ because of the springs and dashpots beneath the rail. The detailed derivation process is omitted in this paper. can be expressed by (9) when in the expressions of , , , and is replaced by , which is given by
Motion equations for the overall system are obtained by assembling (41) to (43) into matrix form as follows: is the function of and , as shown in (40). This function can be rewritten as follows: Substituting (46) into (45) yields Equation (47) is a set of complex nonlinear second-order ordinary differential equations. Therefore, a mode superposition method cannot solve the equations. By contrast, a direct time-integration method can handle this problem easily. For example, an explicit integration scheme proposed by Zhai [17] can be used to solve this problem with high efficiency. In addition, the master DOFs can be reselected and the coefficient matrices in (47) can be recomputed. The procedure is the same as that in the moving load case.
4. Verifying the Models
4.1. Verifying the Moving Load Model
The validity of the model is verified by comparing the dynamic response of the condensed model obtained using DCM, IRS, and SEREP with those of the full FE and analytical models. The original model is introduced in Section 2.1. The analytical model has a motion that is governed by partial differential equations that employ mode superposition and central difference methods to solve the problem. The parameters used in the simulation are given in Table 1. The time response curves of the four models are compared in Figure 7.
In Figure 7, the dynamic beam deflection response curve predicted by the full FE model nearly coincides with the analytical result. The three condensed models provide results that are close to the analytical result, as shown in Figure 7. This figure demonstrates the validity of the three condensed models in the given situation. Regarding accuracy and efficiency, the models condensed by IRS and SEREP provide high-precision results that are close to that of the analytical model. The computation times of the three methods are approximately the same. In terms of reliability, the requirement of SEREP is more rigid than that of IRS. The two parameters involved in SEREP, that is, the numbers of retained modes and total beam elements, are limited within a small range and must be precisely chosen to ensure high accuracy. Consequently, SEREP is less suitable to handle this problem than IRS, and IRS is demonstrated to be more feasible than DCM and SEREP to solve the moving point load problem.
4.2. Verifying the Wheel-Rail Contact Model
The prediction of the moving body responses can be considered to check the reliability of the introduced theory and the developed computer program. For example, passing a rail joint at a constant speed can be compared with the response of the analytical model, in which the rail motion is ruled by partial differential equations, and the first several models are kept. Formulating the analytical model is theoretically the same as that of the Zhai model [16]. The rail joint irregularity model adopted in this case is presented in [18]. This model is used as an excitation source to check the dynamic response of the models.
The parameters used in the simulation are listed in Table 2. The dynamic response calculated by the three models is presented in Figure 8.
(a)
(b)
As shown in Figure 8, the dynamic shock response calculated by IRS is in good agreement with that of the analytical model. The prediction of impact contact force by IRS is slightly higher than that of the analytical model. The relative error of the maximum contact force based on the analytical solution for IRS is 8.4%. The maximum deflection of the rail at the impact point calculated by IRS is nearly the same as the analytical result but with a slight error. The fluctuation of the dynamic response curve is inconsistent with the analytical result.
For DCM, the amplitude of the simulated maximum high-frequency contact force is low, which underestimates the high-frequency vibration amplitude of wheel-rail interaction. The vibration period slightly lengthens. The relative error of the maximum contact force is . However, the rail displacement response curve is not in good agreement with the analytical solution.
SEREP exhibits nearly the same accuracy as IRS. However, the number of retained modes in SEREP must also be precisely chosen because a large or small number of retained modes can make the result inaccurate. Hence, SEREP is inappropriate for practical applications. In addition, the computation time of SEREP is 39.4 s, which is considerably greater than the computation time of IRS (i.e., 7.6 s). This difference is attributed to the resolution of the eigenvalue problem when the mass body moves on to a new element. Consequently, SEREP is meaningless because a model reduction method intends to reduce computation time.
Whether in a moving load or a moving body case, IRS can provide results that are close to that of the analytical solution with fast speed. By contrast, DCM cannot efficiently predict the dynamic response of a moving body case, and its accuracy for the moving load case is lower than that of IRS. SEREP is also computationally inefficient and inconvenient for engineering applications. Therefore, the condensed wheel-rail contact model obtained using IRS can be employed as a valid and highly efficient tool to simulate the dynamic response of the wheel-rail interaction problem.
5. Discussion on Computational Efficiency
5.1. Computational Efficiency of the Moving Load Model
The main objective of a model reduction technique in implementing this problem is to reduce computation time. Based on the discussion in the previous section, the model condensed by IRS performs the best among the three condensation methods. Thus, the extent to which the condensed model obtained using IRS can reduce computational burden must be investigated. The comparison of computation times between the original and condensed models with different (i) beam lengths or (ii) numbers of beam elements adopted is shown in Table 3. and denote the computation times of the original and condensed models, respectively.
Table 3 lists the computation times of the simulations by the original and condensed moving load models at different parameters and by the percentage decrease of the computation time between them. The spaces without data are situations with low accuracy. This issue is discussed in Section 6. The condensed model can reduce approximately three quarters of the computation time of the original model when the number of elements is 30 and the beam is 5 m long. A large number of elements indicate the great extent to which the condensed model can reduce time. However, the results vary slightly as beam length changes.
5.2. Computational Efficiency of the Wheel-Rail Contact Model
The computational accuracy of the wheel-rail contact model or the moving body model is shown in Table 4.
Solving the original wheel-rail contact model requires considerable time because of its relatively large DOFs. Data are not shown in Table 4. This table presents two cases with different rail lengths. The two cases provide nearly the same simulated dynamic results. Computation time can be reduced by 67.8% and 47.7% compared with the analytical model. High calculation accuracy is maintained, as discussed in Section 4.2. The proposed model can provide good simulation results with minimal DOFs involved. Therefore, the condensed FE model obtained using IRS can be employed as a highly efficient tool to simulate wheel-rail interaction behavior.
6. Discussion on the Moving Point Load Model
Based on the conclusions drawn in the preceding section, the highly accurate condensed model by IRS is discussed as follows.
The accuracy of the model may be influenced by several factors, for example, (i) the number of elements used to mesh the beam, (ii) beam length, and (iii) the velocity of the moving load.
The response predicted by the condensed FE model (IRS) is compared with that of the analytical solution by discretely varying the value of the factors to find their optimal value. The influence of such factors on accuracy is also investigated and illustrated in Figure 9. A relative error is cited as the parameter that describes model accuracy. This parameter is defined as the ratio of the absolute value of the deflection response difference between the solutions of the two models to the deflection calculated by the analytical model, and is given by where and are the dynamic beam deflection values calculated by the analytical model and the condensed model (IRS), respectively.
(a)
(b)
(c)
The normalized traveling time of the moving load is shown at the horizontal axis in Figure 9. The scales 0 and 1 represent the beginning and end times, respectively. In Figure 9, the large errors that occur at the beginning and end times are the result of the minimal value of beam deflection calculated by the analytical model because this value is close to zero and functions as the denominator of the relative error expression, which can make it computationally divergent.
Figure 9(a) shows the relative error with 4, 7, 10, and 13 beam elements for the same frequency and velocity. The accuracy of the solution is highest when the beam is meshed into 7 or 10 beam elements. The accuracy of the model with the other elements is lower because of two reasons. First, the low mesh density in the moving load-type problem cannot properly reflect the deformation between the nodes and the continuous nature of the beam. Second, the high mesh density can overcome the aforementioned defects, but the ratio of the number of master DOFs to that of the original model is small. Consequently, the representation of the condensed model for the original model is affected. The solution generally becomes accurate when the ratio is high, except in situations in which the number of DOFs in the original model is too small.
Figure 9(b) shows the relative errors with beam lengths of 2, 5, 10, and 15 m for the same number of elements and velocity. High accuracy is obtained when the beam is 5 m or 10 m long. When beam length increases, accuracy decreases.
Figure 9(c) shows the relative errors with moving load velocities of 100, 200, 300, and 400 km/h for the same number of elements and frequency to construct the model. The highest accuracy is generated when the load moves at a velocity of 200 km/h to 300 km/h.
The number of elements adopted is the most significant influence factor. The optimal number of elements adopted depends on the beam length and velocity of the moving load. The numerical example shows that a long beam length permits the minimal number of total elements to be adopted.
Therefore, the influence of the combined effects of these three factors on the accuracy of the condensed model is further investigated to determine the optimal number of elements required in practical applications, as shown in Figure 10. The three coordinate axes are beam length, moving load velocity, and element number. The size of the dots represents the average relative error, , which is defined as follows: where ) is the relative error at the th time step. is the total number of time steps that is determined as follows: is the coefficient that is defined as follows: and is the number of with a value of 1. The at both sides of the beam is not considered because the numerical error caused by the small value of beam deflection is too large, such that the calculated relative error becomes meaningless, as discussed earlier.
In Figure 10, the size of the dots is proportional to the average relative error, , which is calculated using (49). Small dots denote an average relative error or high calculation accuracy. The graph shows that when the element number is approximately 8 or 10, computational accuracy is generally the highest, which provides a similar result to the aforementioned statement.
The graph is practically significant because it can help determine the optimal value of the element number when a specific beam and the moving load velocity are provided. For example, if the beam is 10 m long and the load is moving at a velocity of 300 km/h, then the highest calculation accuracy is achieved when the beam is meshed into 10 beam elements that correspond to the smallest dot along the vertical auxiliary line.
7. Conclusion
To improve computational efficiency, three dedicated model reduction techniques are introduced to condense the FE model for the dynamic response analysis of a beam-type structure that is subjected to a moving load. The same methodology is used to reduce the size of the moving body problem as a practical model to study the wheel-rail contact behavior in railway engineering.
Applying the model reduction method on a moving load case is different from applying it on an immobile load case in several aspects. First, the master DOFs of the condensed model must be reselected when the load travels from one element to another, while the number of it remains unchanged. Correspondingly, the coefficient matrices should be recalculated based on the rearrangement of the coordinates in the response vector of the original model. Second, the value of the nonzero entries in the force vector varies periodically as the load travels along the beam, and their position changes as the load moves on to the next beam element. Third, to guarantee the continuity of operations, the initial conditions of the condensed model must be redefined as the load moves on to a new element because of the reselection of the master DOFs.
To demonstrate the validity of the proposed model, numerical examples are provided to compare the simulated solution of the condensed FE models obtained by three model reduction techniques, namely, DCM, IRS, and SEREP, with the analytical solution. The results indicate that the condensed model obtained using IRS generally performs the best among the three condensation models and can provide reliable results with fast computational speed under given conditions.
The computational efficiency of the condensed model obtained using IRS is investigated. In the moving load case, the condensed model can reduce computation time by approximately three quarters of that of the original model on the premise that it ensures a sufficiently accurate solution. In the moving body case, the computational burden of the full FE method is too high. By contrast, the condensed model can reduce computation time by more than half of that of the analytical model with excellent accuracy.
A deficiency of the proposed moving load model is that low accuracy is exhibited when the ratio of the number of master DOFs to that of the original model is low and the beam lengthens. However, the quantitative relation among these variables is difficult to achieve. This low-accuracy problem may be solved by increasing the number of master DOFs or adding several DOFs to the immobile nodes in the master coordinates. This problem will be discussed in a future work.
The model reduction technique is only applied to the beam element in this study. However, this technique can theoretically be generalized and adapted to any element type. Accordingly, further work must be conducted before this methodology can be considered as suitable for broad applications.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
This work was supported by the National Natural Science Foundation of China under Grant no. U1234208 and the Open Project Foundation of the State Key Laboratory of Traction Power under Grant no. TPL1310. The authors would like to extend their gratitude to all the sponsors for their generous financial aid.