# Efficient model order reduction for the dynamics of nonlinear multilayer sheet structures with trial vector derivatives.

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 [A.sub.joint], which is identical to the spatial distribution of the joint. However, due to small structural displacements there are areas Ac inside the joint where the surfaces are in contact and other areas [A.sub.g] where the surfaces are not in contact. Consider

[A.sub.Joint] = [A.sub.c] (x(t)) + [A.sub.g] (x(t)). (1)

The argument x(t) 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  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  and friction  phenomenon, the FEM leads to a coupled and nonlinear differential equation

[M.bar][??] + [K.bar]x + [f.sup.NL.sub.(x)] = [f.sub.(t)], (2)

where the (n x 1) vector x contains the system's degrees of freedom (DOF). Note, no material nor any viscous damping is considered in (2). The constant (n x n) matrices [M.bar] and [K.bar] represent the structural mass and stiffness. The (n x 1) vector [f.sub.(t)] holds the external and imposed forces while the (n x 1) vector [f.sup.NL.sub.(x)] 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 n, which can be quite large in the case of relevant finite element (FE) models found in the industry. Such FE models may have [10.sup.7] 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 x is approximated by [??] which is limited to a scaled superposition of r time invariant trial vectors. Formally, this can be given as

x [approximately equal to] [??] = [r.summation over (i = 1)][[phi].sub.i][q.sub.i] = [[PHI].bar]q, (3)

where the (n x 1) vector [[phi].sub.i] holds the time invariant ith trial vector and [q.sub.i] the corresponding time varying scaling factor. All trial vectors can be collected as columns in the (n x r) matrix [[PHI].bar] and the scaling factors are collected in the (r x 1) column vector q. The application of (3) on the equation of motion (2) leads to r instead of n DOF, which will be shown in the next chapter in greater detail. Good reduction approaches provide a number of r which is much smaller than n with a small error

