Abstract

Engineering design often involves representation in at least two levels of abstraction: the system-level, represented by lumped parameter models (LPMs), and the geometric-level, represented by distributed parameter models (DPMs). Functional design innovation commonly occurs at the system-level, followed by a geometric-level realization of functional LPM components. However, comparing these two levels in terms of behavioral outcomes can be challenging and time-consuming, leading to delays in design translations between system and mechanical engineers. In this paper, we propose a simulation-free scheme that compares LPMs and spatially discretized DPMs based on their model specifications and behaviors of interest, regardless of modeling languages and numerical methods. We adopt a model order reduction (MOR) technique that a priori guarantees accuracy, stability, and convergence to improve the computational efficiency of large-scale models. Our approach is demonstrated through the model consistency analysis of several mechanical designs, showing its validity, efficiency, and generality. Our method provides a systematic way to compare system-level and geometric-level designs, improving reliability and facilitating design translation.

Graphical Abstract Figure
Graphical Abstract Figure
Close modal

1 Introduction

In the field of engineering design, the initial step toward creating innovative designs often involves system design, which entails conceptualizing and developing functional structures. Once a system is designed, the subsequent step involves developing its 3D geometry for manufacturing purposes. The objective of the system-based geometric design process is to identify at least one geometric realization that matches the system design. In this paper, we restrict our attention to mechanical systems whose representation at the system-level can be given by mass–spring–damper networks.

1.1 Motivation.

The distributed parameter models (DPMs; e.g., 3D geometric assemblies) and lumped parameter models (LPMs; e.g., mechanical mass–spring–damper networks) that describe the same physical system at different levels of abstraction are represented using different languages and semantics that cannot be directly translated into one another. For instance, Dassault Systémes offers two commercial software, dymola for system modeling and solidworks for geometric modeling. The models created in these programs are incompatible and cannot be automatically translated into each other owing to disparities in model types and representations [14]. This gap presents a significant challenge for ensuring consistency between the system models and computer-aided design/engineering (CAD/CAE) models. This paper concentrates on the crucial technical issue of systematically verifying the consistency between the two models.

LPMs use a network of lumped components to represent the structure of engineering systems, such as masses, springs, and dampers, which are governed by a system of ordinary differential equations (ODEs) in terms of variables that vary with time [5,6]. On the other hand, DPMs explicitly consider the geometric and material properties of engineering systems, which are governed by a system of partial differential equations (PDEs) that take into account variables varying with both time and space [7]. Numerical discretization methods can be used to approximate these PDEs by large systems of ODEs upon approximating the spatial continuum in a finite basis [8].

Figure 1 depicts the process of designing a 3D suspension mechanism based on a system design. After creating an LPM in Modelica [9], mechanical parts are used to realize the lumped components and obtain the expected behaviors. The stiffness and damping of the absorber are modeled by a spring–damper pair, without considering the precise geometric realization. Two potential geometric realization options for the absorber are presented in the library. While there might be various options for the absorber, a designer must verify that the selected option behaves as intended by the lumped component (① in the figure) before the assembly process or design qualified parts guided by the three model consistency conditions explained in Sec. 3. Additionally, it is crucial to confirm that the final geometric assembly behaves as intended by the LPM (depicted as ② in the figure) to ensure a qualified design. Currently, the only reliable way to compare design behaviors is to simulate both the LPM and DPM and compare their differential equation solutions a posteriori, which is computationally prohibitive for large-scale models, not to mention the additional challenges related to selecting appropriate time-steps, stability, and convergence. The goal of this paper is to propose a systematic method to check consistency between the system models and CAD/CAE models.

Fig. 1
The role of model consistency analysis in the process of system-based geometric design
Fig. 1
The role of model consistency analysis in the process of system-based geometric design
Close modal

1.2 Contributions.

Our main technical contributions are:

  1. We propose a generalizable definition of consistency between mechanical LPMs and DPMs that considers both model specifications and solutions, taking into account factors such as mass, initial and boundary conditions (ICs/BCs), and the behavior of interest (BoI).

  2. We develop a simulation-free scheme for checking the consistency between LPMs and DPMs based on the definition proposed above. This idea is to compute a priori error bounds between the solutions of both model types by comparing their parameters, circumventing the costly process of solving differential equations.

  3. We adopt a model order reduction (MOR) technique that a priori guarantees accuracy, stability, and convergence in the simulation-free scheme to improve the efficiency of consistency analysis for large-scale models. The MOR enables the analysis of large-scale designs in a computationally efficient manner.

1.3 Outline.

In Sec. 2, we review the prior works on system-based geometric design, physical model solution comparison, and MOR. In Sec. 3, we define the problem of model consistency analysis. In Sec. 4, we present a simulation-free scheme to compare LPMs and DPMs. We illustrate the application of our scheme to various mechanical problems in Sec. 5. Finally, in Sec. 6, we summarize the implications of our proposed scheme and suggest potential future research directions.

2 Related Work

2.1 System-Based Geometric Design.

The system-to-geometry design strategy has been shown by Ulrich [10] to reduce complexity and clarify the problem-solving process. Most system-based geometric design approaches rely on a one-to-one correspondence between the lumped component in the LPM and the geometric part in a solid model repository. For instance, Finger and Rinderle proposed a component database that contains correspondence between ports of bond graphs and geometric parts for a specific class of mechanical designs [11]. Engelson et al. developed an integrated environment that combines geometric design and system modeling tools to assist engineers in constructing and verifying large, moving rigid-body assemblies [12]. However, these approaches may result in unrealistic designs by ignoring function-sharing. Ulrich’s work stands out as an exception [10], as it demonstrates a systematic approach to merge multiple functions, abstracted by different system-level components, in fewer geometric parts. Despite the numerous approaches proposed for system-based geometric design [10,1216], none of them formally introduce the concept of consistency between system and geometric designs, nor do they provide a systematic approach for assessing the validity of the geometric design with respect to the target system behavior.

