Compact thermo-mechanical models for the fast simulation of machine tools with nonlinear component behavior

2.1 Thermo-mechanical behavior of machine tools

To date, the areas of thermal [2], [3], [4] and mechanical [5], [6] simulation have been widely studied. For subsystems like the main spindle unit, interactions between both thermal and mechanical domain were investigated and respective simulation models were derived [7].

For a body of homogeneous and temporally constant material properties (density ρ, heat capacity c p and conductivity λ), the thermal behavior can be expressed analytically, resulting in a partial differential equation for the temperature field T (1). Here, the external thermal load q ˙ (head flux density) has a spatial (δ) and temporal (t) dependency, where δ denotes the three-dimensional position coordinates,

Figure 1

Closed loop of thermo-mechanical interactions in machine tools.

Figure 2

Feed axis of a modern 5-axis vertical milling machine.

Elastostatic equilibrium for solid bodies can be expressed by Navier-Cauchy equation

Again, constant material parameters (Lamé parameters λ and μ) are considered and F represents a volume-distributed force. Although we restrict our considerations to constant material parameters, the approach is not limited to this assumption. Note that in the case of variable parameters methods of MOR for parameter-varying systems, see [8], would have to be applied to (21) in Sec. 3.2.

Combined thermo-elastic behavior of solid bodies is often assumed to be linear. Although machine tools mainly consist of large structural parts, machine components like bearings, linear guideways or ball screws represent functionally relevant elements that usually act as coupling elements between the structural parts. Due to their nonlinear mechanical characteristics, thermal expansion of structural parts, for example, may affect mechanical conditions in machine components as displayed in Fig. 1.

In general, machine components act as significant heat sources in thermally sensitive regions of a machine tool. While bearings of the main spindle, for example, are characterized by the installation in structures of high power density, linear guideways induce heat along large structural parts that significantly contribute to the overall TCP error [9]. In reaction, transient temperature fields cause thermal expansion and thermal stress that ultimately result in reactive forces. These forces, in turn, influence contact pressings in machine components affecting their frictional behavior that initially caused the heat generation.

In order to investigate the effects of cross-domain behavior, an FE-based coupled model for a feed axis of a 5-axis machine tool is designed (Fig. 2).

Besides two structural parts ( P 1 and P 2) that represent the slide and the headstock, the assembly group contains two linear guideways, arranged in parallel, with two carriages each, i. e. four machine components in total ( C 1 to C 4 ). While the structural parts are modeled via FEM, the components are represented as spring damper elements with nonlinear, load-dependent characteristics. A more detailed description is given below.

2.2 Modeling of structural parts

The given thermo-elastic partial differential equations in (1) and (2) are solved numerically. By discretizing in space and time, for each structural part of the assembly, a matrix equation

with B P u t denoting the external influences on the system, can be derived.

Here, the heat capacity matrix C P only contains non-zero entries for thermal degrees of freedom (DOFs). In contrast, the matrix K P is composed of a thermal heat conduction part K th P , an elastic stiffness part K el P and a part considering thermal elongation K el , th P . Due to the different time scales of the thermal and mechanical system, dynamic effects in the mechanical part are negligible. Thus, the mechanical system is modeled as static, and also the coupling from elastic to thermal degrees of freedom is neglected, in this paper. The resulting structure of each system matrix is visualized in Fig. 3, highlighting the asymmetric structure.

Figure 3

System matrices C P 1 and K P 1 .

The machine structure is meshed with tetrahedral, quadratic SOLID227 coupled-field elements using ANSYS. As a result, for each node four DOFs

follow. The vector Θ P contains structural DOFs Δ P (translatoric displacements δ in x-, y-, and z-direction) and one thermal DOF T P (nodal temperature).

2.3 Modeling of machine components

On the one hand, linear guideways are designed to ensure backlash- and friction-free motion in feed direction (z-direction in Fig. 2). On the other hand, the other DOFs except the linear feed direction are supposed to be locked. Since dynamic effects are not considered in this paper, linear guideways are modeled as nonlinear spring elements with 7 DOFs

In addition to the considered DOFs for structural parts (4), Θ c does also contain respective rotations φ around all three coordinate axes in ϕ c .

Since the stiffness matrix

refers to the geometric center, only direct stiffness entries on the main diagonal are included.

The scalar entries K δ x δ x to K T represent mechanical elasticity and thermal heat conductivity of the respective machine component. Since these values depend on the conditions in each rolling contact inside the component, they are updated explicitly based on previous calculation results. In general, contact stiffness can be calculated analytically, with Hertzian theory [10], based on the relative positioning of the contact partners to each other. Furthermore, the equilibrium of external forces and internal contact forces of each roller element is calculated iteratively in each time step. In [11], respective analytic-numerical models for machine components like linear guideways, ball screws and bearings were developed. Based on these existing models, the nonlinear behavior of the components is calculated a priori for different operating points and provided in tables. Explicitly calculated nonlinear stiffness of guideways can be expressed as follows:

One spring element of seven DOFs connects two so-called condensation nodes of same dimension each. Thus, the resulting stiffness matrix describing one machine component has the dimension 14×14 and the structure

2.4 General assembly of the coupled system

Based on thermal and mechanical modeling of structural parts and machine components, the coupled system can be assembled. In this procedure, the following steps are performed:

1. Definition of boundary conditions (BC) and couplings between structural parts

2. Formulation of multi-point constraints (MCs)

3. Assembly of coupled system.

Furthermore, condensation node moments M C are put in equilibrium with moments coming from coupling forces F n P with corresponding lever arms Δ n P , expressed as the cross-product

In matrix notation,

provides the connection of condensation and coupling DOFs.

Based on this formulation, the global coupling matrix L for one structural part and one condensation node to be coupled can be assembled as follows:

By applying dual assembly based on theory of dynamic substructuring [12], the system represented by the uncoupled stiffness matrix K is coupled by multiplying with the matrix L from both sides. A priori, K is assembled by stacking the individual subsystems K i P and K j C diagonally,

Although forces F C and moments M C of condensation nodes, as well as coupling forces F P , were added as additional DOFs to the system of equations, only respective condensation node forces and moments remain, hereinafter referred to as Lagrange multiplier λ. The overall vector of DOFs Θ now consist of structural part DOFs Θ P , machine component DOFs Θ C and Lagrange multiplier λ,

The considered feed axis is composed of two parts and four components. The coupled stiffness matrix is

Fig. 4 shows the global stiffness matrix with its three different sub-blocks. The first, by far largest block represents part-DOFs with constant entries. The second block with condensation-DOFs contains variable entries that need to be updated every time step. However, this block only covers 56 DOFs while the assembled model has 340,034 DOFs in total (see Tab. 1). In addition, block three realizes the coupling of part and component blocks. For static assemblies without relative motion, this block is also constant. In case of a relative motion, this block needs to be updated based on the current position of the feed axis. This could be achieved through the discretization of the rails of linear guideways into discrete sections. Then, for each section a multi-point constraint has to be modeled and the corresponding matrix L MC has to be set up. As a result, a coupling block K P i C j containing connections to respective coupling nodes exists for each rail section and has to be updated position-based within the model.

Figure 4

Structure of the coupled stiffness matrix.

3.1 MOR for linear time-invariant (LTI) systems

We aim to separate the linear part of the above simulation model and reintegrate a much smaller surrogate into the coupled simulation. Hence, here we consider the representation of a linear time-invariant (LTI) system in generalized state-space form

with constant system matrices E , A ∈ R n × n , which describe the dynamics of the system, the input map B ∈ R n × m , the output map C ∈ R p × n and the feedthrough matrix D ∈ R p × m , respectively. The vector x ∈ R n denotes the state, which includes the DOFs in the considered model and the outputs y ∈ R p contain values in points of interest, like temperatures in sensor positions or the TCP-displacement. The vector u together with the input map B describe the external influences on the system, e. g., in our case, the heat fluxes induced to the linear parts from the ambient temperature or the nonlinear components of the machine.

The goal of our projection-based MOR is the computation of truncation matrices V , W ∈ R n × r , which restricts the model to an r-dimensional subspace with r ≪ n, while the essential information on the input-to-output dynamics of the system is preserved. The ROM is of the form

with the reduced matrices E ˜ : = W T E V, A ˜ : = W T A V, B ˜ : = W T B and C ˜ : = C V, the reduced state x ˜ ( t ) ∈ R r and the approximated outputs y ˜ ( t ) ≈ y ( t ) ∈ R p . The output error | | y − y ˜ | | is required to be small in a suitable norm.

There are several well known MOR techniques for LTI systems (Σ), which differ in the way the truncation matrices V and W are determined. Besides the energy-based reduction methods, like balanced truncation (BT), which we use to demonstrate our approach in this paper, also methods based on Krylov subspaces, e. g. moment matching techniques and rational interpolation, enjoy great popularity in the field of MOR. Another class of MOR methods are the training-based approaches, which rely on repeated system simulations, like proper orthogonal decomposition (POD) and reduced basis (RB) methods. More information on the various MOR techniques, especially concerning the basic principles, algorithms, stability and error estimation can be found in, e. g. [13], [14], [15], [16] and the references therein. In principle, all these methods could be used for the reduction of the linear part in our framework.