e = [bar.x - [??]], (4)

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 [[PHI].bar] 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 [[PHI].bar] for linear systems. However, in nonlinear systems it is more challenging to find a proper subspace [[PHI].bar]. 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 [[PHI].bar] (see  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 [x.sub.1] to [x.sub.m] which are obtained by time integration of the full system for the time steps [t.sub.1] to [t.sub.m]. 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.  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.  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  some general comments are made on this issue and  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 [[PHI].bar] (see ). 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  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 . Segalman  suggested at the same time an approach for localized nonlinearities, which is a nonjoint focused generalization of . A comparison of the proposed method with the one of  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 [x.sub.L] and a time varying portion [DELTA]x, such that

x = [x.sub.L] + [DELTA]x. (5)

The vector [x.sub.L] can be interpreted as an operating point and the vector [DELTA]x as the vibrations around this point. In the very common case that a structure vibrates around its undeformed state, the vector [x.sub.L] is simply the zero vector. Examples where [x.sub.L] is not trivial are prestressed screwed joints or leaf springs. As a consequence of the latter decomposition the stiffness matrix [K.bar] simply holds the stiffness at [x.sub.L]. Rearranging (2) yields

[M.bar][??] + [K.bar]x = f(t) - [f.sup.NL.sub.(x)] (6)

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 r trial vectors. Inserting equation (3) into (5) and premultiplying with [[[PHI].bar].sup.T] delivers

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (7)

with the (r x r) reduced mass matrix [??] = [[PHI].bar].sup.T][M[PHI].bar] and the (rxr) reduced stiffness matrix [??] = [[PHI].sup.T][K[PHI].bar]. The used subspace [[PHI].bar] needs to fulfill two criteria which influence each other. The subspace [[PHI].bar] should provide an approximated displacement x so that accurate nonlinearities can be computed and, as a consequence, the resulting error e (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

[[PHI].bar] = [[[[PHI].bar].sub.V] [[[PHI].bar].sub.S]], (8)

where the (n x v) matrix [[[PHI].bar].sub.V] contains v vibration mode shapes and the (n x s) matrix Os contains [[[PHI].bar].sub.S] trial vectors which are displacement shapes due to static loads. These s static displacement shapes are typically related to s nonzero entries in the external force vector f(t). According to the pseudoforce concept motivated by the structure of (6), one could suggest preserving all DOF which are involved in the vector [f.sup.NL.sub.(x)], 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 [f.sup.NL.sub.(x)]. In that case the resulting quality cannot be guaranteed without a convergence analysis with respect to the number of vibration mode shapes collected in [[[PHI].bar].sub.V]. It is shown below that this approach leads to a poor convergence, which means that the number of v is high.

2.3. Trial Vector Derivatives. The idea is to extend the trial vector base (8) in the form of

[[PHI].bar] = [[[[[PHI].bar].sub.V] [[[PHI].bar].sub.S] [[[PHI].bar].sub.T]]], (9)

where the (nxt) matrix [[[PHI].bar].sub.T] contains the so-called trial vector derivatives (TVD). The TVD approximate the difference between [[[PHI].bar].sub.V] and [[[PHI].bar].sub.S] computed with a nonlinear stiffness matrix with respect to [[[PHI].bar].sub.V] and [[[PHI].bar].sub.S] computed with a linearized stiffness matrix (8). Therefore, the computation of [[[PHI].bar].sub.T] can be performed using a finite difference scheme as Slaats et al. proposed in .

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 [q.sub.1] to [q.sub.r]. The state dependent trial vector [[phi].sub.i](i [less than or equal to] (v + s)) can be expressed by a Taylor series expansion as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (10)

The (n x 1) vector [[phi].sub.i][|.sub.q = 0] holds the ith trial vector computed at the point where the structure is linearized; see (5). This vector is already part of [[[PHI].bar].sub.V] and [[[PHI].bar].sub.S]. Note that the expression [[partial derivative][[phi].sub.i]/[partial derivative]q] holds an n by (v + s) 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

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (11)

The n by [(v + s).sup.2] matrix [[[PHI].bar].sub.T] contains the first order terms of (10) applied to all trial vectors of [[[PHI].bar].sub.V] and [[[PHI].bar].sub.S]. That means, the trial vector base is increased by its squared column dimension and the number of generalized coordinates collected in the vector q 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 [[[PHI].bar].sub.T]; see . The space spanned by [[[PHI].bar].sub.T] can be approximated by a reduced space [[[??].bar].sub.T] which contains most of the important information. This issue is addressed below, where a POD based transformation (reduction) of [[[PHI].bar].sub.T] into [[[??].bar].sub.T] 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 i. A derivation of both sides with respect to the scaling factor of mode number j gives

[partial derivative]/[partial derivative][q.sub.j] [([K.bar] - [[OMEGA].sup.2.sub.i][M.bar]])[[phi].sub.i]] = 0. (12)

The application of the product rule and the assumption that the mass matrix is constant lead to

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (13)

where the (n x 1) vector [partial derivative][[phi].sub.i]/[partial derivative][q.sub.j] holds the TVD of mode i with respect to mode j and is hereon denoted as [[phi].sub.i,j]. When (13) is premultiplied with [[phi].sup.T.sub.i] the sensitivity of the eigenvalue [[OMEGA].sup.2.sub.i] with respect to [q.sub.j] can be computed (see ). Note that the symmetry of [K.bar] and [M.bar] has to be considered in this step along with the mass and stiffness orthogonality of the modes. Consider

[partial derivative][[OMEGA].sup.2.sub.i]/[partial derivative][q.sub.j] = [[phi].sup.T.sub.i][[partial derivative][K.bar]/[partial derivative][q.sub.j]][[phi].sub.i]. (14)

Substitution of the according term in (13) with (14) delivers the TVD:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (15)

It has been observed in the literature  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

[[phi].sub.i,j] = [[K.bar].sup.-1][[partial derivative][K.bar]/[partial derivative][q.sub.j]][[phi].sub.i]. (16)

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 j gives

[partial derivative]/[partial derivative][q.sub.j]([K.bar][[phi].sub.i]) = [[partial derivative]/[partial derivative][q.sub.j]][[phi].bar] (17)

Considering the product rule and the fact that [f.sub.i] is state independent delivers

[[phi].sub.i,j] = [[K.bar].sup.-1] [[partial derivative][K.bar]/[partial derivative][q.sub.j]] [[phi].sub.i] (18)

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