2.2 Solution Comparison Between the Lumped Parameter Model and the Distributed Parameter Model.

Comparing the solutions of LPMs and DPMs is a common requirement in model conversion problems, where the two models must be converted to each other with tolerable differences in their respective solutions. In computer graphics, for instance, deformable objects like cloth fabrics and soft tissues are often converted to lumped mass-spring models for faster simulations due to the simplicity and efficiency of LPMs [1719]. Gelder [20] developed a lumped-spring element based on the geometric angle and length information of 2D linear triangular finite elements, while Baudet et al. [21] extended the method to rectangular finite elements. In both cases, the solutions of the converted LPM and the original finite element model can be directly compared. Natsupakpong and Çavuşoğlu [22] proposed a method for converting 3D deformable objects from finite element models to LPMs by minimizing the difference of the stiffness matrices of the two models. However, the dimension of the LPM solution is often much smaller than that of the DPM after spatial discretization, making direct comparison challenging. In other words, a one-to-one correspondence between the solution dimensions of the two models is typically not guaranteed. This limits the applicability of existing methods to the model solution comparison problem addressed in this paper.

2.3 Model Order Reduction.

The field of MOR has seen the development of various techniques [2327] to approximate a given “full-order” model (FOM) in a numerically efficient and stable manner while preserving certain desired properties. The balanced truncation (BT) method, for instance, removes weakly controllable and observable states from the FOM [28], but can be computationally expensive for large-scale models [29]. The rational Krylov subspace (RKS) method, on the other hand, approximates the FOM transfer function by matching a few significant terms of its Taylor series expansion [30], but does not have a standard rule for choosing the frequency, which is the frequency around which the Taylor series expansion is made. To address this limitation, the iterative rational Krylov algorithm (IRKA) was developed, which iteratively updates the expansion frequency using the reflection of poles of the updated reduced order model (ROM) about the imaginary axis at each iteration step until the difference between the transfer functions of FOM and ROM is minimized [31].

However, the IRKA algorithm cannot ensure the error converges to a local minimum [32]. To overcome this limitation, the CUmulative REduction (CURE) scheme was proposed, which adaptively chooses the expansion frequency and incrementally increases the scale of the ROM by monotonically decreasing the norm of the error transfer function to zero through an accumulation process [33]. Furthermore, a stability-preserving, adaptive rational Krylov (SPARK) algorithm [33] was developed to maintain model stability and is usually embedded in the CURE scheme to generate a family of stable ROMs whose orders are increased by sequential accumulation in a single MOR process. This embedded CURE scheme with SPARK algorithm (hereafter abbreviated by SPARK + CURE) was used in Ref. [34] to generate a family of physically interpretable multi-fidelity surrogate LPMs for physical systems governed by PDEs.

The SPARK+CURE method has several key advantages, namely: (1) automatic search for proper expansion frequencies, (2) preserved model stability, (3) guaranteed error convergence, and (4) a priori error bound. This method is particularly suitable for large-scale mechanical systems, as it can significantly improve the time efficiency of model consistency analysis. In Table 1, we compare the important properties of the BT, RKS, IRKA, and SPARK + CURE methods, the advantages of the latter are in bold. For our proposed simulation-free scheme (Sec. 4), we adopt the SPARK + CURE method to ensure accuracy, stability, and convergence in our model consistency analysis.

Table 1

Comparison of four popular MOR methods

BTRKSIRKASPARK+CURE
Numerical efficiencyLowHighHighHigh
A priori error boundYesNoNoYes
Auto order decisionNoNoNoYes
Maintain stabilityYesNoYesYes
Guarantee convergenceYesNoNoYes
BTRKSIRKASPARK+CURE
Numerical efficiencyLowHighHighHigh
A priori error boundYesNoNoYes
Auto order decisionNoNoNoYes
Maintain stabilityYesNoYesYes
Guarantee convergenceYesNoNoYes

3 Problem Formulation

This section outlines the problem formulation for the model consistency analysis procedure, which involves two model types: the LPM, which provides a system-level description of physical systems, and the DPM, which incorporates the spatiotemporal distribution of materials. The definition of these models and the relationship between their model specifications, behaviors, and ICs/BCs are crucial for formulating the problem, and will be discussed below.

3.1 Definitions.

We focus on linear time-invariant (LTI) translational mechanical LPMs that can be modeled as networks of interconnected lumped components such as masses, springs, and dampers. To uniquely determine the system behavior, initial mass displacements, velocities, and source terms must be specified, which give rise to a system of ODEs in the time domain. We assume that any physically meaningful mass–spring–damper network can be realized by at least one DPM, which has a continuous material distribution and geometry embedded in 3D space. The DPM can be divided into pieces such that there is a one-to-one correspondence between the spatial integration of mass density and lumped masses. The effective stiffness and damping of the DPM correspond to those of the mass–spring–damper network, and the spatial integration of ICs/BCs and body effects of different pieces correspond to the ICs and source terms of lumped masses. The DPM behavior is described by PDEs defined over a region of spacetime. To compute the effective stiffness and damping, solving PDEs is usually unavoidable. The solution generally involves spatial discretization and numerical methods such as finite difference, element, and volume analysis [8].

