Model order reduction of thermo-mechanical models with parametric convective boundary conditions: focus on machine tool
Pablo Hernández-Becerro Daniel Spescha Konrad Wegener
math.NA cs.CE cs.NA









mymainaddress,mysecondaryaddress]Pablo Hernández-Becerromycorrespondingauthor
[mycorrespondingauthor]Corresponding author
hernandez@inspire.ethz.ch
mymainaddress]Daniel Spescha


mysecondaryaddress]Konrad Wegener



[mymainaddress]inspire AG, Technoparkstrasse 1, 8005 Zurich, Switzerland
[mysecondaryaddress]Institute of Machine Tools and Manufacturing (IWF), ETH Zurich, Leonhardstrasse 21, 8092 Zurich,  Switzerland



This paper presents a parametric Model Order Reduction (MOR) method for weakly coupled thermo-mechanical Finite Element (FE) models of machine tools and other similar mechatronic systems. This work proposes a reduction method, Krylov Modal Subspace (KMS), and a theoretical bound of the reduction error. The developed method addresses the parametric dependency of the convective boundary conditions using the concept of system bilinearization. Additionally, this paper investigates the coupling between the reduced-order thermal system and the mechanical response. A numerical example shows that the reduced-order model captures the response of the original system in the frequency range of interest.




Parametric Model Order Reductionthermo-mechanical modelmachine tools








§ INTRODUCTION


Machine tools, such as milling machines, grinding machines, or lathes, are complex mechatronic systems and key components in the manufacturing process. Increasing the precision of machine tools enables to manufacture more accurate parts required in many engineering applications. Among the different error sources limiting the precision of machine tools, thermally induced deviations are the main contributor to geometric errors in manufactured parts, as stated by Mayr et al. <cit.>. During the manufacturing process, the temperature distribution of machine tools is inhomogeneous and varies over time, leading to deformations of the structural parts. Internal and external sources are responsible for the time-varying temperature distribution, as summarized by Bryan <cit.>. Internal sources refer to heat losses at the machine elements during the manufacturing process, such as friction at the bearings. External sources refer to the surroundings of the machine tool, such as fluctuations of the environmental temperature of the workshop.

Virtual prototypes are a great asset to improve the thermo-mechanical behavior of machine tools, as explained by Altintas et al. <cit.>. They enable to test different design alternatives virtually before a physical prototype is available. Virtual prototypes are numerical models that describe the physical processes leading to errors in manufactured parts. Thermo-mechanical models of machine tools are based on the finite element (FE) discretization of the heat transfer and elasticity equations. However, discretized thermo-mechanical models of machine tools have a large scale due their geometrical complexity, requiring a high computational effort. Many engineering applications use surrogate models, a computationally efficient model reproducing the physical behavior of the high-fidelity model. Among the different techniques to create surrogate models, this work focuses on projection based model order reduction (MOR). MOR projects the original model into a lower dimensional subspace containing the most relevant information of the response of the system. As stated by Benner et al. <cit.>, the main benefit of MOR is that it retains the systems structure enabling the traceability of the dynamical evolution. Additionally, the use of the underlying system structure facilitates the derivation of theoretical error bounds of the surrogate model. In comparison to other surrogate modeling techniques, MOR is an intrusive method, i.e. it needs to access and modify the system matrices. 

Several works apply MOR to thermo-mechanical models of machine tools and other mechatronic systems. Herzog et al. <cit.> presented a thermo-mechanical model of a machine tool column. They concentrated on the problem of optimal sensor placement, i.e. selecting the optimal points to place temperature sensors in order to predict the thermally induced deviations. Finding the optimal location of the temperature sensors is an iterative process with a large amount of model evaluations. Thus, the authors proposed a reduced thermo-mechanical model based on proper orthogonal decomposition (POD). The thermo-mechanical model described the thermal response of a machine tool column to localized heat sources, representing the losses of a drive and friction of the ballscrew nut. Herzog et al. considered that the heat loads had a constant position over time, resulting in a time-invariant optimal sensor placement algorithm. However, relative movements between the machine tool axes lead to moving thermal contacts and heat loads. Therefore, efficient thermo-mechanical models of machine tools need to consider the position dependency of the thermal and mechanical response.

Lang et al. <cit.> compared two methods to trace the position dependency of the thermal behavior of reduced thermal models of machine tools. In their first approach, Lang et al. divided the trajectory of the path into several segments and calculated for each of them a reduction basis by means of balanced truncation (BT). During the movement of the machine tool, the model switches to the local reduced system that is the closest to the actual position. The second reduction approach focused on a global reduction method, i.e. one projection matrix valid for all the positions of the axes. The proposed reduction method is a parametric iterative rational Krylov algorithm (IRKA), which for each parameter sample selects the optimal expansion points and tangential directions in a ℋ_2 optimal sense. Both reduction approaches showed similar accuracy compared to the full model, while the time required for reduction of the parametric IRKA was considerably smaller than for the local basis with BT. However, the amplitude of the thermal response of machine tools decays at higher excitation frequencies. Therefore, MOR approaches capturing the thermal response of the system in frequency range of interest are more efficient. Lang et al. considered the position dependency of the temperature response. However, the stiffness of the machine tool varies at different positions of the axes. Thus, it is required to couple efficiently the temperature distribution to a mechanical model that can reproduce the position-dependent mechanical response.    

Convective boundary conditions describe the interaction between the structural parts of the machine tool with the surrounding fluid media. Convective heat exchange is the main driving mechanism of the thermal response of the machine tools. It describes the heat exchange due to fluctuations of the environmental temperature <cit.> or introduction of cooling fluid into the working space <cit.>. The heat transfer coefficient (HTC) is the proportionality constant between the convective heat flux and the temperature difference between the fluid and structure. The modification of the conditions of the external fluid flow results in a variation of the HTC. Therefore, efficient thermo-mechanical models need to enable the modification of the HTC. In the context of projection-based MOR, the HTC is a model parameter that needs to be traced after reduction, i.e. parametric MOR techniques are required.

Benner et al. <cit.> provided a comprehensive review of parametric MOR methods. The parametric reduction approaches are classified into two groups, namely local or global reduction bases. A local reduction basis is constructed evaluating the system response at several samples of the parameter space. In order to switch the value of the parameter, interpolating between reduced systems is required. There are several alternatives in the literature to perform this interpolation, such as interpolating between the locally reduced system matrices <cit.>, between the local subspaces <cit.>, or between the transfer functions <cit.>. Parametric reduction with local bases provides efficient models valid for several discrete samples of the parameters. However, the convective boundary conditions of thermal models of machine tools vary continuously. Therefore, parametric MOR approaches with global reduction bases are more suitable for the models under investigation. These methods create a single set of projection bases for the values of the parameter space. 

Several reduction approaches with global bases rely on bilinearization of the systems equations. These methods express the parametric system in a bilinear form, which is a special type of non-linear systems. Phillips <cit.> used functional series expansion for the reduction of bilinear systems with rational interpolation. Continuing the work of Phillips, Bai and Skoogh <cit.> generalized the construction of the reduction basis matching a desired number of moments of the bilinear system. Breiten and Damm <cit.> proposed a reduction method for bilinear systems with Krylov subspace methods including expansion points at different values than zero. Benner and Breiten <cit.> presented a reduction approach for bilinear systems based on IRKA, named BIRKA. The authors proved that BIRKA is optimal in ℋ_2 for parameter-varying systems. Bruns and Benner <cit.> applied BIRKA to create a heat transfer model of an electrical motor. They investigated the parameters concerning the thermal contacts and convective heat transfer.