[f.sub.i,j] = [[partial derivative][K.bar]/[partial derivative][q.sub.j]][[phi].sub.i]. (19)

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 maybe the computation of [partial derivative][K.bar]/[partial derivative][q.sub.j]. It is discussed in Section 3 how [partial derivative][K.bar]/[partial derivative][q.sub.j] can be computed when nonlinearities due to contact and friction inside a joint are considered.
```
Algorithm 1: Computation of TVD for a given reduction base
[[phi].sub.1] to [[phi].sub.[upsilon] + s].

Input: [[phi].sub.1],
..., [[phi].sub.i], ...,[[phi].sub.[upsilon]+s],[K.bar]

Output:
[[[PHI].bar].sub.T] = [[[phi].sub.1,1],
...,[[phi].sub.1,[upsilon]+s], ...,[[phi].sub.i,1], ...,
[[phi].sub.1,[upsilon]+s], ..., [[phi].sub.1,[upsilon]+s,1], ...,
[[phi].sub.1,[upsilon]+s+[epsilon]+s]]

(1) [F.bar] = [] (2) for j = 1 : (v + s) (3) determine [partial
derivative]K/[partial derivative][q.sub.j] (4) for i = 1 : (v + s)
(5) [f.sub.i,j] = [partial derivative][K.bar]/[partial
derivative][[phi].sub.i] (6) [F.bar] = [[F.bar][f.sub.i,j]] (7) end
(8) end (9) use FE software to compute [bar.[PHI]]T based on
[[K[PHI].bar].sub.V], = [F.bar]
```

2.6. Comment on a Singular Stiffness Matrix. As outlined in the introduction, the literature offers different approaches for a proper reduction base [[PHI].bar]. 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 ). 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 [(v + s).sup.2] 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 [[[PHI].bar].sub.T] 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 [[[??].bar].sub.T] which contains most of the information of [[[phi].bar].sub.T]. For some computations in  just the [[PHI].bar.i,i] 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  a strategy is presented, which is based on the idea that when the according modes [[phi].sub.i] and [[phi].sub.j]. are important, the TVD [[phi].sub.i,j] is important as well. In the present paper a general strategy for the determination of [[[??].bar].sub.T] based on [[PHI].bar] 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 . For the sake of readability, the main points of POD are briefly reviewed next.

Let us assume p independent (n x 1) vectors [y.sub.1] to [y.sub.p] which are gathered in the (n x p) matrix [Y.bar]. POD of rank g delivers g orthonormal (m x i) vectors [u.sub.1] to [u.sub.g], which approximate the space spanned by [bar.Y] optimal in a Euclidean sense. The vectors [u.sub.1] to [u.sub.g] are named proper orthogonal modes (POM) and can

be gathered in the (n x g) matrix [U.bar]. The POMs maximize a function J in the form of

J([u.sub.1], ...,[u.sub.g]) = [g.summation over (i=1)][p.summation over (j=0)][([y.sup.T.sub.j][u.sub.i]).sup.2] [right arrow] max (20)

so that

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (21)

(see  for a proof). The POM [u.sub.1] to [u.sub.g] can be found by using the first g eigenvectors of

[[Y.bar].sup.T][Y.bar][[??].sub.i] = [[lambda].sub.i][[??].sub.i] for i = 1 ... g (22)

and a subsequent transformation

[u.sub.i] = 1/[square root of [[lambda].sub.i]][Y.bar][[??].sub.i] for i = 1 ... g (23)

A slightly different POD method with a weighted inner product maximizes a function / in the form of

J ([u.sub.1], ...,[u.sub.g]) = [g.summation over (i=1)][p.summation over (j=1)][([y.sup.T.sub.j][A.bar][u.sub.i]).sup.2] [right arrow] max (24)

so that condition (21) is fulfilled as well. The (n x n) matrix A has to be symmetric and positive semidefinite. The POM for the weighted POD can be found by using the first g eigenvectors of

[[Y.bar].sup.T][AY.bar][[??].sub.i] = [[lambda].sub.i][[??].sub.i] for i = 1 ... g (25)

and a subsequent transformation along equation (23). The second POD variant can be transformed into the first one if the (n x n) identity matrix is substituted for the matrix [A.bar]. From a mechanical point of view, the characteristics of the matrix [A.bar] 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 g trial vectors which optimally approximate the deformation energy caused