Based on the above definitions, we can define consistency between an LPM and a DPM as the satisfaction of the following three conditions:

  • (C1) The total lumped mass value of the LPM matches the spatial integration of the mass density of the DPM.

  • (C2) The spatial integration of the ICs/BCs, as well as the body effects of the DPM, match the ICs and source terms of the LPM.

  • (C3) The BoI of the DPM matches the BoI of the LPM in terms of spatial integration.

In other words, for an LPM and a DPM to be consistent, the mass, ICs, source terms, and BoI of the LPM must be equivalent to the mass distribution, initial and boundary conditions, and BoI of the DPM, respectively.

3.2 A Mechanical Example.

Figure 2 depicts an example of an LPM and its corresponding DPM realization, where the geometry of a 3D suspension system is considered. The LPM is composed of a mass, two springs, and a damper, and is subjected to a time-varying external force f(t). The network is connected to the ground via a spring. On the other hand, the DPM is an assembly of several solid parts that have linear-elastic material distributed in 3D space. The top surface of the assembly is subjected to time-varying pressure p(x,t), where x is the spatial 3D coordinates. The entire assembly is fixed to the ground. To ensure consistency between the LPM and DPM, we require that the following conditions are satisfied: the mass of the DPM matches the total mass of the LPM, i.e., m=ρ(x)dV where ρ(x) is the material density and dV represents an infinitesimal volume element; the external load on the LPM matches the total force on the DPM, i.e., f(t)=p(x,t)dS where dS represents an infinitesimal surface element; the pre-specified displacement u¯(x,t) of the LPM (which is fixed to the ground) matches that of the DPM, i.e., w¯(t)=(1/Sc)u¯(x,t)dS=0, where Sc is the contact surface area between the tire and the ground; and the solved displacement of the LPM matches that of the DPM, i.e., w(t)=(1/Sp)u(x,t)dS, where Sp is the top surface area of the DPM assembly.

Fig. 2
An LPM, its DPM realization, and the model specification match
Fig. 2
An LPM, its DPM realization, and the model specification match
Close modal

3.3 Formulation.

To analyze consistency between LPM and DPM, one has to verify if the conditions (C1) through (C3) are satisfied.

For practical purposes, we introduce an upper-bound (i.e., inequality constraint) to relax the strict equality constraint of BoI in (C3). Noting that an LPM is a ROM of a DPM, its solution can be at best an approximation of the DPM solution, and hence, it is not reasonable to expect a strict equivalence between the solutions of LPM and DPM.

Given an LPM and a DPM with matched mass properties, ICs/BCs, and source terms, our objective is to check if the overall error between pairs of LPM BoI quantities wi(t) and the corresponding DPM BoI quantities ui(t) for i=1,2,,n is less than an acceptable tolerance ϵ>0. Symbolically, this can be represented as i=1nui(t)wi(t)ϵ. Any commonly used vector norm such as L1, L2, and L can be used as the norm . However, for the sake of specificity, we use L to provide a measure of the maximum deviation between the DPM and LPM BoI quantities.

Below, we present a simulation-free scheme that compares the BoIs of a given pair of LPM and DPM in terms of the L norm, for the purpose of model consistency analysis.

4 Simulation-Free Scheme

Assuming that an LPM and a DPM are provided, along with their matching mass properties, ICs/BCs, and source terms, as well as the BoI quantities’ matching form (i.e., spatial integration over the boundary surface), a scheme is presented below for calculating the a priori error between the corresponding BoI quantities of the two models, without requiring explicit numerical solution of differential equations.

We start by deriving the governing equations for both the LPM and the DPM, which are respectively described by a system of second-order ODEs and PDEs, respectively. To enable a comparison between the two models, we discretize the PDEs of the DPM using spatial discretization methods [35], thereby converting them to a system of second-order ODEs. We then transform the given ODEs of the LPM and the resulting ODEs (semi-discretized PDEs) of the DPM into their respective state-space forms through linear transformation [36], as illustrated below:

  • LPM ODEs in state-space form:
    Elx˙l(t)=Alxl(t)+Blhl(t)yl(t)=Clxl(t)
    (1)
  • DPM ODEs (semi-discretized PDEs) in state-space form:
    Edx˙d(t)=Adxd(t)+Bdhd(t)yd(t)=Cdxd(t)
    (2)

A human-readable LPM ODE system (1) can have tens to hundreds of state variables, while the DPM ODE system (2) can have tens or hundreds of thousands, if not millions, of state variables, as a result of semi-discretizing PDEs on a fine mesh. In both equations, x(t) represents a vector of state variables, y(t) represents a vector of BoIs, and h(t) represents the input signal of external forces. The descriptor matrix E, dynamic matrix A, input matrix B, and output matrix C parameterize the LTI systems. The user determines the matrix C based on the BoI. The matrices A through E have the following forms:
A=[0IKR],B=[0F],C=C,E=[I00M]
(3)
where the matrix I is the identity matrix, M, K, and R are the mass, stiffness, and damping matrices, respectively, and F is a vector that contains the external forces. The matrices M, K, and R used in Al through El in (1) are constructed based on the constitutive relationship of lumped components and their adjacency relations. On the other hand, the matrices M, K, and R used in Ad through Ed in (2) are generated from spatial discretization methods [8].

Since both systems of ODEs are LTI, we can match them using linear projections:

  • The source terms are related by Bl(t)hl(t)=ΓnBd(t)hd(t), where Γn represents a projection matrix that maps the discrete form of forces of DPM to the lumped forces of LPM.

  • The initial displacement and velocity terms are related by ΓIxd(0)=xl(0) and ΓIx˙d(0)=x˙l(0), where ΓI denotes the projection matrix that maps the spatially discretized ICs of the DPM to the ICs of the LPM.

  • The BoI terms are related by Cd=Γf, where Γf is a projection matrix that maps the spatially discretized BoI of DPM to the BoI of LPM.