3.2 MOR for fixed relative positions

Figure 5

Sketch of the position of the feed axis model used for the simulations.

We consider the fixed pose in Fig. 5 of the feed axis model described in Section 2. In order to develop a strategy to compute a ROM for this system, we have to analyze, and make use of, the structure. After the spatial discretization by the FEM and the assembly of the coupled system (Section 2.4) we set

to obtain the following coupled system Σ in generalized state-space form

The system consists of a linear thermo-elastic part Σ lin of dimension n 1 , which is coupled to a nonlinear system Σ nl of dimension n 2 through the coupling matrices A 12 and A 21 . A 22 is the nonlinear subsystem consisting of the stiffness of each linear guide and respective condensation DOFs. Since A 22 realizes the coupling of the structural machine parts, it is responsible for relative displacements and rotations of the rigid bodies, on the one hand, and heat exchange between the bodies, on the other hand. In consequence, the impact of flexible bodies is extended by the rigid body behavior. f ( x 2 ) ∈ R n 2 2 × n 2 2 represents the nonlinear part of the stiffness matrix and equals K C j in (15), which is dependent on the DOFs of the multi-point constraints x 2 . In the simulation this part has to be updated after every time step. Thereby for every component submatrix (one submatrix for each carriage) a factor is calculated, which depends on the difference between the displacement DOFs in x-direction of the respective carriage and a product specific parameter representing the linear guideways. The matrix E 11 consists of the heat capacity matrix in the thermal part and zeros in the elastic part, which is considered to be static (see Section 2), and thus E 11 is singular.

To achieve an efficient reduction process, the linear and nonlinear parts of the system are separated, in order to treat each with the most appropriate method. This is done by a decomposition of the coupling blocks A 12 and A 21 in the form A 12 = B lin C nl and A 21 = B nl C lin , analogous to [17]. Both can be achieved by the (economy-size) singular value decomposition (SVD). For large sparse matrices, the required subset of singular values and vectors can be computed efficiently, by first determining the structural rank of the matrix and then running the sparse SVD with that number of desired values and vectors. Afterwards the system can be reformulated as an interconnected system [18] with an explicit coupling in the form of additional internal input matrices B lin , B nl and output matrices C lin , C nl :

Due to the small dimension of the nonlinear subsystem compared to that of the linear one, we decided to reduce only the linear part and preserve the full nonlinear subsystem. In case the dimension of Σ nl gets too large, suitable techniques for the reduction of nonlinear systems can be found e. g. in [14, Sec. 3], [21] and the references therein.

With a grouping of the states by thermal and elastic degrees of freedom, x th = T and x el = [ δ x δ y δ z ] T , and the combined input u el,th = [ u y nl ] T from external inputs u and internal inputs y nl from the nonlinear subsystem, we can write Σ lin as

This representation takes the form of a semi-explicit differential-algebraic equation (DAE) of index 1 [22, Ch. 5.1], and allows for an index-reduction as it was applied to a thermo-elastic machine tool assembly in [23]. From the first block-row of (18), we get

At this point it has to be guaranteed that A el has to be invertible, which will be addressed in more detail in Sec. 3.3. Inserting (20) into the second block row and the output equation (19) leads to a representation of the same system dynamics by a set of ordinary differential equations together with an updated output equation

Now defining E = E th , A = A th , x = x th , B = B th , C = C th − C el A el − 1 A el,th and D = − C el A el − 1 B el , we receive a standard LTI state-space system (Σ). If both thermal and elastic degrees of freedom use the same mesh and ansatz functions in the FE-formulation and for a three-dimensional problem, as in our case, the system dimension in (21) is already reduced to 1 4 of the state-dimension in (18), while the information of the elasticity system is completely captured in the output equation of (21). Note that for a two-dimensional problem the system dimension would be reduced to 1 3 by the index-reduction. Now, we can further reduce (21), which represents the index-reduced formulation of the linear subsystem Σ lin in (17a), using standard MOR methods for LTI systems. Finally, the resulting ROM Σ ˜ lin is recombined with the nonlinear part Σ nl . Here, special attention has to be drawn to the structure of the new feedthrough matrix D, since it also contains parts including the artificial internal inputs and outputs, which have to be considered in the newly coupled system. The specific structure of D can be split into four parts containing the feedthrough for the coupled system, but also describes the influence of the inputs u on the nonlinear states, the interaction of the linear and nonlinear states, as well as the influence of the nonlinear states on the output:

Consequently, a reduced coupled system of the form

with the updated matrices

is obtained.

3.3 Freed rigid body modes