by the space [Y.bar]. The POM without a weighted inner product approximates the displacement space only. The eigenvalue [[lambda].sub.i] is of special interest because it is related to the importance of the vector [u.sub.i] for the minimization of (20) or 24). Note that large eigenvalues indicate important POM. Therefore, a ratio w(g) is defined which can be used as a criterion to determine the number of g. Consider

w(g) = [[lambda].sub.1] + ... + [[lambda].sub.g]/[[lambda].sub.1] + ... + [[lambda].sub.g] + ... + [[lambda].sub.p]. (26)

Note that w(0) = 0 and w(p) = 1. In case of POD with weighted inner product, w(g) provides a value of how much of the available deformation energy is already covered by [u.sub.1] to [u.sub.g]. 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 T VD, the following substitutions need to be performed. Instead of the matrix [Y.bar] the (n x[(v + s).sup.2]) matrix [[[PHI].bar].sub.T] has to be used. It is important to mention that all column vectors of [[[PHI].bar].sub.T] are normalized with respect to the stiffness matrix [K.bar]. 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 [U.bar] corresponds with the one spanned by the (n x g) matrix [[[??].bar].sub.T]. In case of inner weighting the stiffness matrix [K.bar] is used instead of [A.bar].

The finally obtained reduction space is [[PHI].bar] = [[[[PHI].bar].sub.V] [[[PHI].bar].sub.S] [[[PHI].bar].sub.T]].

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 (n x n) matrix [partial derivative][K.bar]/[partial derivative][q.sub.j]. The qualitative content of the latter matrix is the change in stiffness due to a variation induced by trial vector number j. 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, w gap functions [g.sub.1] ... [g.sub.w] can be defined. The gap function [g.sub.i] represents the normal distance of the node-pair i. Following a simple penalty approach, a constant contact stiffness [k.sub.n] and a constant tangential stiffness [k.sub.t] are acting between those node-pairs with a gap function gt smaller than 0 (penetration). If [g.sub.i] [greater than or equal to] 0 (gaping) no stiffness acts between the node-pairs. We assume that the linear modes [[[PHI].bar].sub.V] and [[[PHI].bar].sub.S] 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 [partial derivative][K.bar]/ [[[PHI].bar].sub.T][q.sub.j] contains the stiffnesses [k.sub.n] and [k.sub.t] which act only on those DOF corresponding to nodepairs with a gap function less than zero. The gap functions are evaluated based on the displacements due to the trial vector [[phi].sub.j].

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

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (27)

Note that the superscript "E" denotes an effective stiffness matrix which consists of the linear portion [K.bar] 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

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (28)

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 [[phi].sub.j]. 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 [[phi].sub.i] with respect to trial vector number [[phi].sub.j], namely, [[phi].sub.i,j] and [[phi].sub.i,-j] The (n x 1) vector [[phi].sub.i,-j]. denotes the trial vector derivative based on a negative scaling of trial vector [[phi].sub.j]. In case of joint induced nonlinearities the matrix [[PHI].bar] is of size (n x 2[(v + s).sup.2]). 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 g 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 [[[??].bar].sub.T].