This paper proposes a parametric MOR reduction approach for thermo-mechanical models of mechatronics systems, such as machine tools. The presented MOR method considers the characteristic behavior of these models in order to build the reduction basis, enabling the possibility to trace the parameters describing the convective boundary conditions. The main contributions of this work can be summarized as:


	
  * Developing a novel MOR approach which efficiently represents the thermal behavior of the system in the frequency range of interest. Associated to the MOR method, an a-priori error estimator is proposed in order to have an upper bound of the error in the frequency range of interest.
	
  * Enabling the traceability of the parameters describing the convective boundary conditions after reduction.
	
  * Coupling efficiently the thermal response to the mechanical structural deformations. For any temperature distribution, the reduced model needs to evaluate the thermally induced structural deformations. 





§ FORMULATION OF THE PARAMETRIC THERMO-MECHANICAL MODEL


Thermo-mechanical models of machine tools describe the temperature distribution and the associated structural deformations. The temperature distribution is a continuous function T(t, z) defined at every point z of the domain Ω over time. Based on the energy conservation principle, the heat transfer equation is a partial differential equation (PDE) describing the temporal and spatial evolution of the temperature field as


    ρ c_p Ṫ(t, z) - div(λ T(t, z)) = 0


where c_p is the specific heat capacity, ρ is the material density, and λ is the thermal conductivity. The definition of the boundary and initial conditions ensure a unique solution of the PDE. The Neumann boundary conditions can be interpreted as a heat flux applied to a surface Γ_1. Neumann boundary conditions are defined as


    λ∂ T(t, z)/∂ n = q̇(t, z)


where q̇(t, z) is the heat flux applied on the boundary Γ_1. Convection is represented by Robin boundary conditions. The Robin boundary conditions are defined as 


    λ∂ T(t, z)/∂ n - h(t, z) (T(t, z) - T_ext(t, z)) = 0


where h(t, z) is the heat transfer coefficient and T_ext(t, z) is the external temperature acting on Γ_2. The Robin boundary condition represents the convective heat exchange between an external fluid media and the structure or the thermal contact between two parts. Finally, the initial condition is 


    T(z, 0) = T_0(z)


over the whole domain Ω.

The heat transfer equations describe the spatial distribution and temporal evolution of the  temperature. However, the thermally induced mechanical deviations are the main output of interest. Therefore, a coupled thermo-mechanical model is required. The thermo-mechanical coupling implies that the temperature field affects the structural deformation due to non-zero thermal expansion coefficients. However, the mechanical work resulting from the deformation of the structure does not alter the temperature field. Thus, a weakly coupled mechanical system is needed. The assumption of weak coupling holds for thermo-mechanical models with small strain rates. 

The force balance equation relates in continuum mechanics the stress tensor with the external forces as  


    -div(σ(t, z)) = f


where σ is the stress tensor at a point z in the domain Ω and f is the external force vector. Small deformations and strains are assumed, and thus a linear elastic model is considered. The stress tensor can be decomposed into an elastic part σ_e and a thermal part σ_th. Considering the constitutive equation, they can be expressed as 


    σ_e = E/1+νϵ + E ν/(1+ν)(1-2ν)trace(ϵ)I
 


    σ_th = - E/(1-2ν)α(T-T_ref)I
 

where ϵ is the strain tensor, E is the Young modulus, ν is the Poisson modulus, I is the identity tensor, α is the thermal expansion coefficient, and T_ref is the reference temperature. Under the assumption of small deformations, the strain tensor can be linearized as the gradient of the displacements u_s as


    ϵ = 1/2(∇u_s + ∇u_s ^T)
 

Numerical methods enable the solution of heat transfer and elasticity equations for general complex domains. The FE method is a well-established numerical method to solve the heat transfer and elasticity equations. Bathe <cit.> provides a comprehensive review of FE methods. Several commercial and open-source FE simulation packages are available, being extensively used in computer aided engineering (CAE). The FE-discretization interpolates a scalar field inside an element e with the values at the nodes. The interpolation functions are called antsatz functions n_e(z), defined for any point z inside the domain of the element Ω^e. Considering the interpolation of the temperature field with the ansatz function and the principle of virtual work, the PDE of Equation (<ref>) to (<ref>) is transformed into a system of ordinary differential equations (ODE) for the element e. Assembling for all the elements of the FE-mesh, the ODE for the whole system can be obtained as   


    C_thẋ(t) + K_cond(t) x(t) + K_conv(t) x(t) = q_ext


where C_th is the thermal capacity matrix, K_cond is the thermal conductivity matrix,  K_conv is the thermal convection matrix, q_ext is the thermal heat input vector, and x is the temperature state vector. 

From the different physical parameters describing Equation (<ref>), the value of the HTC are of special interest for thermo-mechanical models of machine tools. The FE-discretization of the Robin boundary condition affects both the convection matrix K_conv and the heat input vector q_ext. Equation (<ref>) defines the HTC, h(t, z), as the proportionality constant between the convective heat flux and the temperature difference between the structure and an external temperature. For the following derivations, it is assumed that the spatial distribution of the HTC does not change over time. Therefore, the HTC can be expressed separating the spatial and temporal dependency as h(t, z) = h(t)w(z). Considering this, the parametric dependency of the convection matrix and heat input vector of an element e can be expressed as 


    K^e_conv(t) = h(t)∫_Γ^e_2^n_e w(z ) n_e^T d z



    q^e_ext = h(t)∫_Γ^e_2^n_e w(z) T_ext(t, z)  dz
 

Thermal FE models usually discretize the convective boundary condition into a finite number of boundaries with the same HTC, 


    Γ_2 = ⋃_i=1^n_cΓ_i
 

where n_c is the number of independent boundaries. The HTCs for all individual boundaries, h_i(t), can be arranged in a parameter vector p∈ℝ^n_c × 1. Therefore, the system of ODE of Equation (<ref>) can be expressed in state space representation as a function of p as


    Eẋ(t) = A(p) x(t) + Bu(t)


where E∈ℝ^n × n is the left-hand side system matrix or mass matrix, A(p) ∈ℝ^n × n  is the system matrix, B∈ℝ^n × m is the input matrix, x(t) ∈ℝ^n × 1 is the state vector, and u(t)∈ℝ^m × 1 is the input vector.  Equation (<ref>) considers that the input vector also includes the corresponding value of the HTC. Considering the assumptions leading to Equation (<ref>), the parametric dependency of the state space representation of the Equation (<ref>) can be expressed in an affine form as


    Eẋ(t) = Ax(t) +  ∑_i=1^n_ch_i(t)D_i x(t) + Bu(t)

 
where D_i represents the convection matrix for each boundary Γ_i. The temperature output vector of the system y_th(t) is defined as 


    y_th(t) = C_thermx(t)


where C_therm∈ℝ^p × n is the thermal output matrix. The FE discretization of Equation (<ref>) leads to system of linear equations


    Kx_mech = K_th(x(t) - x_ref) + f_ext


where K is the stiffness matrix, K_th is the thermal coupling matrix, x_ref is the reference temperature vector, f_ext is the external force vector, and x_mech is the displacement state vector. On one hand, the thermal states x are an input of the mechanical system. On the other hand, external mechanical forces, such as preloads at the machine tool elements, affect the response of the mechanical system. Therefore, instead of defining the mechanical deviations as an output of the system of Equation (<ref>), a separate state x_mech is defined. Defining a separate state and system matrix has another practical advantage during model development. It enables the creation and validation of a mechanical model of the structure and then extend it in order to consider thermo-mechanical effects. Therefore, having an dedicated mechanical model is beneficial from a practical point of view. 

The state space representation of the mechanical system is 


    Kx_mech(t) = [  K_th B_ext ][     x(t); u_ext(t) ]
 

where K is the mechanical system matrix or stiffness matrix, B_ext is the input matrix of the external mechanical forces, and u_ext(t) is the input vector of the external mechanical forces.

