Abstract
The mechanical response of multilayer sheet structures, such as leaf springs or car bodies, is largely determined by the nonlinear contact and friction forces between the sheets involved. Conventional computational approaches based on classical reduction techniques or the direct finite element approach have an inefficient balance between computational time and accuracy. In the present contribution, the method of trial vector derivatives is applied and extended in order to obtain a-priori trial vectors for the model reduction which are suitable for determining the nonlinearities in the joints of the reduced system. Findings show that the result quality in terms of displacements and contact forces is comparable to the direct finite element method but the computational effort is extremely low due to the model order reduction. Two numerical studies are presented to underline the method’s accuracy and efficiency. In conclusion, this approach is discussed with respect to the existing body of literature.
1. Introduction
Some examples of permanent joints between two solid bodies are screwed or bolted joints, crimp connections, and spot-welded seams. The construction of such joints normally permits large relative displacements between the contacting bodies. This leads to a time invariant potential contact area , which is identical to the spatial distribution of the joint. However, due to small structural displacements there are areas inside the joint where the surfaces are in contact and other areas where the surfaces are not in contact. Consider The argument represents this structural displacement and indicates that the contacting and gapping areas are displacement dependent. In dynamic cases, this implies an indirect dependency on time. The state dependent contact forces, which avoid penetration, and the friction forces act between the two contacting surfaces. Such friction forces are a function of the local contact forces and the relative in-plane displacements of the contacting surfaces. Due to the joint’s construction, such relative displacements are small and are, therefore, often referred to as “microslip” or “small sliding” displacements. These friction forces lead to energy dissipation. In summary, it can be said that such joints lead to nonlinear stiffness and damping effects. The latter observations are well documented in the literature (see [1–7] for more detailed information on the subject of experiments and simulations). Depending on some parameters, such as the spatial extension of the joint, the aforementioned nonlinearities can have an effect—both minor and major—on the structural response (see [8] for an investigation on this issue with two generic structures). In some cases it is of significant importance for the response quality to regard the nonlinear characteristic of the joint in terms of contact and friction. Examples of this are structures with large joints, like car bodies or leaf springs. This is demonstrated in the paper as well. However, even if the global effects of the nonlinearities in a joint are negligible, it can happen that the local ones are not. An example of this is the lifetime prediction of a weld spot. There is a significant difference on the local stress situation around a weld spot, whether or not the reinforcing effect due to the contact between the surrounding sheets is taken into prior consideration (see [9]).
One possible approach to arrive at a mathematical model which considers the contact and frictional phenomena is the use of the finite element method (FEM). In case of a penalty formulation for both the contact [11] and friction [12] phenomenon, the FEM leads to a coupled and nonlinear differential equation where the () vector x contains the system’s degrees of freedom (DOF). Note, no material nor any viscous damping is considered in (2). The constant () matrices and represent the structural mass and stiffness. The () vector holds the external and imposed forces while the () vector contains the nonlinear forces due to contact and friction inside a joint. A direct time integration of (2) is in principle possible, but it is generally considered an inefficient operation due to the system’s dimension , which can be quite large in the case of relevant finite element (FE) models found in the industry. Such FE models may have 107 DOF and this number is still growing. Even though the resulting quality may be good, the computational power required can be impractical.
In order to shrink the number of DOF in (2), model order reduction via projection is possible. The displacement vector is approximated by which is limited to a scaled superposition of time invariant trial vectors. Formally, this can be given as where the () vector holds the time invariant th trial vector and the corresponding time varying scaling factor. All trial vectors can be collected as columns in the () matrix and the scaling factors are collected in the () column vector . The application of (3) on the equation of motion (2) leads to instead of DOF, which will be shown in the next chapter in greater detail. Good reduction approaches provide a number of which is much smaller than with a small error where x is the system response of (2) and is the response of the reduced system. A proper vector norm is symbolized by . For linear systems, where the stiffness matrix is constant, a number of suggestions for a proper trial vector base can be found in the literature (see [13–15]). Comparisons can be found, among others, in [16–18]. Properly applied, the latter methods provide a proper reduction base for linear systems. However, in nonlinear systems it is more challenging to find a proper subspace . The literature offers some choices for model reduction of nonlinear systems, namely, the conservation of all DOF which are involved in nonlinearities, proper orthogonal decomposition (POD), and trial vector derivatives.
Nonlinearities which are restricted to certain (time invariant) areas of a structure are often referred to as “local nonlinearities.” By this definition, joints are such local nonlinearities. Some publications in the literature suggest preserving all the nodal FE DOF which are involved in nonlinearities (see [19–24]). The consequence for jointed structures would be that each FE DOF in a joint would lead to a corresponding DOF in the reduced system. With joints, this approach would lead to several hundreds or even thousands of DOF, which then need to be considered in the time integration of the reduced system. The computational effort would, again, be very high. Some other publications suggest the use of POD for the construction of a proper subspace (see [25] for an overview on POD). Note that in Section 2 additional literature is given concerning POD and the concept is briefly discussed. The basic idea of the so-called POD-snapshot method is to determine the most important subspace of a given space in terms of a Euclidian distance. The given space is spanned by all result vectors to which are obtained by time integration of the full system for the time steps to . An obvious disadvantage of this approach is a necessary simulation of the full system in order to obtain the reduction base. Another option for the construction of a proper subspace is the use of trial vector derivatives. In the literature, the name modal derivative is more common because the term “mode” is often used for any kind of trial vectors. However, in this paper the name “mode” is only used for trial vectors which stem from an eigenvalue problem and the term “trial vector” is used for any kind of displacement shape which is used to form the model reduction base. The basic idea is that for nonlinear systems such as (2) trial vectors are somehow a function of the state. In a first step, a subspace is constructed based on the linear portion of (2). In a second step, the subspace is enriched by an approximation of the difference to those trial vectors, which would have been obtained when the full nonlinear (2) would have been considered. The publication of Slaats et al. [26] gives a very useful introduction to this topic. For the sake of clarity, some parts of the theory are reviewed in Section 2 below. Tiso et al. [27] published a strategy for regarding the nonlinear bending stretching coupling of shells based on trial vector derivatives. Still, an open issue is the a-priori selection of the important trial vector derivatives. This is necessary because trial vector derivatives lead to many potential candidates for the subspace enrichment containing a lot of redundant information. In [27] some general comments are made on this issue and [28] contains a strategy for an “optimal selection.” However, in the present paper, a new strategy will be presented in order to determine the final subspace extension based on trial vector derivatives.
For the particular case, where the nonlinearity is caused by joints, the author has already introduced another approach for the computation of a proper subspace (see [10]). In the latter work, a given subspace based on a linearized structure is enriched by so called Joint Interface Modes. In principle, a meaningful subspace of all DOF inside a joint is determined. This subspace is used as a problem oriented extension of classical reduction bases. Gaul and Becker [29] confirmed the superior convergence rate of the latter Joint Interface Modes (JIM) in comparison to other methods known from the literature. An implementation of this approach is the commercially available software product MAMBA [30]. Segalman [31] suggested at the same time an approach for localized nonlinearities, which is a nonjoint focused generalization of [10]. A comparison of the proposed method with the one of [10] is given in Section 5, which is devoted to a discussion of the proposed method with respect to the already published literature.
The present work is devoted to the application of trial vector derivatives to jointed structures. The proposed strategy leads to an a-priori subspace, which delivers a level of accuracy comparable to the FEM without losing the advantages of model order reduction. Other advantages are as follows.(i)Quick convergence with respect to the number of trial vectors.(ii)Automatic a-priori selection of the necessary trial vectors without a simulation of the full system.(iii)Simple and efficient computation of the trial vectors by utilizing commercially available FE software.(iv)No separate trial vectors in normal and in tangential direction of the contact surfaces are needed.
To the best knowledge of the authors, the following points are new contributions to the literature.(i)Application of trial vector derivatives to jointed structures.(ii)A model order reduction strategy for the efficient simulation of multilayered sheet structures.(iii)A clear a-priori strategy for the automatic computation and selection of the necessary trial vectors using an energy criterion and POD. This is not restricted to joints and, therefore, is a general contribution to the whole issue of modal or trial vector derivatives.
The present work is organized as follows. In the next section the theory of the trial vector derivatives is shortly reviewed and an energy based criterion is introduced in order to determine the important directions—out of all available trial vector derivatives. In the subsequent section, the theory is applied to jointed structures followed by some numerical examples in order to underline the methods’ accuracy and efficiency. Finally, the approach introduced herein is discussed with respect to the already published literature.
2. Trial Vector Derivatives as a General Method for the Generation of A-Priori Reduction Base for Nonlinear Systems
2.1. Nonlinear Mechanical Systems
In the present work we restrict ourselves to structures which can be characterized by an equation of motion in the form of (2). In a next step, the vector x is decomposed into a time invariant portion and a time varying portion , such that The vector can be interpreted as an operating point and the vector as the vibrations around this point. In the very common case that a structure vibrates around its undeformed state, the vector is simply the zero vector. Examples where is not trivial are prestressed screwed joints or leaf springs. As a consequence of the latter decomposition the stiffness matrix simply holds the stiffness at . Rearranging (2) yields which can be interpreted as a linear system which is loaded by external and fictive forces due to the nonlinearities. This is a well-known approach in the literature and the fictive forces are sometimes called “pseudoforces”; see exemplarily [32, 33].
2.2. Model Order Reduction via Projection
As mentioned in (3) the possible displacements are restricted to a weighted superposition of trial vectors. Inserting equation (3) into (6) and premultiplying with delivers with the () reduced mass matrix and the () reduced stiffness matrix . The used subspace needs to fulfill two criteria which influence each other. The subspace should provide an approximated displacement so that accurate nonlinearities can be computed and, as a consequence, the resulting error (4) should be small. For linear systems, the so-called component mode synthesis (CMS) has in recent years become a reliable standard tool (see [13–15, 34]). The subspace used is a combination of two groups of trial vectors where the () matrix contains vibration mode shapes and the () matrix contains trial vectors which are displacement shapes due to static loads. These static displacement shapes are typically related to s nonzero entries in the external force vector . According to the pseudoforce concept motivated by the structure of (6), one could suggest preserving all DOF which are involved in the vector , something which is borne out of the literature (see [19–24]). This approach is promising in case of local nonlinearities with a small number of DOF involved. With distributed nonlinearities, such as joints, the latter approach would lead to a very large number of static modes and the advantage of model reduction is lost. Another “brute force” solution would be to ignore the nonlinearities introduced by . In that case the resulting quality cannot be guaranteed without a convergence analysis with respect to the number of vibration mode shapes collected in . It is shown below that this approach leads to a poor convergence, which means that the number of is high.
2.3. Trial Vector Derivatives
The idea is to extend the trial vector base (8) in the form of where the () matrix contains the so-called trial vector derivatives (TVD). The TVD approximate the difference between and computed with a nonlinear stiffness matrix with respect to and computed with a linearized stiffness matrix (8). Therefore, the computation of can be performed using a finite difference scheme as Slaats et al. proposed in [26].
In general, a trial vector is a function of the stiffness and mass matrix used. Due to the state dependency of the stiffness matrix, a certain trial vector is state dependent as well. Because of the projection (3) the state is a function of the trial vectors weighting factors to . The state dependent trial vector () can be expressed by a Taylor series expansion as The () vector holds the th trial vector computed at the point where the structure is linearized; see (5). This vector is already part of and . Note that the expression holds an by () matrix. One important consideration is to interpret the first order trial vector derivatives as independent trial vectors themselves and add them to the model reduction base. Consider The by matrix contains the first order terms of (10) applied to all trial vectors of and . That means, the trial vector base is increased by its squared column dimension and the number of generalized coordinates collected in the vector increases by the same size. This is, in general, too much for an efficient model reduction scheme. Practical examples show that there is a lot of redundant information in ; see [28]. The space spanned by can be approximated by a reduced space which contains most of the important information. This issue is addressed below, where a POD based transformation (reduction) of into is proposed.
2.4. Trial Vector Derivatives for Vibration Modes (See [26, 35])
When vibration modes are used as trial vectors, the TVD can be computed based on the eigenvalue problem which is given for mode number . A derivation of both sides with respect to the scaling factor of mode number gives The application of the product rule and the assumption that the mass matrix is constant lead to where the () vector holds the TVD of mode with respect to mode and is hereon denoted as . When (13) is premultiplied with , the sensitivity of the eigenvalue with respect to can be computed (see [35]). Note that the symmetry of and has to be considered in this step along with the mass and stiffness orthogonality of the modes. Consider
Substitution of the according term in (13) with (14) delivers the TVD:
It has been observed in the literature [26] that the terms related to inertia can be neglected, also a subject of further investigations in this paper. Doing this leads to a simplified equation in the form of
2.5. Trial Vector Derivatives for Static Displacement Shapes (See [26, 35])
A derivation of both sides of the characteristic equation for statics with respect to the scaling factor of trial vector number gives Considering the product rule and the fact that is state independent delivers which is formally identical with (16).
Equations (16) and (18) can be interpreted in such a way that a TVD is the static response due to the fictive force
This may yield an efficient algorithm for the computation of the TVD when commercially available FE software is used (see Algorithm 1). The advantage of Algorithm 1 is that it uses the strength of such commercially available FE software, which is the computation of static displacements based on a defined number of forces. If several load cases are used for one computational run, the stiffness matrix needs only to be decomposed just once, meaning that a huge number of such load cases can be treated efficiently. However, the challenge may be the computation of . It is discussed in Section 3 how can be computed when nonlinearities due to contact and friction inside a joint are considered.
| ||||||||||||||||||||||
2.6. Comment on a Singular Stiffness Matrix
As outlined in the introduction, the literature offers different approaches for a proper reduction base . Some of them lead to a regular stiffness matrix and some of them do not. The fixed boundary CMS (Craig/Bampton) is an example of a method dealing with an invertible stiffness matrix while the free boundary CMS leads to a singular one (see [14]). Note that a singular stiffness matrix (which corresponds to a floating structure) is not a drawback for the computation of TVD. For such structures the FE software needs to perform a so-called “inertia relief” computation when line 9 of Algorithm 1 is executed.
2.7. Selection of the Final Reduction Base Using Proper Orthogonal Decomposition
Using the TVD provides additional trial vectors for the model reduction. This is by far too much for efficient model reduction and it is known from the literature that the space spanned by contains a lot of redundant information; see [27, 28]. This can lead to numerical problems and is the prime motivation behind the search for a lower dimensional reduction space which contains most of the information of . For some computations in [27] just the TVD have been regarded. In other words, just the change of a trial vector due to its own variation is considered, while the dependency on the other modes is neglected. It can be seen therein that this strategy sometimes leads to good results, but often it does not. In [28] a strategy is presented, which is based on the idea that when the according modes and are important, the TVD is important as well. In the present paper a general strategy for the determination of based on is given. This strategy is based on proper orthogonal decomposition (POD). Literature on POD can be found, among others, in [25, 37]. A detailed derivation of POD can be found in [38]. For the sake of readability, the main points of POD are briefly reviewed next.
Let us assume independent () vectors to which are gathered in the () matrix . POD of rank delivers orthonormal () vectors to , which approximate the space spanned by optimal in a Euclidean sense. The vectors to are named proper orthogonal modes (POM) and can be gathered in the () matrix . The POMs maximize a function in the form of so that (see [38] for a proof). The POM to can be found by using the first eigenvectors of and a subsequent transformation A slightly different POD method with a weighted inner product maximizes a function in the form of so that condition (21) is fulfilled as well. The () matrix has to be symmetric and positive semidefinite. The POM for the weighted POD can be found by using the first eigenvectors of and a subsequent transformation along equation (23). The second POD variant can be transformed into the first one if the () identity matrix is substituted for the matrix . From a mechanical point of view, the characteristics of the matrix are the same as those of the stiffness matrix of a linear FE model. In that context, the POM based on a weighted inner product can be interpreted as those trial vectors which optimally approximate the deformation energy caused by the space . The POM without a weighted inner product approximates the displacement space only. The eigenvalue is of special interest because it is related to the importance of the vector for the minimization of (20) or (24). Note that large eigenvalues indicate important POM. Therefore, a ratio is defined which can be used as a criterion to determine the number of . Consider Note that and . In case of POD with weighted inner product, provides a value of how much of the available deformation energy is already covered by to . Both POD approaches are investigated and discussed in the numerical examples below. In order to apply POD to the problem of reducing the space of available TVD, the following substitutions need to be performed. Instead of the matrix the () matrix has to be used. It is important to mention that all column vectors of are normalized with respect to the stiffness matrix . This is important because the TVD represent just “directions” and the vector’s length is not of interest. Without this normalization, a vector’s magnitude would influence the selection of the “important” POM. The space spanned by corresponds with the one spanned by the () matrix . In case of inner weighting the stiffness matrix is used instead of .
The finally obtained reduction space is .
3. Application of Trial Vector Derivatives to Jointed Structures
The application-relevant part of the general strategy outlined in the chapter above is the computation of the () matrix . The qualitative content of the latter matrix is the change in stiffness due to a variation induced by trial vector number . In cases where nonlinearities induced by joints are used alongside a penalty contact and a penalty friction model (see [11, 12]) this matrix can be obtained as follows. Note that for the sake of simplicity a coincident mesh of the two contact surfaces is assumed and, therefore, it is sufficient to consider a node-to-node contact. Assuming that w FE node-pairs are involved in the joint contact, gap functions can be defined. The gap function represents the normal distance of the node-pair . Following a simple penalty approach, a constant contact stiffness and a constant tangential stiffness are acting between those node-pairs with a gap function smaller than 0 (penetration). If (gaping) no stiffness acts between the node-pairs. We assume that the linear modes and are computed based on the undeformed reference configuration, which means that no contact stiffness is acting between the joint node-pairs. In such a case, the matrix contains the stiffnesses and which act only on those DOF corresponding to node-pairs with a gap function less than zero. The gap functions are evaluated based on the displacements due to the trial vector .
There is another particularity related to nonlinearities caused by the contact and friction in joints. Up to now, we have assumed that all derived quantities are differentiable. This is not the case in terms of the simple contact law under consideration, because Note that the superscript “” denotes an effective stiffness matrix which consists of the linear portion plus the additional stiffness which is added due to contacting joint areas. At that point it would be possible to construct a differentiable contact law using polynomials of higher order. This approach is not constructive because it does not account for the “strong asymmetry” of the nonlinearity under consideration. “Strong asymmetry” means that In other words, the change in the stiffness matrix is completely different in case of a positive or in case of a negative scaling of the trial vector . The mechanical explanation for a particular node-pair is that when a positive scaling leads to a negative gap function and a high stiffness, the negative scaling leads to a positive gap function and no stiffness, and vice versa. Therefore, it is advisable to compute two modal derivatives for a trial vector with respect to trial vector number , namely, and . The () vector denotes the trial vector derivative based on a negative scaling of trial vector . In case of joint induced nonlinearities the matrix is of size (). See Algorithm 2 for the final flow. Note that this flow is valid for structures which are linearized around their undeformed reference configuration. The particular choice of is problem dependent and is discussed in Section 4.
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
4. Numerical Examples
In this section, two practical examples are presented. The first one deals with friction bar containing one joint and the second one is a multilayer sheet structure consisting of three sheets overlapping each other. For both examples, a static and a dynamic computation is presented. The focus is on the convergence of the result with respect to the number of trial vectors considered in .
For both examples, the matrices and were computed according to the Craig/Bampton method [39]. For this particular method, the content of the matrix is named “Fixed Boundary Normal Modes” and the content of the matrix is named “Constrained Modes.” In both examples and are computed with respect to the undeformed reference condition and using zero stiffness between the contacting surfaces.
In both examples the stiffness and mass matrix of the FE model were computed by MSC.NASTRAN [40]. The matrices were imported into Scilab [41] and all subsequent computations have been done directly in Scilab or in its toolbox Xcos [42].
Due to the small sliding assumption and the coincident meshing a node-to-node contact and friction model is applied. The contact forces were computed following a simple penalty model with a linear stiffness of 104 N/mm. The friction forces, if necessary, were computed as given in [43]. The latter model requires two parameters, namely, the friction coefficient and an in-plane stiffness when a node-pair is in contact. For the friction coefficient a value of 0.25 is used and the in-plane stiffness is set to 104 N/mm. The contact and friction forces are computed for all involved nodes pairs and the force vector () is constructed according to the physical DOF. This physical force vector is then projected into the reduction space by a premultiplication with the matrix .
For the time integration the Scilab [41] toolbox Xcos [42] has been used. In Figure 1 a screenshot of the Scilab-Xcos model can be seen. The parameters necessary for the time integration are listed in Table 1.
4.1. Friction Bar
An FE model of the structure under consideration is depicted in Figure 2. Two metallic sheets with a dimension of 125 mm × 25 mm × 0.5 mm are connected via two beams. The extension of the joint is 85 mm × 25 mm. The length of the overlapping area (the joint) is indicated in Figure 2 with a thick black dashed line.
The two sheets and the two beams, which can be interpreted as two weld spots, are modeled out of steel. On the left hand side, the structure is fixed mounted. On the free end, all nodal DOF are attached to a rigid bar. The free node of the rigid bar is in the middle of the sheet. This is indicated in Figure 1 by a solid dark grey line and a dot. The load for the static and dynamic investigations consists of a force in the direction () and a torque around the x axis (). The entire model has 906 DOF and 648 DOF involved in the joint contact.
4.1.1. Trial Vector Derivatives
As mentioned before, the reduction base for the linear system is computed in accordance with the Craig/Bampton method [39]. The six DOF of the free end are defined as so-called “boundary DOF.” This leads to six static deflection shapes for the matrix . The two most relevant “Constraint Modes” (deflection along and rotation around ) are visualized in Figure 3.
For the matrix the first ten “Fixed Boundary Normal Modes” are considered. The eigenfrequency of the 10th mode is 1337 Hz. Figure 4 contains a visualization of the first three modes (=columns) of . It can be seen that the first two columns represent bending modes and the third column represents a torsion mode around the x axis. Note, again that the “Constraint Modes” as well as the “Fixed Boundary Normal Modes” have been computed based on a linear structure where no contact and friction are considered.
In the next step, the trial vector derivatives according to the theory outlined in Sections 2 and 3 are computed. In case of 16 modes the procedure ends up with 512 TVD. Figure 5 contains 6 arbitrarily selected TVD. As mentioned in the literature TVD typically contain a degree of redundant information. This can already be seen in Figure 5 because and or and look quite similar.
In order to get a suitable small number of trial vectors for the joint, the POD method, as outlined in Section 2, is applied to the TVD. In this example, both methods, POD with and without weighting, are applied. Figure 6 contains visualizations of the first four columns of the matrix . The visualized vectors are based on weighted POD.
The question of which one of the two POD approaches is preferable will be answered later, when the results are evaluated.
4.1.2. Static Response Computation
The static response is computed for N and Nmm. For this computation, no friction is considered. For the convergence analysis, a () vector is constructed. The vector contains w entries, namely, to , and () maintains the contact force at FE node-pair . As convergence criteria the Euclidean norm of the vector is evaluated. Figure 7 contains a convergence analysis with respect to the number of considered joint vectors in the matrix (=g). Note that the columns of this matrix are donated as “joint mode” later on. For the reference solution the contact forces are computed based on a fully nonlinear FE analysis using all nodal DOF. The dashed and the dotted curves represent the Euclidian norm of the contact force vector following the proposed approach. In order to underline the necessity of such modes, an additional computation is performed without joint modes. Instead of additional joint modes the number of vibration modes is increased. The results of this approach are documented in the dashed and dotted lines.
It can be seen in Figure 7 that 20 additional joint modes lead to a good representation of the contact force. The convergence in cases where vibration modes were only used is not that good. Note that there is no significant difference when POD with or without weighting is used. From this point of view, both computational approaches qualitatively lead to the same result. Figure 8 contains the function for both approaches. The scalar value g represents the number of trial vectors finally considered in .
A comparison of Figures 7 and 8 indicates that in case of weighted POD, the function corresponds better with the convergence of the structure. This is not a proof, but an observation which can be made for other structures as well. The latter observation can be used as an a-priori estimation for the number of necessary modes. For the friction bar the solution can be considered as converged when 25 additional joint modes are taken into consideration. This corresponds with and that means that 99.94% of the available strain energy of all TVD is considered in the final reduction base. Based on this observation we propose a number of such that The actual limit for the function is problem dependent and probably needs some practical experience using the proposed method.
4.1.3. Dynamic Response Computation without Friction
For the dynamic response analysis the former loads are applied as step functions to the structure. The step function is implemented as a half wave cosine with a frequency of 10 Hz.
Figure 9 contains the evolution of the norm of the contact force vector with respect to time. It can be seen that a matrix with 25 columns delivers excellent results which corresponds with the convergence analysis of the static computation.
As mentioned in Section 2 the trial vector derivatives for vibration modes can be computed with or without neglecting the inertia effect. Figure 10 contains the evolution of the Euclidian norm of the contact force vector for joint modes according to TVD following equation (15) and those according to (16) where the inertia effects are neglected. It can be seen that the consideration of the inertia terms in (15) does not lead to a better convergence. This is in accordance with the literature; see [26]. Therefore, it is advisable to compute the TVD based on (16) which is much simpler and holds some advantages with respect to commercial FE codes (see Section 2).
4.1.4. Dynamic Response Computation with Friction
In case of a metallic structure the dry friction inside a joint is typically the dominating source of damping. The reader is referred to the introduction in Section 1 for more explanations and literature quotations on that issue. In the cited literature it can be read that a Coulomb, like dry friction model which is applied between all FE node pairs, catches most of the relevant effects.
Figure 11 contains a convergence analysis and it can be seen that the consideration of no additional joint modes is not an option. The explanation can be found in Figure 12 where the Euclidian norm of the () friction force vector is evaluated with respect to time. The vector contains friction forces in the and directions of the FE node-pairs being involved in the joint contact.
Figure 12 indicates that an accurate joint flexibility is essential for the result quality when friction is regarded. When no additional joint vectors are regarded, the damping due to friction is extremely overestimated. The difference between 10 and 25 additional joint modes is already small and a converging behavior very similar to the example without friction can be observed. In cases of friction, it was not possible to compute a reference solution using all FE nodal DOF.
Note that up to now the TVD have been computed without the consideration of the tangential stiffness when the contact is closed. In other words, the matrix only contains the change in stiffness due to the contact in case of a closed gap. In a next step the trial vector derivatives are computed regarding a normal and an in-plane stiffness in case of a closed contact. The result is not plotted in a diagram because there is no notable difference to the curves in Figure 12. In other words, the in-plane displacement in the joint is dominated by the flexibility in the joint normal direction while the contributions due to the in-plane flexibility can be neglected.
4.1.5. Computational Efficiency
All simulations are performed with the same XCos [42] block diagram (see Figure 1). The simulations just differ in the particular use of the stiffness and mass matrix and the number of DOF. While the reduced models require a reduced stiffness and mass matrix the full FE model needs the matrices obtained by the FE software. It can be reported that the time integration with 0, 10, and 25 additional joint modes needs 0.07%, 0.12%, and 0.25% CPU time in comparison to the full FE model. In other words it can be said that the use of 25 additional joint modes delivers an accuracy which is comparable to that of a full nonlinear FE simulation without losing the advantages of model reduction.
4.2. Multilayer Sheet Structure
A 3D view and a front view of the FE model of the structure under consideration are shown in Figure 13. Three metallic sheets of the dimensions 135 mm × 25 mm × 0.5 mm are connected via three beams. Due to its construction the structure contains two joints with an extension of 135 mm × 25 mm. The three sheets and the three beams, which can be interpreted as two weld spots, are made out of steel. On the left hand side, the structure is fixed mounted. On the free end, all nodal DOF are attached to a rigid body. The free node of the rigid body is in the middle of the sheets. This is indicated in Figure 13 by a white “spider” in the enlarged picture of the free end. The load for the static and dynamic investigations consists of a force in the direction () and a torque around the x axis (). The entire model has 834 DOF and all of them are involved in the joint contact.
In order to demonstrate a special advantage of the proposed method in contrast to an already existing method, the structure is artificially subdivided into two regions which are named joint area 1 and joint area 2 (see Figure 13). This subdivision is obtained by displacement constraints between these two areas. As a consequence, joint area 2 is in fact not present for the problem under consideration. This is done to demonstrate an advantage of the proposed method in contrast to the one published in [10]. The latter mentioned advantage is discussed in detail in Section 5.
As in the previous example, the reduction base for the linear system is computed according the Craig/Bampton method [39]. The six DOF of the free end are defined as so-called “boundary DOF.” Consequently, the matrix which holds the static deflections shapes contains six trial vectors. The number of vibration modes which are stored in the matrix is 10.
4.2.1. Trial Vector Derivatives
The trial vector derivatives and the final used joint modes are computed along the theory given in Sections 2 and 3. Note that POD with weighting has been used to get the final joint modes out of all TVD. In order to get an impression of the joint modes the first six columns of the matrix are visualized in Figure 14.
What all TVD have in common is that there are no deformations in joint area 2, since neither the trial vectors in nor the selected vibration modes in lead to deformations in joint area 2. The first vibration mode with displacements in joint area 2 is mode number 28.
4.2.2. Static Response Computation
The static response has been computed for N and Nmm. In this computation, no friction is considered. As in the example before, the Euclidian norm of the contact force vector is used for the convergence analysis.
Figure 15 contains a convergence analysis with respect to the number of additionally used vectors. One curve contains the Euclidian norm of the contact force in case of additional joint modes based on TVD. The other curve is obtained when no joint modes are used but ordinary vibration modes. For the reference solution the contact forces are computed based on a fully nonlinear FE analysis using all nodal DOF. It can be seen that the computations with joint modes converge quicker to the reference solution and it can be assumed that the use of 10 joint modes leads to acceptable results. This is confirmed in the next subsection when dynamic investigations are performed.
As observed in the first example the function can be used as an a-priori indication of the number of necessary joint trial vectors (=g), see Figure 16.
The lower bound of suggestion (29) would lead to a matrix with 15 trial vectors which delivers acceptable results for all static and dynamic investigations.
4.2.3. Dynamic Response Computation including Contact and Friction
For the dynamic response analysis, the former loads are applied as step function to the structure. The step function is implemented as a half wave cosine with a frequency of 100 Hz.
Figures 17 and 18 contain the evolution of the Euclidian norm of the contact force vector and the friction force vector with respect to time. It can be seen in both figures that a matrix with 20 columns delivers acceptable results. This corresponds with the prior static convergence analysis. Due to the rigid body at the end of the beam and the spot welds, the friction forces and, therefore, the energy dissipation are very low.
4.2.4. Computational Efficiency
Again, all simulations are performed with the same XCos [42] block diagram (see Figure 1) with different numbers of DOF and different stiffness and mass matrices. For the multilayer sheet structure the time integration with 0, 20, and 60 additional joint modes needs 0.42%, 0.44%, and 0.49% of the CPU time which is needed for the full FE model. The balance between computation efficiency and accuracy is again excellent.
5. Discussion
As mentioned previously, the literature offers several options for the model reduction of nonlinear mechanical systems. This chapter is devoted to a comparison of the proposed method with existing methods found in the literature.
5.1. Preservation of FE DOF Which Are Involved in Nonlinear Phenomena
Some publications in the literature suggest preserving all the nodal FE DOF which are involved in nonlinearities; see [19–24]. This approach is promising with respect to accuracy but not efficient in terms of computational effort when it is applied to jointed structures because a lot of FE DOF are necessary to model the joint surfaces. A “brute force” application of this approach to the first example (friction bar) would lead to 628 additional DOF, while the proposed TVD based approach needs only 25 additional trial vectors (DOF) for a similar result quality. The gap between the proposed approach and that known from the literature increases with the complexity of the involved joints.
5.2. Proper Orthogonal Decomposition
The POD method can be directly used for the model reduction of nonlinear systems; see [25] for a review on POD. As descripted in the introduction in Section 2, the matrix includes in such a case structural displacement shapes which are obtained as solutions of the time integration of the fully nonlinear system for the time instances . POD can be used to obtain a reduction base which approximates the original space in the sense of a minimal Euclidean distance.
In comparison to the proposed approach, two disadvantages of POD are obvious.(i)In order to get a reduction base via POD at least one time integration of the fully nonlinear system is necessary. This is time consuming and the obtained reduction base is strongly related to the space spanned by the results of the time integration of the full system. In other words, applying other loads in the reduced system as they are applied in the full system may lead to considerable errors.(ii)POD minimizes a Euclidian distance. The relative displacements inside the gap are typically much smaller as the overall structural displacements which would, therefore, dominate the latter Euclidian distance. However, those relative gap displacements are significant for the nonlinear stiffness and damping due to contact and friction. From this perspective it seems questionable whether a direct application of POD is possible in case of jointed structures.
5.3. Joint Interface Modes
The author of this work has already introduced another approach for the computation of a proper reduction base for jointed structures (see [10]). In this work, a given subspace is enriched by so-called Joint Interface Modes. In principle the structure is statically condensed to the DOF of the joints. This condensation leads to a reduced mass and stiffness matrix which are used for a subsequent eigenvalue problem. The obtained eigenvectors are used as Joint Interface Modes. Gaul und Becker [29] confirmed the superior convergence rate in comparison to other methods, known from the literature. The commercially available software product MAMBA [30] is an implementation of the latter approach.
The theory of the proposed reduction base is completely different to the one found in [10] because it is a result of a power series expansion. The advantages are as follows.(i)For complex joints and structures the convergence in terms of additional trial vectors is faster. Figure 19 contains a convergence analysis of the static computation of the multilayer sheet structure. It can be seen that the proposed approach leads to trial vectors with better convergence. Two reasons can be observed by a close look at the trial vectors. The first one is that all joints in [10] are considered separately from each other. For a particular example the contact between the middle and the lower sheet is considered completely separate from the contact between the middle and the upper sheets. This is not the case for the proposed approach and this is closer to the full system because the three sheets cannot deform independently from each other. This can be seen in the left figure of Figure 20 where the third and the 19th trial vectors computed along [10] are visualized. It shows that two sheets deform but the third does not. This does not represent the construction of the beam where the spot welds enforce a coupling between all three sheets. The second reason is that the trial vectors obtained by [10] are completely independent of the initial reduction base . In [10] a subspace for all available joint DOF is computed regardless of the importance with respect to the initial trial vector base . As a consequence, joint trial vectors are computed for joint areas, which are less meaningful for a particular initial reduction base than other areas. In the second example of this work, joint area 2 is artificially added to the structure in order to illustrate this point. In fact, joint area 2 is meaningless for the problem under consideration and, therefore, no trial vector computed along the proposed method shows deformations in joint area 2. This is to be expected, because no mode in the initial reduction base leads to deformations in that area. This is not true for the joints modes computed along [10]. The right plot of Figure 20 reveals that the modes with a number higher than 19 have deformations in that area, even though it is not necessary. Note that this may have a significant impact on the convergence in case of geometrical complex joints, like those found in car bodies.(ii)The proposed algorithm is easier to implement because for the implementation of [10] a Guyan reduction [44] of all nonjoint DOF to the joint DOF is needed. This is a remarkable computational effort for large FE structures.(iii)In [10] the joint trial vectors in normal and in-plane directions of the contact surface have to be computed separately. This is not necessary for the proposed approach even if the importance of such in-plane displacements is not yet clear; see the comments on Figure 12.
6. Conclusion
In this paper a new reduction base for an accurate and efficient model reduction of jointed structures is presented. Based on an initial reduction base of a linear structure so-called trial vector derivatives are computed. In order to get an accurate and small reduction base, proper orthogonal decomposition with inner weighting has been applied to the trial vector derivatives. This transformation is not limited to jointed structures. It seems to be a promising strategy whenever trial vector derivatives lead to a high dimensional space with a lot of redundant information. As a byproduct, an a-priori estimator of the required number of joint trial vectors is given. All convergence analyses indicate that just a few of these joint trial vectors significantly improve the system response because the joint nonlinearities in terms of contact and friction can be considered accurately.
The reduction of trial vectors and, therefore, degrees of freedom has significant impact on the time that is required for the time integration. Although the quality of the result is comparable to that of a nonlinear FE computation, the reduction in CPU time is more than a factor of hundred for both examples.
In the future, applying the proposed approach to very complex structures like car bodies is planned. It is expected that a full dynamic simulation can be done considering all contact and friction forces in the joint areas all over the car body.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgment
This project was supported by the program “Regionale Wettbewerbsfähigkeit OÖ 2010–2013,” which is financed by the European Regional Development Fund and the Government of Upper Austria.