For both examples, the matrices [[[phi].bar].sub.V] and [[[phi].bar].sub.S] were com puted according to the Craig/Bampton method . For this particular method, the content of the matrix [[[phi].bar].sub.V] is named "Fixed Boundary Normal Modes" and the content of the matrix [[[phi].bar].sub.S] is named "Constrained Modes." In both examples [[[phi].bar].sub.V] and [[[phi].bar].sub.S] 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 . The matrices were imported into Scilab  and all subsequent computations have been done directly in Scilab or in its toolbox Xcos .
```
Algorithm 2: Construction of reduction base for jointed structures.

Input: FE-Model, M, K
Output: [[PHI].bar] = [[[[PHI].bar].sub.S]
[[[PHI].bar].sub.V] [[[??].bar].sub.T]]
%%Computation of trial vectors for the linear system
(1) Compute s trial vectors based on static deflection shapes
[[[PHI].bar].sub.V]
(2) Compute v trial vectors based on vibration modes
[[[PHI].bar].sub.V]
%%Computation of trial vector derivatives
(3) [F.bar] = []
(4) for j = 1 : (v + s)
(5)       for k = -1:2:1
(6)                  [partial derivative][K.bar]/[partial
derivative][q.sub.k.j] = zeros(", n)
(7)                  for z = 1 : w
(8)                         compute g(z) due to k[[phi].sub.i]
(9)                         if (0(z) <= 0)
(10)                                 add [k.sub.n] and [k.sub.t] to
[partial derivative][K.bar]/
[partial derivative][q.sub.k.j]
(corresponding to DOF of FE
node pair z)
(11)                         end
(12)              end
(13)              for i = 1 : (v + s)
(14)                        [f.sub.i,k x j] = [[partial derivative]
[K.bar]/[partial derivative][q.sub.k x
j]][[phi].sub.i]
(15)                        [F.bar] = [[F.bar] [f.sub.ik x j]]
(16)              end

(17)      end
(18) end
(19) use FE software to compute [[[PHI].bar].sub.T] based on
[[K[PHI].bar].sub.T], = [F.bar]
(20)normalize each column of [[[PHI].bar].sub.T] with respect to
[K.bar]
%%POD
(21) if (POD with inner weighting)
(22)      solve eigenvalue problem [[[PHI].bar].sub.T.sup.T]
[[K[PHI].bar].sub.V][[phi].sub.T,i] =
[[lambda].sub.i][[phi].sub.T,i]
(23) else
(24)      solve eigenvalue
problem [[[PHI].bar].sub.T.sup.T][[[PHI].bar].sub.T]
[[phi].sub.T,i] = [[lambda].sub.i][[phi].sub.T,i]
(25) determine g
(26) [[[??].bar].sub.V] = []
(27) for i = 1 : g
(28)         [[[??].bar].sub.T] = [[[[??].bar].sub.T] [1/[square
root of [[lambda].sub.i]] [[[PHI].bar].sub.T]
[[phi].sub.T,i]]
(29) end
```

<169_TB001>

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 [10.sup.4] N/mm. The friction forces, if necessary, were computed as given in . The latter model requires two parameters, namely, the friction coefficient [mu] and an in-plane stiffness when a node-pair is in contact. For the friction coefficient [mu] a value of 0.25 is used and the in-plane stiffness is set to [10.sup.4] N/mm. The contact and friction forces are computed for all involved nodes pairs and the force vector ([f.sup.NL.sub.(x)]) is constructed according to the physical DOF. This physical force vector is then projected into the reduction space by a premultiplication with the matrix [[[PHI].bar].sup.T].

For the time integration the Scilab  toolbox Xcos  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 x 25 mm x 0.5 mm are connected via two beams. The extension of the joint is 85 mm x 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 y direction ([F.sub.y]) and a torque around the x axis ([M.sub.x]). 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 . The six DOF of the free end are defined as so-called "boundary DOF." This leads to six static deflection shapes for the matrix [[[PHI].bar].sub.S]. The two most relevant "Constraint Modes" (deflection along y and rotation around x) are visualized in Figure 3.