In order to complete the state space representation, an output vector is defined as y_mech(t) 


    y_mech (t) = C_mechx_mech


where C_mech is the output matrix. The output of the mechanical model are the displacements at several points of interest, such as the relative deviations between tool and workpiece. 



§ FORMULATION OF THE REDUCED-ORDER MODEL


The FE-discretization of thermo-mechanical models of machine tools leads to the system equations describing the thermal transient behavior. In order to increase the computational efficiency of the models, this work proposes to create surrogate models by projection-based MOR. MOR methods construct a subspace containing the most relevant information of the high fidelity model. 

Let  𝒱 be the reduced subspace of dimension r ≪ n, where n is the dimension of the original system. Let V∈ℝ^n × r, such that 𝒱 = span (V). The reduced state of the system x̃(t) is defined as


    x(t) ≈Vx̃(t)


The different MOR approaches propose different methods to construct the reduction basis V. Once the reduction subspace is defined, the thermal system is projected into the reduction subspace substituting Equation (<ref>) into Equation (<ref>) and (<ref>). The system can be expressed as


    EVẋ̃̇(t) = AVx̃(t) +  ∑_i=1^n_ch_i(t)D_i Vx̃(t) + Bu(t)




    ỹ(t) = C_thermV_C̃_thermx̃(t)


By enforcing the Petrov-Galerkin condition, a basis W can be found such that


    W^TEV_Ẽẋ̃̇(t) = 
    W^TAV_Ãx̃(t) +
    ∑_i=1^n_ch_i(t) W^TD_i V_D̃_ix̃(t) +
    W^TB_B̃u(t)


leading to the projected system matrices Ẽ, Ã, D̃_i,  B̃, and C̃_therm and the following reduced system


    Ẽẋ̃̇(t)  = Ãx̃(t) + ∑_i=1^n_ch_i(t)D̃_i x̃(t) + B̃u(t)



    ỹ(t) = C̃_thermx̃(t)


Similar considerations can be made for the mechanical system of Equation (<ref>) and (<ref>). 


This work considers one-sided projection, i.e. W = V. As explained by Antoulas <cit.>, the main advantage of one-sided projection is that it ensures the preservation of the stability after projection.    



§ KRYLOV MODAL SUBSPACE (KMS) REDUCTION


The goal of this work is to create a reduced system that reproduces the thermal response for any value of the parameters describing the convective heat exchange between the structure and the surrounding fluid. Equation (<ref>) presents an affine representation of the parametric dependency, separating the system matrix into two terms. The first step considers the thermal part of the system that does not depend on the parameter. A reduction basis is created such that the reduced system reproduces the response of the system where the values of the parameters are set to a fixed value. This section introduces the MOR approach and an error bound for the error between the original and the reduced-order model. In a second step, the projection basis is extended in order to enable the modification of the HTC after reduction, which is the focus of the next section. 

The projection basis needs to capture the most relevant part of the response of the system. The thermal models of mechatronic systems, such as machine tools, have a certain characteristic behavior compared to other general first order systems. Firstly, the amplitude of the thermal response of the system decays at higher excitation frequencies. Therefore, the reduced system only needs to approximate the thermal response of the system in a low frequency range. Secondly, the steady state response is relevant in order to describe the behavior of the system. Thus, the reduced system needs to match the steady state temperature field of the original system for the different input and output combinations. 

This work proposes a reduction method that considers the characteristic behavior of thermal models of machine tools to construct the projection basis. On one hand, moment matching Krylov subspace methods approximate the response of the system around an expansion point s_e. The expansion point is placed at a low frequency close to zero in order to match the steady state response of the system. On the other hand, the eigenvectors of the system up to a certain frequency are included in the projection basis. Therefore, the thermal reduced system reproduces the thermal behavior of the system in the frequency range of interest. Spescha <cit.> introduced the concept of KMS reduction with application in structural dynamics. This paper extends the KMS method to first order systems and presents a theoretical bound of the reduction error.

Let 𝒱_k ⊂ℝ^n be the Krylov subspace with one expansion point s_e, such that 


    𝒱_k = span((s_eE-A)^-1B)


An orthonormal basis V_k of 𝒱_k  can be constructed, such that 𝒱_k = span (V_k). The original system of Equation (<ref>) can be projected into the subspace 𝒱_k. The reduced system matches the response of the system around the expansion point s_e. If an expansion point s_e is chosen close to 0, the reduced system matches the steady state response. 

Let 𝒱_μ⊂ℝ^n be the truncated modal subspace. The subspace 𝒱_μ is the span of the first μ eigenvectors of the system of Equation (<ref>), defined as 


    𝒱_μ = span([ ϕ_1 ϕ_2   … ϕ_μ ])


where ϕ_i is an eigenvector associated to the eigenvalue α_i such that (s_eE - A)ϕ_i = α_i ϕ_i. The system matrices A and E, defined in Equation (<ref>), (<ref>), (<ref>), and (<ref>), are negative semi-definite, resulting in all real non-positive eigenvalues. Therefore, the eigenvectors form an orthonormal basis V_μ of 𝒱_μ, such that 𝒱_μ = span (V_μ). The original system of Equation (<ref>) can be projected into the subspace 𝒱_μ. The reduced system approximates the dynamic thermal behavior in a certain frequency range. 

The KMS reduction projects the system of Equation (<ref>) by means of a projection matrix V. The basis V spans linear subspace 𝒱⊂ℝ^n, i.e. 𝒱 = span (V), such that


    𝒱 = 𝒱_μ + 𝒱_k = {x∈ℝ^n  /  ∃v_1 ∈𝒱_μ v_2 ∈𝒱_k x = v_1 +v_2 }


The subspace 𝒱 captures the information about both the steady state response and the thermal transient behavior of the system. Therefore, this projection basis satisfies the requirement for accurate approximation of the thermal behavior of mechatronic systems.



 §.§ A-priori error estimator


The reduced system needs to represent accurately the thermal response of the system in the frequency range of interest, i.e. for all ω∈ [0 , ω_max]. The upper bound of the frequency range of interest, ω_max, determines how many eigenvectors, μ, are to be included in the KMS reduction basis. Therefore, this work proposes an a-priori error estimator for the KMS method that relates the maximum eigenfrequency, ω_μ = |α_μ|, of the KMS basis with the frequency range of interest. 

The error estimator presented in this section considers that the input and output matrix of Equation (<ref>) and (<ref>) are the same, i.e. B = C^T_therm, leading to one-sided projection. For the derivation of the error estimator of the KMS method, the Krylov subspace defined in Equation (<ref>) has an expansion point s_e ∈ℝ close to zero and a single iteration. Before introducing the error estimator, some definitions and preliminary results are introduced. 

Firstly, a suitable error definition is required. Different error values are proposed in the literature, as summarized by Benner et al. <cit.>. Let E(jω) be the absolute reduction error frequency response function (FRF) defined for each frequency ω as


    E(jω) = H(jω) - H̃(jω)
  

where H(jω) and H̃(jω) are the FRF of the original and reduced system respectively. E(jω) is a matrix of dimension p (number of outputs) by m (number of inputs). Let e_ij(jω) be the relative error for the ith input and jth output combination as 


    e_ij(jω) =  h_ij(jω) - h̃_ij(jω)  / h_ij(jω)


where h_ij(jω) is the element in ith row and jth column of H(jω), and h̃_ij(jω) is the element in ith row and jth column of H̃(jω). 

Secondly, some remarks about the subspaces associated to the KMS 𝒱 are required. The subspace 𝒱_ν⊂ℝ^n can be defined as the subspace of the remaining modes not included in 𝒱_μ, i.e.

    𝒱_ν = span([ ϕ_μ+1 ϕ_μ+2     …   ϕ_n ])