By substituting the projection in (a) into (1), we obtain a revised state-space form:
Elx˙l(t)=Alxl(t)+Blhd(t)yl(t)=Clxl(t)
(4)
This new form indicates that the source term of the LPM can be replaced with the source term of the DPM. This substitution is crucial for computing an upper-bound for the L error between yl(t) and yd(t), as demonstrated below.
Given that we are working with LTI systems, we can solve (4) and (2) using Laplace transforms [37]. The solutions take the following forms:
yl(s)=(Cl(sElAl)1Bl)Gl(s)hd(s)
(5a)
yd(s)=(Cd(sEdAd)1Bd)Gd(s)hd(s)
(5b)
where s represents the complex frequency variable, while Gl(s) and Gd(s) correspond to the transfer functions of the LPM and semi-discretized DPM, respectively.
When two models with the same input hd(t) are compared, the maximum difference between their outputs yd(t) and yl(t) is upper-bounded [38] as follows:
maxt[0,)yd(t)yl(t)Gd(s)Gl(s)H20hd(t)2dt
(6)
as long as hd(t) is a finite energy input (i.e., 0hd(t)22dt<), which is the case for most engineering problems. The formula to compute Gd(s)Gl(s)H2 is given in (7), where j is the imaginary unit and ω is the frequency in radians per unit time.
Gd(s)Gl(s)H2=(12π|Gd(jω)Gl(jω)|dω)1/2
(7)

Notably, the integration term in (6) is constant for a given source term hd(t), as the BC of DPM is time-invariant. This means that as the norm Gd(s)Gl(s)H2 approaches zero, the maximum difference between yd(t) and yl(t) also approaches zero.

In essence, the deviation between the BoIs of LPM and DPM in the time domain can be approximated by the deviation between Gd(s) and Gl(s) in the frequency domain using the H2 norm. An important observation is that the computation of Gd(s)Gl(s)H2 involves solving algebraic equations instead of ODEs, as discussed in Ref. [39]. This means that we can avoid solving ODEs altogether. To speed up the computation of these algebraic equations, we use the ROM of the DPM, as shown in the equation below.

Erx˙r(t)=Arxr(t)+Brhd(t)yr(t)=Crxr(t)
(8)

In order to enhance the computational efficiency of solving large-scale algebraic equations arising from DPMs, a MOR method is employed. Specifically, we adopt the SPARK + CURE method from Ref. [33], which generates a small-scale surrogate model for the DPM. Due to space constraints, please refer to Panzer’s original Ph.D. thesis [33] for more details about this method. This surrogate model has the same source terms ud(t) as the DPM in (2), and its transfer function is denoted by Gr(s). The SPARK + CURE method guarantees an upper-bound Gd(s)Gl(s)H2ε¯1 between the transfer functions of the DPM and the surrogate model. Unlike the spatial-discretized DPM, we assume that the given LPM in (1) has a small scale, so we do not apply MOR to it.

Because both the given LPM and the surrogate model of the semi-discretized DPM have a small state space, computing the value ε¯2 of Gr(s)Gl(s)H2 is time-efficient. With these two values, we can use a triangular inequality to obtain an upper-bound ε¯ for the error between the transfer functions Gd(s) of the DPM and Gl(s) of the LPM as follows:
Gd(s)Gl(s)H2=Gd(s)Gr(s)+Gr(s)Gl(s)H2Gd(s)Gr(s)H2+Gr(s)Gl(s)H2ε¯1+ε¯2=ε¯
(9)
If the upper-bound ε¯ is within the given tolerance, we can view the BoI of DPM and the BoI of LPM as similar.
The equation above provides an upper-bound for the absolute error. To obtain an upper-bound for the relative error, we can use the equation below, where the transfer function Gl(s) of the LPM is selected as the reference.
Gd(s)Gl(s)H2Gl(s)H2ε¯1+ε¯2Gl(s)H2=ε¯rel
(10)

Figure 3 provides a visual representation of the five-step simulation-free scheme proposed in this study, where an LPM with masses, springs, and dampers and a DPM with 3D suspension are used as examples. The approach is outlined step-by-step as follows. Given the ODEs of the LPM and the PDEs of the DPM:

  1. Use a spatial discretization (e.g., finite element [35]) method to convert the PDEs into a system of ODEs.

  2. Convert the ODEs of both models to state-space form.

  3. Use the SPARK + CURE method to reduce the scale of the DPM and ensure a priori guaranteed H2 error ε¯1 between the DPM and its surrogate model.

  4. Calculate the exact H2 error ε¯2 between the surrogate DPM and the LPM.

  5. Estimate the upper-bound of the H2 error ε¯ between the transfer functions of the given LPM and semi-discretized DPM using the triangular inequality and verify whether it falls within the acceptable tolerance.

Fig. 3
Simulation-free scheme to compute the error bound between the BoIs of a pair of LPM and DPM
Fig. 3
Simulation-free scheme to compute the error bound between the BoIs of a pair of LPM and DPM
Close modal

5 Applications

In this section, we demonstrate the practical application of the proposed simulation-free scheme by using it to compare the solutions of the LPM and the DPM for two mechanical designs: a bracket (Fig. 4(b)) and a frame (Fig. 4(c)). Both designs have linear isotropic material properties, and their mass properties, ICs/BCs, and the form of the BoI quantities are pre-specified to ensure consistency between the two models.

Fig. 4
LPM and DPMs of bracket and frame: (a) LPM, (b) DPM of the bracket, and (c) DPM of the frame
Fig. 4
LPM and DPMs of bracket and frame: (a) LPM, (b) DPM of the bracket, and (c) DPM of the frame
Close modal