An important fact, which has to be considered in the context of decoupling an interconnected system containing a mechanical part, is the problem of freed rigid body modes. In our case, the nonlinear part of the system represents the carriages, which connect the slide with the headstock. Thus, in the considered use case, decoupling the linear and nonlinear part of the system also leads to 6 rigid body modes in the headstock model.

Due to zero eigenvalues, A el is not invertible, but the invertibility is a prerequisit for the application of (20). In the following we want to discuss some strategies to tackle this problem.

An obvious idea is a change in the FE-model itself based on an additional fixing of the concerned assemblies. However, this could be disadvantageous due to an excessive stiffness in the FE-model and restrictions in the modeling process.

A second possibility results from a numerical perturbation of the model, leading to a shift of the eigenvalues. We follow this strategy in the numerical experiments in this contribution. To this end, we introduce an ϵ-perturbation of the form A ˆ el = A el − ϵ I in (18), analogous to [24], in our case using ϵ = 10 − 5 . In that case we have to condone the loss of error bounds of the MOR method as a drawback.

A third opportunity will be part of future investigations. Here a change from the stationary elasticity model to a dynamic second-order elasticity model including the mass matrix enables the treatment of the rigid body modes as part of MOR methods used in the context of elastic multibody simulation, see [25], [26].

Data and code availability

All MATLAB codes for execution of the experiments reported here are available as a code package authored by J. Vettermann, J. Saak and A. Steinert [28]. The data matrices for a test model are provided differing in system dimension and used parameters to the here presented model and thus generating slightly different results.

4.1 Results

Figure 6

Relative and absolute errors of the relative x-displacement between shoes and rails of the linear guideways and the related temperature difference for the ROM for the fixed position sketched in Fig. 5 simulated as fully linear system.

Figure 7

Relative and absolute errors of the relative x-displacement between shoes and rails of the linear guideways and the related temperature difference for the ROM for the fixed position sketched in Fig. 5 simulated with active nonlinearity.

Figure 8

Comparison of the results for the FOM and the ROM at the TCP for the fixed position sketched in Fig. 5 simulated with active nonlinearity.

The generation of the ROM and all simulations have been executed on a compute server, whose hardware and software specifications are described in Tab. 2, at TU Chemnitz. Since we did not have exclusive access to the machine, timings below are expected to be representative, but may have been affected by other computations on the same machine.

From a physical point of view, relevant comparison variables are the relative x-displacement between shoes and rails of the linear guideways and the absolute z-displacement of the TCP. While the former effect has significant influence on the component stiffness and thus on the absolute displacement under static load, the latter effect directly affects the achievable working accuracy.

We observe a very good conformity between the results of the original full order model (FOM) and the ROM both in the linear and nonlinear case, as Fig. 6, 7 and 8 show exemplary for the position pictured in Fig. 5. The simulation was executed with the implicit Euler method for a time span of 100 s with a constant time step size of 0.1 s, thus we have a total number of 1,000 time steps. The plots in Fig. 6 and 7 show the errors of the difference between the temperature DOFs as well as the displacement DOFs in x-direction of the carriages. A good match of these displacement values is of particular importance, due to the fact, that the update of the nonlinear part of the system depends on them. With average relative displacement errors of order 10 − 2 and thus a reduction error smaller than the expected FE-modeling error, the model accuracy is preserved. We aim at a good approximation of the displacement and temperature results in the ROM at our point of interest — the TCP. A reliable approximation especially of the TCP displacement in the ROM is crucial for, e. g. real-time error correction of the thermally induced TCP error during the manufacturing process. As Fig. 8 shows, a very good reproduction of the FOM trajectory with the ROM simulation with a relative error of order 10 − 3 can be achieved for the z-displacement, the temperature results are even better with a relative error of order 10 − 4 .

In Tab. 3 the order of the models as well as the times needed for the reduction and simulation process for the investigated simulation models are summarized. Note that we achieve a speedup of the simulation time in the fully linear case, for FOM and ROM, by the reuse of the LU-decomposition in the linear solver for the step, which is possible due to the linearity of the system and the use of a constant time step size in the integration step. This explains the large difference between the simulation times for the FOM in the linear and nonlinear cases.

We see that the simulation of the ROM comes up with a formidable speedup of the simulation times compared to the FOM. Even in the case that just one simulation is required, the sum of the time needed for the reduction and simulation process of the ROM is smaller than one simulation cycle of the FOM. Furthermore, a comparison of the storage amount of the FOM and ROM might be an important fact, which plays a major role if a model has to be saved directly on the numerical control system of the machine tool. We can not only observe a drastical reduction of the simulation time, but also of the memory requirements for the ROM. While the FOM needs 212.81 MB, the ROM just requires 0.23 MB.