Due to the properties of the system matrix, the subspace 𝒱_μ is the orthogonal complement of 𝒱_ν, such that ℝ^n = 𝒱_μ⊕𝒱_ν. Let Φ be a matrix whose columns are the eigenvectors ϕ_i of the system of Equation (<ref>), normalized to the capacity matrix such that Φ^TEΦ = I. The system of Equation (<ref>) and (<ref>) can be expressed in modal coordinates x = Φx_m as


    Iẋ_m(t) = Ωx_m(t) + Φ^TBu(t) =  Ωx_m(t) + B_mu(t)


    y_th(t) = C_thermΦx_m(t) = C_mx_m(t)


where Ω = diag(α_1, …, α_n) is a diagonal matrix with the n eigenvalues α_k of the system. The system matrix can be expressed as block matrix, such that


    [ I_μ   0;   0 I_ν ][ ẋ_μ; ẋ_ν ] = [ Ω_μ   0;   0 Ω_ν ][ x_μ; x_ν ]+[ B_μ; B_ν ]u(t)


    y_th = [ C_μ C_ν ][ x_μ; x_ν ]



Let 𝒱_ν k⊂ℝ^n be the Krylov subspace with one moment around the expansion point s_e ∈ℝ of the original system projected into the subspace 𝒱_ν, such that 


    𝒱_ν k = span((s_eI_ν-Ω_ν)^-1B_ν)


The definition of 𝒱_ν k leads to a preliminary result expressed in Fact <ref>. 


	
	The KMS 𝒱 is the direct sum of 𝒱_μ and 𝒱_ν k.	
	
    𝒱 = 𝒱_μ⊕𝒱_ν k

	



	
	The Krylov subspace of the system in modal coordinates Equation (<ref>) is 
	
	
    𝒱_k = span((s_eI-Ω)^-1B_m)

	
	As s_eI+Ω is a diagonal matrix, the Krylov subspace of the original system can be expressed as 
	
	
    𝒱_k = 𝒱_μ k + 𝒱_ν k = span((s_eI_μ-Ω_μ)^-1B_μ) + span((s_eI_ν-Ω_ν)^-1B_ν)

	
	The 𝒱_μ k is a subspace of 𝒱_μ, i.e. 𝒱_μ k⊂𝒱_μ. Therefore, 𝒱 = 𝒱_k +𝒱_μ = 𝒱_μ k + 𝒱_ν k + 𝒱_μ = 𝒱_ν k + 𝒱_μ. 
	
	In addition, 𝒱_ν k is also a subspace of 𝒱_ν, i.e. 𝒱_ν k⊂𝒱_ν. Since 𝒱_ν∩𝒱_μ = 0, it follows that 𝒱_ν k∩𝒱_μ = 0. This implies that 𝒱_ν k is the orthogonal complement of 𝒱_μ, i.e. 𝒱 =  𝒱_ν k⊕𝒱_μ
 
	



Fact <ref> relates the subspaces 𝒱, 𝒱_μ, and 𝒱_ν k, stating that any vector in 𝒱 can be decomposed uniquely into two terms, one in 𝒱_μ and another one in 𝒱_ν k. This property is useful to separate the error of the reduced system, as shown in Fact <ref>.   


	
	The error e_ij(j ω) of Equation (<ref>) can be separated as the product of two terms	
	
    e_ij(j ω) = -e_μ_ij(jω )e_ν k_ij(jω )

	
	being e_μ_ij(jω ) and e_ν k_ij(jω ) defined as 
	
    e_μ_ij(j ω) = h_ij(jω ) - h̃_μ_ij(jω )/h_ij(j ω )

	
    e_ν k_ij(j ω ) = h̃_ν k_ij(j ω ) - h̃_ν_ij(j ω )/h̃_ν_ij(j ω )

	
	where h̃_μ_ij(jω ) is the element in ith row and jth column of the FRF of the system projected into 𝒱_μ, h̃_ν_ij(jω ) is the element in ith row and jth column of the FRF of the system projected into 𝒱_ν, and h̃_ν k_ij(jω ) is the element in ith row and jth column of the FRF of the system projected into 𝒱_ν k.
	
	




	The Fact <ref> enables to express the FRF of the reduced system as h̃_ij(j ω)  = h̃_μ_ij(j ω) + h̃_ν k_ij(j ω). Additionally, the transfer function of the original system can be expressed as h_ij(j ω)  = h̃_μ_ij(j ω) + h̃_ν_ij(j ω). Substituting in the error definition 
	
    e_ij(j ω) = h̃_μ_ij(j ω) + h̃_ν_ij(j ω) - h̃_μ_ij(j ω) - h̃_ν k_ij(j ω)/h_ij(j ω)

	
	Multiplying the previous expression by h̃_ν_ij(j ω)/h̃_ν_ij(j ω), the following is obtained
	
	
    e_ij(j ω) = h̃_ν_ij(j ω)  - h̃_ν k_ij(j ω)/h_ij(j ω)h̃_ν_ij(j ω)/h̃_ν_ij(j ω)

	
	Substituting h̃_ν_ij(j ω) in the numerator by h_ij(j ω)  - h̃_μ_ij(j ω) and reorganizing the terms, the error separation of Equation (<ref>) is obtained as
	
	
    e_ij(j ω) = -h̃_ν k_ij(j ω) - h̃_ν_ij(j ω)/h̃_ν_ij(j ω)·h_ij(j ω)  - h̃_μ_ij(j ω)/h_ij(j ω)

	

	