For the matrix [[[PHI].bar].sub.V]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 [[[PHI].bar].sub.V]. 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 [[phi].sub.1,2] and [[phi].sub.l,4] or [[phi].sub.1,3] and [[phi].sub.1,5] 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 [[[??].bar].sub.T]. 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 [F.sub.y] = -0.7 N and [M.sub.x] = -20 Nmm. For this computation, no friction is considered. For the convergence analysis, a (1 x w) vector [f.sub.c] is constructed. The vector contains w entries, namely, [f.sub.c,1] to [f.sub.c,w] and [f.sub.c,i] (1 [less than or equal to] i [less than or equal to] w) maintains the contact force at FE node-pair i. As convergence criteria the Euclidean norm of the vector fc is evaluated. Figure 7 contains a convergence analysis with respect to the number of considered joint vectors in the matrix [[[??].bar].sub.V] (=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 w(g) for both approaches. The scalar value g represents the number of trial vectors finally considered in [[[??].bar].sub.t].

A comparison of Figures 7 and 8 indicates that in case of weighted POD, the w(g) 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 w(g) = 0.9994 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 g such that

0.95 [less than or equal to] w(g) [less than or equal to] 0.9999. (29)

The actual limit for the function w(g) 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 [f.sub.c] with respect to time. It can be seen that a matrix [[[??].bar].sub.t] 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 [f.sub.c] 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 . 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 (2w x 1) friction force vector [f.sub.R] is evaluated with respect to time. The vector [f.sub.R] contains friction forces in the x and z directions of the w 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 [partial derivative][K.bar]/[partial derivative][q.sub.j] 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  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 x 25 mm x 0.5 mm are connected via three beams. Due to its construction the structure contains two joints with an extension of 135 mm x 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 y direction (Fy) and a torque around the x axis ([M.sub.x]). 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 . 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 . The six DOF of the free end are defined as socalled "boundary DOF." Consequently, the matrix [[[PHI].bar].sub.S] which holds the static deflections shapes contains six trial vectors. The number of vibration modes which are stored in the matrix [[[PHI].bar].sub.V] is 10.

4.2.1. Trial Vector Derivatives. The trial vector derivatives [[[??].bar].sub.T] and the final used joint modes [[[??].bar].sub.V] 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 [[[??].bar].sub.T] 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 O[[[PHI].bar].sub.S] nor the selected vibration modes in [[[PHI].bar].sub.V] 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 [F.sub.y] = -50 N and [M.sub.x] = -2500 Nmm. In this computation, no friction is considered. As in the example before, the Euclidian norm of the contact force vector [f.sub.c] 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 w(g) 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 [[[??].bar].sub.T] 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 [f.sub.c] and the friction force vector [f.sub.R] with respect to time. It can be seen in both figures that a matrix [[[??].bar].sub.T] 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  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 respecttoaccuracybut notefficient in termsofcomputational 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  for a review on POD. As descripted in the introduction in Section 2, the matrix [Y.bar] includes in such a case z structural displacement shapes which are obtained as solutions of the time integration of the fully nonlinear system for the time instances 0, [t.sub.1], [t.sub.2], ..., [t.sub.z-1]. POD can be used to obtain a reduction base which approximates the original space (= Y) 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 ). 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  confirmed the superior convergence rate in comparison to other methods, known from the literature. The commercially available software product MAMBA  is an implementation of the latter approach.

The theory of the proposed reduction base is completely different to the one found in  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  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  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  are completely independent of the initial reduction base [[[[PHI].bar].sub.V] [[[PHI].bar].sub.S]] .In  a subspace for all available joint DOF is computed regardless of the importance with respect to the initial trial vector base [[[[PHI].bar].sub.V] [[[PHI].bar].sub.S]]. As a consequence, joint trial vectors are computed for joint areas, which are less meaningful for a particular initial reduction base [[[[PHI].bar].sub.V] [[[PHI].bar].sub.S]] 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 [[[[PHI].bar].sub.V] [[[PHI].bar].sub.S]] leads to deformations in that area. This is not true for the joints modes computed along . 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  a Guyan reduction  of all nonjoint DOF to the joint DOF is needed. This is a remarkable computational effort for large FE structures.

(iii) In  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 socalled 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.

http://dx.doi.org/10.1155/2014/913136

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 Wettbewerbsfahigkeit OO 2010-2013," which is financed by the European Regional Development Fund and the Government of Upper Austria.