5.1 Model Specifications.

The same LPM topology is utilized for both mechanical designs, as shown in Fig. 4(a), but with different parameters for the lumped components. The values are provided in Table 2, using SI base units. Both lumped masses have zero initial displacements and velocities.

Fig. 5
Meshes of bracket and frame: (a) mesh of bracket and (b) mesh of frame
Fig. 5
Meshes of bracket and frame: (a) mesh of bracket and (b) mesh of frame
Close modal
Table 2

Parameters of the LPM used for the bracket and the frame designs

m1m2k1k2r1r2f1
LPM for bracket3.8465×1053.512×1033.316×1044.688×1031.4697×1052.9052×1030.005
LPM for frame7.997×1056.9139×1044.8561×1082.2308×1082.8102×1071.5075×1052186.56
m1m2k1k2r1r2f1
LPM for bracket3.8465×1053.512×1033.316×1044.688×1031.4697×1052.9052×1030.005
LPM for frame7.997×1056.9139×1044.8561×1082.2308×1082.8102×1071.5075×1052186.56

Two mechanical parts, namely, a bracket and a frame, have been designed to implement the LPM instances. The DPMs for both designs have zero initial displacements and velocities. We discretize the bracket and the frame using second-order tetrahedral finite elements (Fig. 5) to generate DPM ODEs. The top surfaces of the bracket and the frame are subjected to pressures that vary with time. The time-varying characteristics of the pressures are depicted in Figs. 6(a) and 6(b) for the bracket and the frame, respectively.

Fig. 6
Pressures on bracket and frame: (a) pressure on bracket and (b) pressure on frame
Fig. 6
Pressures on bracket and frame: (a) pressure on bracket and (b) pressure on frame
Close modal

In order to compare the LPM and DPM models of the two designs, it is necessary to ensure that their mass properties, ICs/BCs, and BoI quantities match. The requirements for this matching are explained in (C1)–(C3) in Sec. 4, and an example of the matching process is provided in Sec. 3.2. To maintain consistency, we will follow the same matching procedure as the example in Sec. 3.2.

5.2 Surrogate Distributed Parameter Model of the Bracket.

Here, we demonstrate the effectiveness of our simulation-free approach in enhancing the efficiency of model comparison for large-scale models. We used tetrahedral finite elements to discretize the bracket’s geometry at three different resolutions (only one sketch is shown in Fig. 7(a)), resulting in three sets of state-space equations with a state variable count ranging from approximately 3000 to 6000. We then applied the SPARK + CURE method to generate three sets of ROMs that have a priori relative H2 error bounds just below 1%.

Fig. 7
Comparison of the simulation results and computation time between DPM and ROMs—the bracket example: (a) case 1, (b) case 2, (c) case 3, and (d) time comparison
Fig. 7
Comparison of the simulation results and computation time between DPM and ROMs—the bracket example: (a) case 1, (b) case 2, (c) case 3, and (d) time comparison
Close modal

To assess the quality of the ROMs, we used the backward Euler method [7] to simulate the DPM and all its surrogate models, using the same solver and time-steps. The simulation results are presented in Figs. 7(a) through 7(c), where the abscissa denotes computation time and the ordinate denotes the average displacement of the bracket’s top surface. Due to limited space, only a few ROM curves are labeled. The solver employed in the simulations was the “sparse state space” toolbox in Matlab [40], which can analyze high-dimensional dynamical systems with state-space dimensions of O(104) or more. This solver preserves system matrix sparsity, allowing it to use computationally efficient operations, such as sparse lower–upper decompositions [41], that would otherwise be infeasible or time-consuming. It is important to note that the simulations were only used to illustrate the accuracy of the surrogate models and were not necessary for our proposed model comparison method.

5.3 Error Analysis.

As the order of ROM increased, which corresponds to an increase in the number of state variables of ROM, the simulated response of the ROM more closely approached that of the DPM for all cases. To quantify the error introduced by the MOR process, we calculated two different measures, as shown in Fig. 8: the relative H2 error of the transfer functions and the root mean square error (RMSE) of the simulated response, whose formulas can be found in Ref. [42]. The SPARK + CURE method guaranteed the strictly monotonic decay of the bound of the a priori relative H2 error, as shown in Figs. 8(d)8(f). However, it did not provide a rigorous a priori bound on RMSE, as shown in Figs. 8(a)8(c), although an overall reduction pattern was commonly observed.

Fig. 8
RMSE and a priori relative H2 error bound—the bracket example: (a) RMSE of case 1, (b) RMSE of case 2, (c) RMSE of case 3, (d) error bound of case 1, (e) error bound of case 2, and (f) error bound of case 3
Fig. 8
RMSE and a priori relative H2 error bound—the bracket example: (a) RMSE of case 1, (b) RMSE of case 2, (c) RMSE of case 3, (d) error bound of case 1, (e) error bound of case 2, and (f) error bound of case 3
Close modal
For all three cases, the RMSE was around O(109), two orders of magnitude smaller than the DPM output of O(107). The a priori relative H2 error bound showed that if the order of ROM was larger than 36, 38, and 34 in the three cases, respectively, the exact relative H2 error would be less than 102.15060.707%, 102.18590.652%, and 102.17970.661%. By selecting the ROM of order 34 generated from case 3 as the surrogate DPM for the bracket, we computed the a priori relative H2 error bound between the discretized DPM and the LPM using (10) as follows:
Gd(s)Gl(s)H2Gl(s)H2ε¯1+ε¯2Gl(s)H2=0.0463=ε¯rel
(11)