The result of Fact <ref> enables the separation of the error in two terms. The goal is to find an bound for the error e_ij(j ω). The first step for an error bound of the KMS method is to determine an upper bound for the term e_ν k_ij(j ω). The following theorem shows a theoretical error bound for all frequencies of the term e_ν k_ij(j ω), provided that the input and the output are the same, namely i = j. 


	The magnitude of the error e_ν k_ii(j ω) defined in Equation (<ref>) is bounded by e_est(jω) defined as
	
	
    | e_est(jω) | = ω^2 + s_e^2/ω ^2 + w_m^2 > | e_rk_ii(jω) |

	
	for all ω∈ [0 , ∞] such that ω∈ℝ given that ω_m < ω_mo+1 and that the input and the output are the same, i.e. i = j. 
	
	



	The FRF of the system projected into 𝒱_ν, is
	
	
    h̃_ν_ii(jω) = ∑_k=1^νb_k^2/(jω + ω_k)

	
	where ν = n - μ is the dimension of the subspace 𝒱_ν, ω_k are the absolute value of the eigenvalues of the system, and b_k corresponds to the kth element of the ith column of B_ν. The subspace 𝒱_ν k can be defined according to Equation (<ref>) as range(v) where	
	
	
    v^T = [ b_1/ω_1+s_e b_2/ω_2+s_e           … b_k/ω_k+s_e           … b_r/ω_r+s_e ]

	
	Projecting the system of Equation (<ref>) into 𝒱_ν k, the following system is obtained
	
	
    ẽẋ̃̇ + ãx̃ =  b̃ u(t)

	
	where ẽ, ã, and b̃ take the following values 
	
	
    ã = ∑_k=1^νb_k^2/(ω_k+s_e)^2ω_k

	
	
    ẽ = ∑_k=1^νb_k^2/(ω_k+s_e)^2

	
	
    b̃ = ∑_k=1^νb_k^2/ω_k+s_e = c̃

	
	The transfer function of the reduced system h̃_ν k is
	
	
    h̃_ν k_ii(j ω) = b̃^2/j ωẽ +ã =(∑_k=1^νb_k^2/ω_k+s_e)^2/j ω∑_k=1^νb_k^2/(ω_k+s_e)^2+ ∑_k=1^νb_k^2/(ω_k+s_e)^2ω_k

	
	Firstly, some properties of the error e_ν k_ii(j ω) need to be discussed. According to Equation (<ref>), the poles of the transfer function of the real error, e_ν k_ii(s) with s ∈ℂ, are the poles of h̃_ν k_ij(s), the poles h̃_ν_ii(s), and the zeros of h̃_ν_ii(s). Given that the system matrices are positive semi-definite, all the poles of h̃_ν_ii(s) are real. Additionally, the poles of h̃_ν_ii(s) are real, as it can be seen from Equation (<ref>). Furthermore, the zeros of the transfer function h̃_ν_ii(s) are real numbers. This fact can be shown by contradiction. Let s = α± jβ be zeros of the transfer function h̃_ν_ii(s), such that β > 0. The zero of the transfer function needs to satisfy that 
	
	
    h̃_ν_ii(α + jβ) = ∑_k=1^νb_k^2/(α + jβ + ω_k) = ∑_k=1^νb_k^2 /((α +ω_k)^2 + β^2) (α +  ω_k - jβ )  = 0

	
	Given that β > 0 and b_k^2 /((α +ω_k)^2 + β^2)≥ 0, the transfer function is only zero for β = 0, reaching a contradiction. Therefore, the zeros of transfer function h̃_ν_ii(s) are real. Thus, the poles of e_ν k_ii(j ω) are all real. This leads to a smooth FRF response, without any resonance frequency. This property is relevant for the derivation of the error bound.
	
	In order to proof that the proposed estimator of Equation (<ref>) bounds the error for all frequencies, several counterexamples are created, which are illustrated in Figure <ref>. Figure <ref> shows the case where the magnitude of the error, e_ν k_ii(jω), is higher than the magnitude of the estimator, e_est(jω), at high frequencies. Figure <ref> depicts the case where the poles of the error at lower frequencies than the poles of the estimator. The third counterexample in Figure <ref> illustrates the case where the slope of the error is lower than the slope of the estimator. The slope of the error is related with the number of zeros at low frequency. The error estimator has two zeros at the expansion point s_e. Thus, it needs to be shown that the actual error at least two zeros at the expansion point, which is placed close to zero. Figure <ref> shows graphically that if the other three counterexamples are not true, the error estimator is an upper bound of the FRF of the of the real error for the whole frequency range. 
	
		
	
	
	Therefore, it needs to be proven that the counterexamples of Figure <ref> are not possible following the next steps:
	
	
		
  * lim_ω→∞| e_est(j ω) |> lim_ω→∞| e_ν k_ii(j ω) |
		
  * All the poles of the actual error e_ν k_ii(j ω) are at a higher frequency than the ones of the error estimator e_est(j ω)
		
  * e_ν k_ii(j ω) has at least two zeros at ω close to zero
	
	
	The first condition refers to the magnitude of the FRF at infinity, which is associated to the Markov parameters. From the definition of the error estimator lim_ω→∞| e_est(jω) | = 1. The magnitude of e_ν k_ii(j ω) at infinity can be expressed as 
	
	
    lim_ω→∞| e_ν k_ii(jω) | = lim_ω→∞|h̃_ν k_ii(j ω)/h̃_ν_ii(j ω)  -1| = |lim_ω→∞h̃_ν k_ii(j ω)/h̃_ν_ii(j ω)  -1| < 1

	
	
	In order to evaluate the limit of the error, the limit lim_ω→∞h̃_ν k_ii(j ω)/h̃_ν_ii(j ω) needs to be evaluated. 
	
	
    lim_ω→∞h̃_ν k_ii(j ω)/h̃_ν_ii(j ω) = 1/j ωb̃^2/ẽ/1/j ω∑_k=1^νb_k^2 = b̃^2/ẽ/∑_k=1^νb_k^2

	
	where b̃^2/ẽ can be expressed according to Equation (<ref>) as
	
	
    b̃^2/ẽ = (∑_k=1^νb_k^2/ω_k+s_e)^2/∑_k=1^νb_k^2/(ω_k+s_e)^2

	
	Equation (<ref>) shows that the limit lim_ω→∞h̃_ν k_ii(j ω)/h̃_ν_ii(j ω) is a positive real number. Therefore, the condition on the magnitude of the error of Equation (<ref>) can be expressed as 
	
	
    0 < lim_ω→∞h̃_ν k_ii(j ω)/h̃_ν_ii(j ω) <2

	
	In fact, an upper bound for the term b̃^2/ẽ can be found as follows
	
	
    b̃^2/ẽ < (∑_k=1^νb_k^2/ω_k+s_e)^2/1/ω_1+s_e∑_k=1^νb_k^2/ω_k+s_e =  ∑_k=1^νω_1+s_e/ω_k+s_e b_k^2 < ∑_k=1^νb_k^2

	
	
	considering that ω_1 ≤ω_k for all k ∈ [1, ν]. The upper bound for the term b̃^2/ẽ implies that the condition of Equation (<ref>) is satisfied. Therefore, the magnitude of the actual error is bounded by 1, as stated in Equation (<ref>). This completes the first step of the proof, showing that the magnitude of the error estimator is an upper bound of the magnitude of the actual error, e_ν k_ii(j ω), at infinity.
	
	
	The second step in this proof is ensuring that the poles of the actual error are at a higher frequency than the ones of the error estimator. Let e_ν k_ii(s) the transfer function of the actual error, where s ∈ℂ. The error estimator, e_est(s), defined in Equation (<ref>) has two poles at ω_m. The poles of the error e_ν k_ii(s) are the poles of h̃_ν k_ii(s), the poles h̃_ν_ii(s), and the zeros of h̃_ν_ii(s), which need to be at a lower frequency than ω_m. 
	
	The poles of h̃_ν_ii(s) are ω_1 …ω_ν. By definition, these poles are at a higher frequency than ω_m. 
	
	The zeros of h̃_ν_ij(s) are real numbers, as discussed in the beginning of this proof. Furthermore, it needs to be shown that these zeros are at higher frequencies than ω_m, which is proven by contradiction. Assume that ω_0 is a zero of the transfer function, i.e. h̃_ν_ii(s=-ω_0)=0, such that 0<ω_0 <ω_1. Substituting this into Equation (<ref>), the transfer function is 
	
	
    h̃_ν_ii(s=-ω_0) = ∑_k=1^νb_k^2/(-ω_0 + ω_k)

	
	However, given that 0≤ω_1 ≤ω_ν the all the terms of the summation are strictly positive. Thus,  h̃_ν_ii(s=-ω_0)≠ 0, reaching a contradiction. Therefore, it is shown that the zeros of the transfer function of the reduced system h̃_ν_ii(s) are at a higher frequency than ω_m. 
	
	The pole h̃_ν k_ii(s), ω_ν k_ii, can be written as
	
    ω_ν k_ii = ã/ẽ = ∑_k=1^νb_k^2/(ω_k+s_e)^2ω_k/∑_k=1^νb_k^2/(ω_k+s_e)^2

	
	which can be understood as a weighted summation, where the weight factors are b_k^2/(ω_k+s_e)^2. This shows that the poles of the transfer function satisfy ω_1 ≤ω_0 ≤ω_r. The location of ω_0 depends on the static gains,  b_k. The more controllable and observable a mode ω_k is, the closer ω_0 is to ω_k. Thus, the pole of the transfer function of the reduced system h̃_ν k_ii(s) is at a higher frequency than ω_m. This concludes the second step, ensuring that the poles of the actual error e_ν k_ii(s) are at a higher frequency than the ones of the error estimator.

	
	The third step needs to ensure that the transfer function of the actual error, e_ν k_ii(s), has at least two zeros close to zero. The zeros of the actual error are related with the zeros h̃_ν k_ii(s) - h̃_ν_ii(s). Due to the properties of moment matching Krylov reduction, the reduced-order system matches at least two moments of the transfer function at the expansion point, s_e, as explained by Antoulas <cit.>. By choosing the expansion point sufficiently low, the h̃_ν k_ii(s) - h̃_ν_ii(s) has two zeros close to zero. Therefore, the transfer function e_ν k_ii(s) has also two zeros close to zero. 
	
	This proof shows that the error e_ν k_ii( j ω) has a bounded magnitude at infinity, the poles of the error are a higher frequency than the poles of the estimator, and the error has at least two zeros close to zero frequency. The combination of these three conditions proves that the magnitude of the error e_ν k_ii( j ω) is bounded by the estimator | e_est(jω) |.
	
	