References

 L. Gaul and R. Nitsche, "The role of friction in mechanical joints," Applied Mechanics Reviews, vol. 54, no. 2, pp. 93-105, 2001.

 L. Gaul and J. Lenz, "Nonlinear dynamics of structures assembled by bolted joints," Acta Mechanica,vol. 125, no. 1-4, pp. 169181, 1997.

 W. D. Iwan, "A distributed-element model for hysteresis and its steady-state dynamic response," ASME Journal of Applied Mechanics, vol. 33, pp. 893-900, 1966.

 D. D. Quinn and D. J. Segalman, "Using series-series Iwan-type models for understanding joint dynamics," Journal of Applied Mechanics, Transactions ASME, vol. 72, no. 5, pp. 666-673, 2005.

 Y. Song, C. J. Hartwigsen, D. M. McFarland, A. F. Vakakis, and L. A. Bergman, "Simulation of dynamics of beam structures with bolted joints using adjusted Iwan beam elements," Journal of Sound and Vibration, vol. 273, no. 1-2, pp. 249-276, 2004.

 D. J. Segalman, "Modelling joint friction in structural dynamics," Structural Control and Health Monitoring, vol. 13, no. 1, pp. 430-453, 2006.

 Y. Song, D. M. McFarland, L. A. Bergman, and A. F. Vakakis, "Effect of pressure distribution on energy dissipation in a mechanical lap joint," AIAA Journal, vol. 43, no. 2, pp. 420-425, 2005.

 W. Witteveen, H. Irschik, H. Riener, M. Engelbrechtsmiiller, and A. Plank, "An efficient mode based approach for the dynamic analysis of jointed and local damped structures: joint interface modes," in Proceedings of International Conference on Noise and Vibration Engineering (ISMA'08), pp. 1815-1824, Leuven, Belgium, 2008.

 M. Breitfuss and W. Witteveen, "Efficient mode based approach for the dynamics of jointed FE structures," in Proceedings of the 8th International Conference on Structural Dynamics (EURODYN '11), G. De Roeck, G. Degrande, G. Lombaert, and G. Muller, Eds., Leuven, Belgium, July 2011.

 W. Witteveen and H. Irschik, "Efficient mode-based computational approach for jointed structures: joint interface modes," AIAA Journal, vol. 47, no. 1, pp. 252-263, 2009.

 T. A. Laursen, Computational Contact and Impact Mechanics, Chapter 3. 3. 2, Springer, Berlin, Germany, 2002.

 T. A. Laursen, Computational Contact and Impact Mechanics, Chapter 3. 4. 3.1, Springer, Berlin, Germany, 2002.

 A. K. Noor, "Recent advances and applications of reduction methods," Applied Mechanics Reviews, vol. 47, no. 5, pp. 125-146, 1994.

 R. J. Craig, "A review of time-domain and frequency-domain component mode synthesis methods," The International Journal of Analytical and Experimental Modal Analysis, vol. 2, no. 2, pp. 59-72, 1987.

 Q. Q. Zu, Model Order Reduction Techniques, Springer, London, UK, 2004.

 P. Koutsovasilis and M. Beitelschmidt, "Comparison of model reduction techniques for large mechanical systems: AA study on an elastic rod," Multibody System Dynamics, vol. 20, no. 2, pp. 111-128, 2008.

 W. Witteveen, "On the modal and non-modal model reduction of metallic structures with variable boundary conditions," World Journal of Mechanics, vol. 2, no. 6, pp. 311-324, 2012.

 B. Besselink, U. Tabak, A. Lutowska et al., "A comparison of model reduction techniques from structural dynamics, numerical mathematics and systems and control," Journal of Sound and Vibration, vol. 332, pp. 4403-4422, 2013.

 Z.-Q. Qu, "Model reduction for dynamical systems with local nonlinearities," AIAA Journal, vol. 40, no. 2, pp. 327-333, 2002.

 M. I. Friswell, J. E. T. Penny, and S. D. Garvey, "Using linear model reduction to investigate the dynamics of structures with local non-linearities," Mechanical Systems and Signal Processing, vol. 9, no. 3, pp. 317-328, 1995.

 T. J. Royston and R. Singh, "Periodic response of mechanical systems with local non-linearities using an enhanced Galerkin technique," Journal of Sound and Vibration, vol. 194, no. 2, pp. 243-263, 1996.

 R. H. B. Fey, D. H. Van Campen, and A. De Kraker, "Long term structural dynamics of mechanical systems with local nonlinearities," Journal of Vibration and Acoustics, Transactions of the ASME, vol. 118, no. 2, pp. 147-153, 1996.

 J. H. Gordis and J. Radwick, "Efficient transient analysis for large locally nonlinear structures," Shock and Vibration, vol. 6, no. 1, pp. 1-9, 1999.

 P. Hagedorn and W. Schramm, "On the dynamics of large systems with localized nonlinearities," Journal of Applied Mechanics, Transactions ASME, vol. 55, no. 4, pp. 946-951, 1988.

 G. Kerschen, J.-C. Golinval, A. F. Vakakis, and L. A. Bergman, "The method of proper orthogonal decomposition for dynamical characterization and order reduction of mechanical systems: an overview," Nonlinear Dynamics, vol. 41, no. 1-3, pp. 147-169, 2005.

 P. M. A. Slaats, J. de Jongh, and A. A. H. J. Sauren, "Model reduction tools for nonlinear structural dynamics," Computers and Structures, vol. 54, no. 6, pp. 1155-1171, 1995.

 P. Tiso, E. Jansen, and M. Abdalla, "Reduction method for finite element nonlinear dynamic analysis of shells," AIAA Journal, vol. 49, no. 10, pp. 2295-2304, 2011.

 P. Tiso, "Optimal second order reduction basis selection for nonlinear transient analysis," in Modal Analysis Topics Volume 3, Conference Proceedings of the Society for Experimental Mechanics Series, pp. 27-39, Springer, NewYork, NY, USA, 2011.

 L. Gaul and J. Becker, "Damping prediction of structures with bolted joints," Shock and Vibration, vol. 17, no. 4-5, pp. 359-371, 2010.

 http://mamba.ecs.steyr.com/.

 D. J. Segalman, "Model reduction of systems with localized nonlinearities," Journal of Computational and Nonlinear Dynamics, vol. 2, no. 3, pp. 249-266, 2007.

 J. A. Stricklin and W. E. Haisler, "Formulations and solution procedures for nonlinear structural analysis," Computers and Structures, vol. 7, no. 1, pp. 125-136, 1977

 R. Villaverde and M. M. Hanna, "Efficient mode superposition algorithm for seismic analysis of non-linear structures," Earthquake Engineering & Structural Dynamics, vol. 21, no. 10, pp. 849-858, 1992.

 D. De Klerk, D. J. Rixen, and S. N. Voormeeren, "General framework for dynamic substructuring: history, review, and classification of techniques," AIAA Journal, vol. 46, no. 5, pp. 1169-1181, 2008.

 M. I. Friswell and J. E. Mottershead, Finite Element Model Updating in Structural Dynamics, Kluwer Academic Publishers, 1995.

 2014, http://help.scilab.org/docs/5.4.1/fr_FR/CVode.html.

 A. Chatterjee, "An introduction to the proper orthogonal decomposition," Current Science, vol. 78, no. 7, pp. 808-817, 2000.

 S. Volkwein, "Model Reduction Using Proper Orthogonal Decomposition," 2013, http://www.uni-graz.at/imawww/volkwein/POD.pdf.

 R. R. Craig and M. C. Bampton, "Coupling of substructures for dynamic analysis," AIAA Journal, vol. 6, no. 7, pp. 1313-1319, 1968.

 M. S. C. NASTRAN, Version 2011. 1. 0, http://simcompanion .mscsoftware.com.

 Scilab, "Version 5. 4. 0," http://www.scilab.org/.

 Scilab, "Version 5. 4. 0," http://www.scilab.org/scilab/features/ xcos.

 C. Garcia Garino and J.-P. Ponthot, "A quasi-coulomb model for frictional contact interfaces. Application to metal forming simulations," Latin American Applied Research, vol. 38, no. 2, pp. 95-104, 2008.

 R. J. Guyan, "Reduction of stiffness and mass matrix," AIAA Journal, vol. 3, no. 2, p. 380, 1965.

Wolfgang Witteveen and Florian Pichler

Department of Mechanical Engineering, Upper Austria University of Applied Sciences, 4600 Wels, Austria

Correspondence should be addressed to Wolfgang Witteveen; wolfgang.witteveen@fh-wels.at

Received 17 September 2013; Revised 19 February 2014; Accepted 11 March 2014; Published 29 April 2014

```
Table 1: Solver parameter for time integration.

Solver kind                      Sundials/CVODE-
BDF-FUNCTIONAL


Integrator absolute tolerance         1e-6
Integrator relative tolerance         1e-8
Tolerance on time                     1e-10
```
Title Annotation: Printer friendly Cite/link Email Feedback Research Article Witteveen, Wolfgang; Pichler, Florian Shock and Vibration Report Jan 1, 2014 11033 Effect of shear on ultrasonic flow measurement using nonaxisymmetric wave modes. Dynamic mechanical properties and constitutive relation of an aluminized polymer bonded explosive at low temperatures. Derivatives (Mathematics) Dynamics Dynamics (Mechanics) Leaf springs Vector spaces Vectors (Mathematics)