Given that the pre-conditions ensure the matches of masses, ICs, and BCs, the LPM and DPM would be considered consistent if the upper-bound of error ε¯=0.0463 is within the user-defined tolerance. To compare the numerical solutions between the LPM and DPM, we plot the errors in Fig. 9(b), while the numerical solutions are compared in Fig. 9(a). The maximum difference between the model solutions over the entire time domain is less than 4.8197×108, which is one order of magnitude smaller than the maximum model output O(106).

Fig. 9
Numerical solution comparison between the LPM and the DPM—bracket example: (a) solution comparison and (b) error
Fig. 9
Numerical solution comparison between the LPM and the DPM—bracket example: (a) solution comparison and (b) error
Close modal

5.4 Time Efficiency Analysis.

Through simulations of the DPM and all ROMs for three cases of the bracket, we have obtained two polynomial curves (Fig. 7(d)) that fit the relation between computation time and the order of DPM for both a posterior solution comparison based on direct ODE simulations and our simulation-free scheme. The curves intersect at a certain point, indicating that if the order of the discretized DPM is lower than 4637, comparing solutions through simulations is more time-efficient. However, if the order is higher than 4637, our simulation-free scheme will be more time-efficient, with an accuracy loss of at most 1%. It should be noted that if a smaller user-selected relative H2 error bound is used (i.e., <1%), a new cross-over point that slightly shifts toward a larger order of DPM will be found.

5.5 Surrogate Distributed Parameter Model of the Frame and Error Analysis.

In a similar manner, we apply the SPARK + CURE approach to the frame design and compare the results of both DPM and ROMs (Fig.10(a)). The computed relative H2-error bound is presented in Fig. 10(b). The figure indicates that if the ROM order is greater than 236, the relative H2 error between the DPM and the surrogate DPM would be less than 102.01680.962%. Based on this, we can compute the a priori H2 relative error bound between the discretized DPM and the LPM as follows:
Gd(s)Gl(s)H2Gl(s)H2ε¯1+ε¯2Gl(s)H2=0.8229=ε¯rel
(12)
Fig. 10
Comparison of the simulation results between DPM and ROM and a priori relative H2 error bound—frame example: (a) solution comparison and (b) relative H2 error bound
Fig. 10
Comparison of the simulation results between DPM and ROM and a priori relative H2 error bound—frame example: (a) solution comparison and (b) relative H2 error bound
Close modal

The agreement between the LPM and DPM can be determined by the user-defined threshold for the relative error. When the threshold is greater than 0.8229, the models are in agreement. If the threshold is smaller than 0.8229, the models are not in agreement. In Fig. 11(a), we compare the simulation results of these two models, and their numerical solutions difference is depicted in Fig. 11(b). It can be observed that the difference is significantly smaller than the displacement of design models in most of the time-steps.

Fig. 11
Numerical solution comparison between LPM and DPM—frame example: (a) solution comparison and (b) error
Fig. 11
Numerical solution comparison between LPM and DPM—frame example: (a) solution comparison and (b) error
Close modal

6 Conclusion

The paper introduces a model consistency concept and presents a simulation-free scheme for comparing the behavior of the LPM and the DPM. By leveraging the SPARK + CURE MOR method, the proposed approach can address the computational challenges associated with large-scale models. Our analysis of two mechanical models demonstrates that our simulation-free scheme provides relatively tight H2 a priori error bounds for behavior comparison. Furthermore, we show that our proposed scheme is more time-efficient than direct numerical simulation if the order of the discretized DPM is above a certain threshold, with increasing efficiency as the DPM order increases. The proposed approach can be applied to various engineering applications to enhance the system-based geometric design and facilitate integration between system-level and geometric designs.

However, we acknowledge that our framework has some limitations. One limitation is that the scheme to compare model solutions does not account for errors caused by spatial discretization of DPMs, and users must select an appropriate spatial discretization method and resolution. Additionally, the CURE scheme can only be applied to reduce the order of dissipative models, and the a priori guarantees provided by SPARK + CURE are only valid for LTI models. There are two ways to extend our proposed scheme to nonlinear models. The first approach is to divide the nonlinear MOR models into piecewise linear MOR models, and then apply the SPARK + CURE method, similar to the trajectory piecewise linear method [43]. The second approach is to use MOR methods designed specifically for nonlinear models, such as the symplectic MOR methods [44], the MOR methods based on the proper orthogonal decomposition [45], the Gaussian process regression [46], etc. However, many of these methods do not offer rigorous a priori error guarantees.

Acknowledgment

This material is based upon work supported by the Defense Advanced Research Projects Agency (DARPA) under Agreement Nos. HR00111990029 and HR00112090065. The responsibility for errors and omissions lies solely with the authors.

Conflict of Interest

There are no conflicts of interest.

Data Availability Statement

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.

References