The result of Theorem <ref> proposes an error bound of the term e_ν k_ii(j ω) of Equation (<ref>). In order to estimate the reduction error, e_ii(j ω), an upper bound of the term e_μ_ii(j ω) of Equation (<ref>) is also required. 


	The magnitude of the error e_μ_ii(j ω) defined in Equation (<ref>) is bounded by 1, i.e.	
	
    | e_μ_ii(jω) |≤ 1

	
	for all ω∈ [0 , ∞] such that ω∈ℝ. 
	
	




	The FRF of the error e_μ_ii(j ω) defined in Equation (<ref>) of Fact <ref> can be expressed as 
	
	
    e_μ_ii(j ω) = h̃_ν_ii(jω )/h_ii(j ω )

	
	considering that h̃_ν_ii(jω ) = h_ii(jω ) -h̃_μ_ii(jω ), as shown in Fact <ref>. In order to prove this fact, the two extreme conditions are considered. Firstly, the case where all modes in 𝒱_ν are not observable and not controllable is analyzed. In this case, the magnitude of h̃_ν_ii is zero. Thus, the magnitude of the error | e_μ_ii(j ω) | is zero for all frequencies, as all the relevant modes describing the response of the system are contained in 𝒱_μ. Secondly, it is considered that all modes in 𝒱_μ are not observable and not controllable. In this case, the original system and the system projected into 𝒱_ν  are the same, i.e. h̃_ν_ii = h_ii(j ω ). Thus, the magnitude of the error | e_μ_ii(j ω) | is 1, as none of the relevant modes describing the response of the system are included in 𝒱_μ. These two cases represent the two extreme conditions. For any other case, some of the modes 𝒱_μ and 𝒱_ν are observable and controllable, leading to a magnitude of the error between 0 and 1. Therefore, | e_μ_ii(jω) |≤ 1 for all frequencies. 


The results of Theorem <ref> and Fact <ref> state that an error bound of e_ii(j ω) for the KMS reduction error is 


    | e_ii(j ω) | = | e_μ_ii(jω ) || e_ν k_ii(jω ) | < | e_ν k_ii(jω) | <  ω^2 + s_e^2/ω ^2 + w_m^2



Theorem <ref> shows a theoretical error bound for the KMS reduction, provided that the same inputs and outputs are considered. This error bound estimates a-priori the error and enables to choose how many modes need to be included in the KMS projection basis. Provided a frequency range of interest [0, ω_max], the error estimator determines that all modes under ω_m need to be included in the reduction basis so that the magnitude of the error of the reduced-order system, | e_ii(j ω) |, does not exceed a certain value ϵ. Section <ref> illustrates the validity of the error estimator with numerical examples for different combinations of the inputs and outputs.   



§ PARAMETRIC KMS REDUCTION


The KMS method creates a reduction basis that captures the most relevant information of the reduction system in the frequency range of interest. This section presents a method to extend the KMS basis to include the possibility to modify the values of the HTC in the reduced system. The system of Equation (<ref>) shows the parametric dependency as an affine representation. This affine formulation of the system can be interpreted as bilinear system. The term ∑_i=1^n_ch_iD_i x(t) can be understood as an input of the system that depends linearly on the state. Bilinear systems are a special class of non-linear systems that are linear in the input and state, but not jointly linear in the input and state.

Before presenting the reduction approach, several considerations are required on the convection matrices D_i of Equation (<ref>). Equation (<ref>) defines the convection matrix of an element e, which is integrated numerically. The evaluation of the integral at the Gauss points leads to off-diagonal terms in K^e_conv. This is called consistent convection matrix. However, it is customary to diagonalize this matrix in order to avoid numerical oscillations close to the convective boundary condition, as pointed out by Bruns <cit.>. Therefore, a diagonal convection matrix D_i is considered, with zeros for all the degrees of freedom (dof) not affected by the convective boundary condition i. The diagonal convection matrix D_i can be understood as a linear map 𝒟_i, whose range has a dimension equal to the number of nodes in the convective boundary,  n_d_i.

Another important consideration for the development of the reduction method is determining the dimension of the subspace spanned by (s_eE-A)^-1D_i. The matrix (s_eE-A)^-1 can be understood as a linear transformation, 𝒜. Given that (s_eE-A)^-1 is invertible, 𝒜 is a bijective linear transformation, i.e. dim(range((s_eE-A)^-1)) = n. Thus,


    dim(range(𝒜∘𝒟_i)) = n_d_i
  

The intersection of the subspaces spanned by the several convection matrices also needs to be discussed. In thermo-mechanical models of mechatronic systems, the different convective boundary conditions affect different surfaces. Thus, it can be stated that Γ_i ∩Γ_j = 0  ∀ i≠ j. Therefore


    range(𝒟_i) ∩range(𝒟_j) = 0  ∀ i≠ j


This can be extended to the range of 𝒜∘𝒟_i as


    range(𝒜∘𝒟_i) ∩range(𝒜∘𝒟_j) = 0  ∀ i≠ j


Let V_KMS be the KMS reduction basis of the system without parametric convective boundary conditions. The basis V_KMS  spans the linear subspace defined in Equation (<ref>). All the states in the KMS are a linear combination of the columns of V_KMS, i.e. x∈span(V_KMS).


In order to consider the consider the state-dependent inputs, D_i x, an iterative approach is applied. The states multiplied by the convection matrix, i.e.  D_iV_KMS, provide a new set of inputs for the reduction. The matrix D_iV_KMS can be added as further inputs to the Krylov subspace of Equation (<ref>), providing a new set of vectors V^1_i. This set of new basis vectors multiplied by the convection matrix, i.e. D_iV^1_i, provide new inputs. They can be added to the Krylov subspace delivering a new set of basis vectors V^2_i. This iterative process can be expressed as


    span(V_KMS) = span((s_eE-A)^-1B) + span(V_μ) 
    V^1_i = (s_eE-A)^-1D_iV_KMS
    …
    V^k_i = (s_eE-A)^-1D_iV^k-1_i
    …
    V^n_me_i = (s_eE-A)^-1D_iV^n_me-1_i


where n_me is the number of iterations. The sum of the subspaces spanned by the matrices V_KMS, V^1_i, …, V^k_i, …, V^n_me_i is the reduced subspace 𝒱_i 


    𝒱_i = span{V_KMS}  + span{ (s_eE-A)^-1D_iV_KMS} + …
    + span{ ( (s_eE-A)^-1D_i  )^n_me-1V_KMS}
 

Comparing the definition of 𝒱_i with the definition of the Krylov subspace, it can be stated that 


    𝒱_i = 𝒦_n_me{(s_eE-A)^-1D_i, V_KMS}
  

For the creation of the subspace 𝒱_i, n_me < n_d_i terms are selected. This process of adding new vectors to the reduction basis V_KMS can be understood as a Krylov reduction of the subspace spanned by (s_eE-A)^-1D_i. The reduced subspace can be calculated for all the n_c convective boundary  conditions. Given 𝒱_i and 𝒱_j for two different convective boundary conditions, it can be stated that 


    𝒱_i∩𝒱_j = span(V_KMS)


which comes from the properties of the convection matrices shown in Equation (<ref>). The reduction subspace 𝒱_p is the addition of subspaces 𝒱_i as 


    𝒱_p =  ∑_i=1^n_c𝒱_i = span(V)


where V is an orthonormal basis of the subspace 𝒱_p. 

The KMS reduction considers certain modes of the system, V_μ, in order to form the projection basis. The KMS method selects modes up to a certain frequency. The error estimator ensures that the error to the original system remains below a given tolerance for the frequency range of interest. However, the system matrix changes with the introduction of the convective boundary conditions, as shown in Equation (<ref>). Thus, the eigenvalues and eigenvectors of the system change. It is required to determine if more modes need to be considered in order to ensure that the error remains bounded below a certain ϵ in the frequency range of interest. 

In order to investigate this, the Weyl's Inequality Theorem is presented as a preliminary result. Theorem <ref> bounds the eigenvalues of the sum of two matrices M + N given the eigenvalues of M and N. 

Weyl's Inequality
	Let M, N be symmetric matrices of dimension ℝ^n and let S = M + N. Let α_1 ≥α_2 ≥…α_n, β_1 ≥β_2 ≥…β_n, and γ_1 ≥γ_2 ≥…γ_n be the eigenvalues of the matrices M, N, and S respectively. Then,
	
	
    γ_j  ≤α_i + β_j-i+1    ∀ i ≤ j

	
 


	See Bhatia <cit.> Chapter 3. 



The result of Theorem <ref> can be used to analyze the properties of the eigenvalues of the system of Equation (<ref>) with increasing values of the HTC h_i.  


	Let the system matrix be A^l_d = A +  ∑_i=1^n_ch^l_iD_i, where h^l_i≥0 is a sample of the ith parameters defining the n_c convective boundary conditions. Let 0 ≥α^l_1 ≥α^l_2 ≥…α^l_n be the eigenvalues of A^l_d. Let the system matrix be A^m_d = A +  ∑_i=1^n_ch^l_iD_i, where h^m_i is another sample of the i parameters such that h^m_i ≥ h^l_i≥ 0. Let 0 ≥α^m_1 ≥α^m_2 ≥…α^m_n be the eigenvalues of A^m_d. Then, α^m_j ≤α^l_j for j =1 … n. 
	



	The system matrices A and D_i are negative semi-definite, i.e. all the eigenvalues are real and non-positive. The result of Theorem <ref> can be particularized for i=j as
	
	
    γ_j  ≤α_j + β_1
 
	
	Additionally, considering a negative semi-definite matrix N it can be stated that 
	
	
    γ_j  ≤α_j + β_1≤α_j
 
	
	as β_1≤ 0. Let δ h_i be the difference between h^m_i and h^l_i for i = 1 … n_c, i.e. Δ h_i = h^m_i - h^l_i. From the statement of this theorem, Δ h_i ≥ 0. The system matrices can be expressed as A^m_d = A^l_d + ∑_i=1^n_cΔ h_iD_i. Applying the particularization of the Theorem <ref> to negative semi-definite matrices of Equation (<ref>), the following result is obtained
	
	
    α^m_j ≤α^l_j

	



The result of Theorem <ref> states that the higher the HTC, the more negative the eigenvalues of the system are. The physical interpretation of this theorem is that the higher the convective heat exchange with the surrounding, the faster the time constants of the system are.  On the other side, a system with low HTC evacuates less efficiently the heat and thus has slower time constants. 


Theorem <ref> can be particularized for the case that the initial HTC are equal to zero, i.e. h^l_i = 0 for all i = 1 … n_c. This is the case of the first iteration of the bilinearization, as shown in Equation (<ref>). The eigenvalues of the system without convective boundary conditions are less negative. Thus, the eigenfrequcies of the selected eigenmodes in V_μ increase in magnitude after including the convection matrices. Therefore, the reduced system with proposed bilinearization remains valid in the frequency range of interest.  

The proposed parametric MOR approach extends the KMS projection basis to include the parametric dependency of the HTC. Therefore, the dimension of the reduced system increases with respect to the non-parametric reduced model. The dimension of the parametric model, r_p, depends on the dimension of the number of convective boundary conditions, n_c, the number of iteration, n_me, and the dimension of the non-parametric reduced system, r,  as


    r_p = r(1+n_men_c)


 



§ THERMO-MECHANICAL COUPLING


The KMS reduction constructs a projection basis that reduces the thermal system of Equation (<ref>). The thermal model determines the transient behavior of the system. However, the main output of interest of the thermo-mechanical model is the mechanical deformation of the structure. Therefore, a reduction basis for the system of Equation (<ref>) and (<ref>) needs to be defined. 

The mechanical model only considers the quasi-static response of the structure. Therefore, the reduced subspace for the mechanical system can be constructed using a Krylov subspace with one moment and one expansion point at a low frequency. The reduced system matches the static gain of the mechanical system. In order to construct the projection basis, the number of inputs of Equation (<ref>) needs to be considered. The number of mechanical inputs is equal to the number of possible independent temperature distributions and the number of external loads. In principle, the number of independent temperatures is n, i.e. the dimension of the original thermal system. However, the reduced thermal model determines the temperature distribution. Thus, the possible temperature distributions are limited to a linear combination of the vectors of the reduction basis V. All the possible linearly independent thermal body forces F_th are


    F_th = K_thV


Therefore, the state space representation of Equation (<ref>) can be expressed in terms of the reduced temperature state x̃ instead of the full state x as


    A_mechx_mech = [ K_thV F_ext ][    x̃; u_ext ]
 

where the number of mechanical inputs is the order of the thermal reduced system and the external forces. The reduction basis V_mech for the mechanical system can be calculated using the Krylov subspace method as


    span(V_mech) = 𝒦_r {(s_eI-A_mech)^-1, (s_eI-A_mech)^-1[ K_thV F_ext ]}
 

where the expansion point s_e is chosen close to zero. The Krylov subspace matches exactly static response of the system of Equation (<ref>) for the different temperature distributions considered in the thermal system projected by V. 



§ NUMERICAL RESULTS



This section illustrates the MOR methods developed in this paper with numerical examples. Figure <ref> shows the geometry and FE-discretization of a 3-axis milling machine. The thermo-mechanical model represents the thermal response of the machine tool under internal heat sources, such as the heat dissipated in the linear drives, and external effects, such as fluctuations of the environmental temperature.




Figure <ref> shows the geometry and FE-mesh of the table of the machine tool, which can move in Y-direction. The simple geometry under consideration leads to an original system with a small number of dofs, i.e. 4157 thermal dofs. The small size of the original system facilitates the evaluation of its thermal response in order to compare it with the reduced system. The convective boundary conditions of the machine tool table are separated into two regions, as shown in Figure <ref>. Each of the convective boundary conditions has an independent value of the HTC, namely HTC_top and HTC_bot. Additionally, the bottom of the table is exposed to a heat flux, representing the heat losses of a linear drive.




The a priori error estimator of Equation (<ref>) facilitates the selection of the parameters for the KMS reductions. Given a maximum error between the reduced and original system, ϵ, and a frequency range of interest, ω∈ [0, ω_max], the error estimator provides the value of ω_m. The reduction basis includes all eigenvectors with eigenfrequencies below the frequency ω_m, i.e. ω_μ < ω_m. For the FE thermal model under consideration, the following parameters are chosen: 


	
  * ϵ =  0.05
	
  * ω_max = 0.01 rad/s