1.
Systemes
,
D.
,
2018
,
Dymola User Manual
,
Dassault Systèmes
,
Vélizy-Villacoublay, France
.
2.
Chen
,
J.
,
Ilies
,
H. T.
, and
Ding
,
C.
,
2022
, “
Graph-Based Shape Analysis for Heterogeneous Geometric Datasets: Similarity, Retrieval and Substructure Matching
,”
Comput. Aided Des.
,
143
, p.
103125
.
3.
Almattar
,
T.
,
2019
,
Learn SOLIDWORKS 2020: A Hands-On Guide to Becoming an Accomplished SOLIDWORKS Associate and Professional
,
Packt Publishing Ltd
.
4.
Chen
,
J.
, and
Ilieş
,
H. T.
,
2020
, “
Maximal Disjoint Ball Decompositions for Shape Modeling and Analysis
,”
Comput. Aided Des.
,
126
, p.
102850
.
5.
Karnopp
,
D. C.
,
Margolis
,
D. L.
, and
Rosenberg
,
R. C.
,
2012
,
System Dynamics: Modeling, Simulation, and Control of Mechatronic Systems
,
John Wiley & Sons
,
Hoboken, NJ
.
6.
Wang
,
R.
, and
Shapiro
,
V.
,
2019
, “
Topological Semantics for Lumped Parameter Systems Modeling
,”
Adv. Eng. Inform.
,
42
, p.
100958
.
7.
Mazumder
,
S.
,
2015
,
Numerical Methods for Partial Differential Equations: Finite Difference and Finite Volume Methods
,
Academic Press
,
Cambridge, MA
.
8.
Mattiussi
,
C.
,
2000
, “The Finite Volume, Finite Element, and Finite Difference Methods as Numerical Methods for Physical Field Problems,”
Advances in Imaging and Electron Physics
,
Elsevier
,
Amsterdam, Netherlands
, pp.
1
146
.
9.
Fritzson
,
P.
,
2011
,
Introduction to Modeling and Simulation of Technical and Physical Systems With Modelica
,
John Wiley & Sons
,
Hoboken, NJ
.
10.
Ulrich
,
K. T.
,
2003
,
Product Design and Development
,
Tata McGraw-Hill Education
,
New York
.
11.
Finger
,
S.
, and
Rinderle
,
J.
,
1990
,
A Transformational Approach to Mechanical Design Using a Bond Graph Grammar
,
Carnegie Mellon University, Engineering Design Research Center
,
Pittsburgh, PA
.
12.
Engelson
,
V.
,
Bunus
,
P.
,
Popescu
,
L.
, and
Fritzson
,
P.
,
2003
, “
Mechanical CAD With Multibody Dynamic Analysis Based on Modelica Simulation
,”
44th Scandinavian Conference on Simulation and Modeling
,
Mälardalen University in Västerås, Sweden
,
Sept. 18–19
.
13.
Prabhu
,
D.
, and
Taylor
,
D.
,
1989
, “
Synthesis of Systems From Specifications Containing Orientations and Positions Associated With Flow Variables
,”
1989 ASME Design Automation Conference
,
Montreal, Quebec, Canada
,
Sept. 17–21
, pp.
273
280
.
14.
Greer
,
J. L.
,
2002
, “Effort Flow Analysis: A Methodology for Directed Product Evolution Using Rigid Body and Compliant Mechanisms,” PhD thesis, The University of Texas at Austin.
15.
Kota
,
S.
, and
Ananthasuresh
,
G.
,
1995
, “
Designing Compliant Mechanisms
,”
Mech. Eng. CIME
,
117
(
11
), pp.
93
97
.
16.
Greer
,
J. L.
,
Jensen
,
D. D.
, and
Wood
,
K. L.
,
2004
, “
Effort Flow Analysis: A Methodology for Directed Product Evolution
,”
Des. Stud.
,
25
(
2
), pp.
193
214
.
17.
Mollemans
,
W.
,
Schutyser
,
F.
,
Van Cleynenbreugel
,
J.
, and
Suetens
,
P.
,
2004
, “
Fast Soft Tissue Deformation With Tetrahedral Mass Spring Model for Maxillofacial Surgery Planning Systems
,”
International Conference on Medical Image Computing and Computer-Assisted Intervention
,
Saint-Malo, France
,
Sept. 26–29
, Springer, pp.
371
379
.
18.
Baraff
,
D.
, and
Witkin
,
A.
,
1998
, “
Large Steps in Cloth Simulation
,”
25th Annual Conference on Computer Graphics and Interactive Techniques
,
Orlando, FL
,
July 19–24
, pp.
43
54
.
19.
Kähler
,
K.
,
Haber
,
J.
, and
Seidel
,
H.-P.
,
2001
, “Geometry-Based Muscle Modeling for Facial Animation,”
Graphics Interface
, Vol. 2001, pp.
37
46
.
20.
Gelder
,
A. V.
,
1998
, “
Approximate Simulation of Elastic Membranes by Triangulated Spring Meshes
,”
J. Graph. Tools
,
3
(
2
), pp.
21
41
.
21.
Baudet
,
V.
,
Beuve
,
M.
,
Jaillet
,
F.
,
Shariat
,
B.
, and
Zara
,
F.
,
2007
, “New Mass-Spring System Integrating Elasticity Parameters in 2d.”
22.
Natsupakpong
,
S.
, and
Çavuşoğlu
,
M. C.
,
2010
, “
Determination of Elasticity Parameters in Lumped Element (Mass–Spring) Models of Deformable Objects
,”
Graph. Models
,
72
(
6
), pp.
61
73
.
23.
Baur
,
U.
,
Benner
,
P.
, and
Feng
,
L.
,
2014
, “
Model Order Reduction for Linear and Nonlinear Systems: A System-Theoretic Perspective
,”
Arch. Comput. Methods Eng.
,
21
(
4
), pp.
331
358
.
24.
Fang
,
D.
,
Huang
,
Z.
,
Zhang
,
J.
,
Hu
,
Z.
, and
Tan
,
J.
,
2022
, “
Flow Pattern Investigation of Bionic Fish by Immersed Boundary–Lattice Boltzmann Method and Dynamic Mode Decomposition
,”
Ocean Eng.
,
248
, p.
110823
.
25.
Qu
,
Z.-Q.
,
2004
,
Model Order Reduction Techniques With Applications in Finite Element Analysis: With Applications in Finite Element Analysis
,
Springer Science & Business Media
,
Berlin, Germany
.
26.
Huang
,
B.
, and
Wang
,
J.
,
2022
, “
Applications of Physics-Informed Neural Networks in Power Systems—A Review
,”
IEEE Trans. Power Syst.
,
38
(
1
), pp.
572
588
.
27.
Fang
,
D.
,
Zhang
,
J.
, and
Huang
,
Z.
,
2023
, “
Modal Analysis on Mechanism of Bionic Fish Swimming by Dynamic Mode Decomposition
,”
Ocean Eng.
,
273
, p.
113897
.
28.
Chahlaoui
,
Y.
,
Gallivan
,
K. A.
,
Vandendorpe
,
A.
, and
Van Dooren
,
P.
,
2005
, “Model Reduction of Second-Order Systems,”
Dimension Reduction of Large-Scale Systems
,
Springer
,
Amsterdam, The Netherlands
, pp.
149
172
.
29.
Besselink
,
B.
,
Tabak
,
U.
,
Lutowska
,
A.
,
Van de Wouw
,
N.
,
Nijmeijer
,
H.
,
Rixen
,
D. J.
,
Hochstenbach
,
M.
, and
Schilders
,
W.
,
2013
, “
A Comparison of Model Reduction Techniques From Structural Dynamics, Numerical Mathematics and Systems and Control
,”
J. Sound Vib.
,
332
(
19
), pp.
4403
4422
.
30.
Lohmann
,
B.
, and
Salimbahrami
,
B.
,
2000
, “
Introduction to Krylov Subspace Methods in Model Order Reduction
,”
Colloquium of Automation
,
Leer, Germany
.
31.
Gugercin
,
S.
,
Antoulas
,
A. C.
, and
Beattie
,
C.
,
2008
, “
H2 Model Reduction for Large-Scale Linear Dynamical Systems
,”
SIAM J. Matrix Anal. Appl.
,
30
(
2
), pp.
609
638
.
32.
Beattie
,
C. A.
, and
Gugercin
,
S.
,
2009
, “
A Trust Region Method for Optimal H 2 Model Reduction
,”
48h IEEE Conference on Decision and Control (CDC) Held Jointly with 2009 28th Chinese Control Conference
,
Shanghai, China
,
Dec. 15–18
, IEEE, pp.
5370
5375
.
33.
Panzer
,
H. K.
,
2014
, “Model Order Reduction by Krylov Subspace Methods With Global Error Bounds and Automatic Choice of Parameters,” PhD thesis, Technische Universität München.
34.
Wang
,
R.
, and
Behandish
,
M.
,
2021
, “
Surrogate Modeling for Physical Systems with Preserved Properties and Adjustable Tradeoffs
,”
AAAI 2021 Spring Symposium on Combining Artificial Intelligence and Machine Learning with Physics Sciences
,
Stanford, CA
.
35.
Bathe
,
K.-J.
,
2006
,
Finite Element Procedures
,
Prentice Hall
,
Hoboken, NJ
.
36.
Atkinson
,
K.
,
Han
,
W.
, and
Stewart
,
D. E.
,
2011
,
Numerical Solution of Ordinary Differential Equations
,
John Wiley & Sons
,
Hoboken, NJ
.
37.
Callier
,
F. M.
, and
Desoer
,
C. A.
,
2012
,
Linear System Theory
,
Springer Science & Business Media
,
Berlin, Germany
.
38.
Antoulas
,
A. C.
,
Beattie
,
C. A.
, and
Gugercin
,
S.
,
2010
, “Interpolatory Model Reduction of Large-Scale Dynamical Systems,”
Efficient Modeling and Control of Large-Scale Systems
,
Springer
,
Amsterdam, The Netherlands
, pp.
3
58
.
39.
Peeters
,
J.
, and
Michiels
,
W.
,
2013
, “
Computing the H2 Norm of Large-Scale Time-Delay Systems
,”
IFAC Proc. Vol.
,
46
(
3
), pp.
114
119
.
40.
Castagnotto
,
A.
,
Varona
,
M. C.
,
Jeschek
,
L.
, and
Lohmann
,
B.
,
2017
, “
sss & sssmor: Analysis and Reduction of Large-Scale Dynamic Systems in Matlab
,”
At-Automatisierungstechnik
,
65
(
2
), pp.
134
150
.
41.
Trefethen
,
L. N.
, and
Bau
,
D.
,
2022
,
Numerical Linear Algebra
,
SIAM
,
Philadelphia, PA
.
42.
Wang
,
R.
,
2021
,
Consistency Analysis Between Lumped and Distributed Parameter Models
,
The University of Wisconsin-Madison
.
43.
Rewieński
,
M.
, and
White
,
J.
,
2006
, “
Model Order Reduction for Nonlinear Dynamical Systems Based on Trajectory Piecewise-Linear Approximations
,”
Linear Algebra Appl.
,
415
(
2–3
), pp.
426
454
.
44.
Peng
,
L.
, and
Mohseni
,
K.
,
2016
, “
Symplectic Model Reduction of Hamiltonian Systems
,”
SIAM J. Sci. Comput.
,
38
(
1
), pp.
A1
A27
.
45.
Chaturantabut
,
S.
, and
Sorensen
,
D. C.
,
2010
, “
Nonlinear Model Reduction Via Discrete Empirical Interpolation
,”
SIAM J. Sci. Comput.
,
32
(
5
), pp.
2737
2764
.
46.
Marzouk
,
Y. M.
, and
Najm
,
H. N.
,
2009
, “
Dimensionality Reduction and Polynomial Chaos Acceleration of Bayesian Inference in Inverse Problems
,”
J. Comput. Phys.
,
228
(
6
), pp.
1862
1902
.