The maximum frequency of interest, ω_max, determines the frequency range containing the most relevant part of the response to the considered thermal inputs. For the thermal model of Figure <ref>, the value of ω_m is 0.04367 rad/s. This results in including 97 modes in the KMS projection basis V_KMS. The KMS reduction requires additionally an expansion point at a low frequency in order to match the steady state response of the system. For this thermal model, the reductions basis V_KMS considers an expansion point s_e = 10^-8 rad/s. The bilinear reduction extends the basis V_KMS with the information about the convective boundary conditions. For the investigated thermal model, the number of parametric convective boundary conditions, n_c, is 2. The number of iterations for the bilinear reduction (m_d) is set to 2. The resulting reduced model has 496 dof, which is considerably smaller than the dimension of the original system.  

In order to evaluate the performance of the reduction method, the transfer functions of the thermal response of the original and reduced systems are compared. The input of the transfer function is the fluctuation of the environmental temperature. The output of the transfer function is the temperature measured at the drive of the Y-axis. The transfer functions are evaluated for different values of the parameters describing the convective boundary conditions. The error between the reduced and the original system is calculated between 10^-5 and 1 rad/s according to Equation (<ref>) for the considered system input and output. Figure <ref> depicts the relative error between original and the system reduced by means of the parametric KMS method. The error is negligibly small at low frequencies, as the chosen expansion point matches the steady state response. At high frequencies, the error increases reaching values close to 1. The reduced system succeeds in reproducing the thermal response of the system in the frequency range of interest, i.e.  up to 0.01 rad/s. Additionally, Figure <ref> shows that the error estimator of Equation (<ref>) is an upper bound of the relative error for all frequencies.



Figure <ref> shows the reduction errors for one single combination of the input, environmental temperature, and output, mean temperature at the drive. Figure <ref> extends the evaluation of the relative error of the reduced and original system to other input and output combinations, fixing the values of HTC_top and HTC_bot to 4 and 8 W/m^2K respectively. For the considered numerical example, the relative errors for different inputs and output combinations show more variability than the errors for Figure <ref>. The proposed error estimator remains an upper bound of the relative error for the transfer functions of the reduction error.    



In order to study the thermal response of the reduced system, it is interesting to compare the eigenvalues of the reduced and the original systems. Let ω_i be the ith eigenvalue of the system of Equation (<ref>) and let ω̃_i be the ith eigenvalue of the system projected to the subspace spanned by the projection basis V of Equation (<ref>). The relative error for the ith eigenvalue can be defined as 


    e_i = |ω̃_i-ω_i /ω_i|


Figure <ref> shows the relative error of the first 90 eigenvalues between the original and reduced system. The eigenvalues of the systems lie in the frequency range of interest, ranging from 7 · 10^-5 to 0.042 rad/s. The relative errors are evaluated for different combinations of the HTC, considering the same values of the HTC parameters as in Figure <ref>. The relative error of the eigenvalues increases at higher mode numbers, as illustrated in Figure <ref>. Theoretically, if the values of the HTC are all zero, i.e. if no convective boundary conditions are applied, the KMS reduction basis matches exactly the eigenvalues of the original system. However, including values of the HTC different than zero results in a difference between the reduced and the original eigenvalues due to the bilinearization process. For the thermal FE model under consideration, the maximum relative error between the eigenvalues of the reduced and original system remain below 2 · 10 ^-5. Therefore, the reduced system with parametric convective boundary conditions succeeds in approximating the eigenvalues of the original system for different combinations of values of the HTC.




The parametric KMS method creates a reduced thermal model, which reproduces the dynamics of the thermal behavior. However, the thermally induced structural deformation is the output of interest of the thermo-mechanical models of machine tools. In order to evaluate the mechanical response of the model of Figure <ref>, the mechanical boundary conditions are defined as follows


	
  * The structure is fixed at the location of the guide carriages 
	
  * The stiffness value is 10^10 N/m in all directions   


Equation (<ref>) defines the projection basis for the thermo-mechanical system. The expansion point s_0 for the Krylov subspace of the mechanical system is 30 rad/s. Starting from an original mechanical system with 12471 dofs, the coupling to thermal parametric reduced system results in a thermo-mechanical reduced system of dimension 532.   

In order to evaluate the accuracy of the reduction, the FRF of the thermo-mechanical the original and reduced system are compared. The thermal input is the fluctuation of the environmental temperature, which corresponds to the thermal input shown in Figure <ref>. The mechanical outputs are the linear displacements in X-, Y-, and Z-direction measured at the center of the machine tool table, where the workpiece is placed. The relative error between the two systems are calculated in the frequency range between 10^-5 and 1 rad/s. Figure <ref> shows the FRFs of the relative error between the original and reduced systems for different values of the HTC, considering the displacement in Y-direction as the output of the FRF. Figure <ref> illustrates that the resulting relative errors between reduced and original system remain below 0.01 for up to the maximum frequency of interest, 0.01 rad/s. Thus, the thermo-mechanical coupled reduced system reproduces accurately the response of the original system for different values of the HTC. 



Figure <ref> extends the evaluation to different outputs, namely the deviations in X- , Y-, and Z-direction.  The mechanical response of the system is evaluated for a value of the parameters HTC_top and HTC_bot of 4 and 8 W/m^2K respectively. Figure <ref> illustrates that the  relative error of the mechanical response remains below 0.01 for the frequency range of interest. Therefore, this numerical example illustrates that the reduced thermo-mechanical system captures the response of the original system. 







§ CONCLUSIONS


This work presents a MOR approach to deal with the parametric dependency of the convective boundary conditions in weakly coupled thermo-mechanical models. The KMS method proposes a reduction basis combining the information of the Krylov subspace basis with an expansion point at a low frequency and the thermal eigenmodes. The reduced system reproduces the steady state response of the system as well as the transient behavior up a to a certain frequency of interest. This method is especially suited for thermal models of machine tools and other similar mechatronic systems, which show a decaying amplitude of the response at higher frequencies. This work proposes an upper bound of the magnitude of the reduction error. The error bound provides an a-priori estimation of the KMS reduction error and determines the maximum eigenfrequecy to be included in the reduction basis so that the error remains below a certain value for the frequency range of interest.

Among the different parameters describing thermal models of machine tools, this work focuses on the HTC, which defines the convective heat exchange between the structure and the surrounding fluid media. This paper proposes a parametric MOR method to enable the modification of the value of the HTC after reduction. The KMS reduction basis is extended using the concept of system bilinearization. The numerical results show that the parametrically reduced model reproduces the thermal response of the original system for different values of the HTC. Additionally, the parametric reduction retains the thermal eigenfrequencies of the original model.  

The main objective of thermal error models of machine tools is the prediction of the thermally induced deviations. Therefore, the model needs to evaluate structural deformations associated to inhomogeneous, time-varying temperature distributions. This work presents a method to couple a reduced thermal and mechanical system. The developed coupling method creates a dedicated reduced mechanical system. On one hand, the mechanical system describes the response to any mechanical input, e.g. preloads or gravity, for any combination of the position of the axes. On the other hand, the mechanical system provides the response to any temperature distribution computed by the reduced system. The numerical results show that the reduced system succeeds in capturing the mechanical response of the original system for several values of the HTC and different outputs. The relative error of the reduction remains below 0.01 for the reduced thermo-mechanical system in the frequency range of interest. 

The developed MOR methods are implemented in a simulation software MORe <cit.>, which is designed to facilitate the development of physical models of machine tools. Future work will investigate the theoretical error bounds of the KMS method for arbitrary expansion points as well as generalize the developed error estimator to different input and output combinations. Future research will develop thermo-mechanical models of the whole machine tool assembly with parametric boundary conditions. These models will evaluate the thermo-mechanical response of machine tools to varying conditions in the surrounding environment.  



§ ACKNOWLEDGMENTS


The authors would like to thank the Swiss Innovation Agency (Innosuisse) and the Swiss National Science Foundation (SNSF) for their financial support.



