Efficient calculation of higher order time history response derivatives by substructuring method

Abstract

Time history response derivatives with respect to the design variables are frequently used in optimization design, damage detection, structural control, etc. This paper proposes a substructuring method for efficient calculation of higher order time history response derivatives of large-scale structures. First, the global structure is disassembled into several small substructures and the substructural displacement is projected onto the range space of a few substructural master eigenvectors. Afterwards, the derivative of substructural master eigenvectors of a few substructures containing the design variables are assembled to form the reduced first, second, and higher order sensitivity equations with a size of the number of master eigenvectors. The equivalent eigenvector which relates the slave and master eigenvectors is derived to compensate for the inertial effect of discarded slave eigenvectors. Finally, the first, second, and higher order time history response derivatives of global large-scale structure are efficiently solved from the reduced sensitivity equations by using Newmark-β method. A numerical one-bay plane frame and a numerical highway bridge are applied to verify the accuracy and efficiency of proposed substructuring method.

Keywords

large-scale structure substructuring method time history response derivative master eigenvector derivative reduced sensitivity equation

Introduction

To investigate the influence of the change of design variables on the structural response, Taylor approximation based sensitivity analysis is used. The first, second, and higher order response derivatives with respect to the design variables are calculated to form the coefficients of Taylor approximation. In the field of structural health monitoring, structural damage identification or finite element model updating, the response derivatives are also used to form the gradient and Hessian matrices to indicate the searching direction (Hou and Xia, 2021; Lin et al., 2020; Mottershead et al., 2011; Weng et al., 2011). There are several methods which have been proposed to calculate the first, second, and higher order response derivatives, such as finite difference method, continuous sensitivity (Ilinca et al., 2007), discrete sensitivity method (Chung et al., 2009), complex variable method (Gomez-Farias et al., 2015), etc. Azqandi and Hassanzadeh (2021) proposed the extended complex variables method to calculate the first and second order derivatives of beam tip displacement with respect to the beam length, where the calculation of second order derivatives is independent of the step size. Lei et al. (2019) proposed a Lagrange direct method to calculate the second order derivative of modal strain energy of undamped structures. To construct the gradient and Hessian matrices, Zhang and Yu (2020) proposed a method to calculate the first and second order derivatives of the complex mode for an asymmetrical damped linear system, and the proposed method was compared with algebraic method, direct method, full-mode method, and truncated modal method (Zhang et al., 2023). In sensitivity based methods, the response derivatives of the analytical model are usually required to be calculated repeatedly. For the structure with tons of degrees of freedom (dofs) and design variables, calculation of the response derivatives of large-scale structure is computationally expensive (Weng et al., 2017).

Substructuring method is advantageous to analyse large-scale structure because it analyses the structure in a piece-wise manner rather than as a whole (De Klerk et al., 2008). Substructuring decoupling method decouples the target substructure from the global structure and the properties of global structure are calculated from the isolated substructure by using the traditional global method (Jalali and Rideout, 2022; Yuen and Huang, 2018; Zhang et al., 2018). Substructuring coupling method calculates the properties of global structure by assembling those of independent substructures by applying constraints on the interface between adjacent substructures (Craig and Bampton, 1968; Hurty, 1965; MacNeal, 1971; Rubin, 1975). Chen et al. (2021) proposed a substructuring coupling method for efficient vibration analysis of cables in cable-stayed bridges, where each cable is regarded as a substructure and the damper as massless interface elements. Pradeepkumar and Nagaraj (2022) applied the substructuring coupling method to reduce the size of complex wind turbine blade model from 8937 super-elements to 12 super-elements. Voormeeren and Rixen (2012) reviewed a family of substructuring decoupling method based on a dual assembly. Wang et al. (2020) predicted the frequency response functions of nonlinear substructure from those of the global structure by using the substructuring decoupling method.

Substructuring method has been combined with sensitivity analysis to alleviate the computational burden of calculating the structural response derivatives. Substructuring decoupling based sensitivity analysis method is widely used for system identification and damage detection (Weng et al., 2020). Koh et al. (1991) firstly applied Kalman filter method to identify substructural variables, in which the interface force is treated as unknown input. Lei et al. (2015) applied extended Kalman filter method to simultaneously identify the stiffness variables, unmeasured ground motion, and unknown interface force of flexible buildings with bending deformation. The time history response derivatives of target substructure are iteratively calculated to indicate searching direction for a weighted global iteration algorithm (Koh et al., 1991; Lei et al., 2015). Weng et al. (2013b) isolated the target substructure from the global structure through the force equilibrium and displacement compatibility, and the eigensensitivity based finite element model updating was performed on the target substructure. Zhu et al. (2013) combined substructuring method and response sensitivity based finite element model updating method to identity interface force and structural damage under moving force, in which the unknown interface force is approximated as Chebyshev polynomials. Hou et al. (2015) proposed a frequency-domain method of substructure identification for local health monitoring using substructure isolation method, where the virtual force distortions are used to model interfaces. Li and Law (2012) proposed a substructural damage identification method without the information of interface forces and interface responses, where the response sensitivity method is used to update the isolated substructure.

The component mode synthesis (CMS) is the most well-known substructuring coupling method, which combines piecewise analysis and model reduction technique. CMS method is classified into fixed interface CMS method (Craig and Bampton, 1968; Hurty, 1965), free interface CMS method (MacNeal, 1971; Rubin, 1975), and mixed interface CMS method (Liew et al., 1996). Kron’s substructuring method is a kind of free interface CMS method, which assembles substructures in a dual way and uses the complete set of free interface modes to approximate substructural displacement (Kron, 1963). Kron’s substructuring method is advantageous to analyze large-scale structures with complex constraints because the free interface modes are calculated directly from substructural stiffness and mass matrices without adding additional constraints on substructures. Kron’s substructuring method has been applied to calculate structural response (Tian et al., 2021), response sensitivity (Weng et al., 2013a), sensitivity-based model updating (Weng et al., 2011), and damage identification (Zhu et al., 2021).

In the existing CMS method, the reduced eigenequation for eigensolution is derived by describing the dynamic response with a few lower order modes (master modes) and the contribution of discarded modes (slave modes) is compensated by residual flexibility. The conventional CMS method only considers the first-order residual flexibility and neglects higher-order residual flexibility (Park and Park, 2004). Weng et al. (2013a) uses the second-order residual flexibility to further consider the inertial effect related to discarded modes, where the unknown eigenvalue-dependent parameter is regarded as an additional coordinate. Kim and Lee (2015) and Kim et al. (2017) use the higher-order residual flexibility to consider the contribution of discarded modes, thus improving the accuracy of the Craig-Bampton and dual Craig-Bampton methods for eigensolutions. Go et al. (2020) proposed a family of Craig–Bampton methods considering the higher order residual flexibility for discarded mode compensation by using O’Callahan’s approximation or adding generalized coordinate vectors containing unknown eigenvalues. An iterative procedure is also employed to consider the contribution of discarded modes accurately (Chung et al., 2021). The aforementioned research focuses on improving the accuracy of substructuring method to calculate eigensolutions and eigensensitivity. Few research focuses on the substructuring method to calculate time history response and response sensitivity with accurate compensation of discarded modes (Gruber, 2019). It is difficult to extend the concept of higher-order residual flexibility to the substructuring method in calculation of time history response and response sensitivity because the damping matrix and external force need to be reduced for the time domain case.

This paper proposes a Kron’s substructuring method to calculate the first, second, and higher order time history response derivatives with respect to the design variables. First, the global structure is disassembled into several small substructures. The derivatives of a few lower eigenvectors with respect to the design variables are calculated for substructures containing the design variables. Afterwards, the substructural eigenvector derivatives are assembled to form the reduced first, second, and higher order sensitivity equations. The equivalent eigenvector is derived to relate substructural lower and higher eigenvectors and equivalent eigenvector related matrix is used to compensate for the inertial effect of discarded higher eigenvectors. Finally, the time history response derivatives of global structure are efficiently solved from the reduced sensitivity equations by using Newmark-β method (Newmark, 1959). As the eigenvector derivatives of substructures that do not contain the design variables are zeroes, the time history response derivatives of global structure are determined by the eigenvector derivatives of a few substructures that include the design variables. Besides, only a small number of substructural eigenvectors are required since the contribution of discarded substructural eigenvectors are accurately considered. The above features ensure the time efficiency and accuracy of the proposed substructuring method to the calculation of first, second, and higher order time history response derivatives meanwhile. The proposed substructuring method is applied to a numerical one-bay plane frame and a numerical highway bridge to verify its accuracy and efficiency.

Substructuring method for time history response

The proportionally damped global structure is divided into N_s substructures. Projecting substructural time history displacement onto the modal space gets

x {(t)}^{p} = [Φ^{p}] q (t)

(1)

where

x {(t)}^{p} = {[x {(t)}^{(1)}, . . ., x {(t)}^{(i)}, . . ., x {(t)}^{(N_{s})}]}^{T}

(1a)

Φ^{p} = D i a g [Φ^{(1)}, Φ^{(2)}, . . ., Φ^{(i)}, . . ., Φ^{(N_{s})}], Φ^{(i)} = [ϕ_{1}^{(i)}, ϕ_{2}^{(i)}, . . ., ϕ_{j}^{(i)}, . . ., ϕ_{n_{i}}^{(i)}]

(1b)

x ( t )^{(
i
)} represents displacement of the ith substructure, q ( t ) stands for the modal coordinate, and

ϕ_{j}^{(i)}

is the jth eigenvector of the ith substructure. t refers to the vibration time.

Applying constraints on the interface dofs between adjacent substructures gets the assembled motion equation of global structure in modal domain (Zhu et al., 2021)

[\begin{array}{l} I^{p} & 0 \\ 0 & 0 \end{array}] {\begin{cases} \ddot{q} (t) \\ \ddot{μ} \end{cases}} + [\begin{array}{c} a_{1} I^{p} + a_{2} Λ^{p} & 0 \\ 0 & 0 \end{array}] {\begin{cases} \dot{q} (t) \\ \dot{μ} \end{cases}} + [\begin{array}{c} Λ^{p} & - {[D Φ^{p}]}^{T} \\ - D Φ^{p} & 0 \end{array}] {\begin{cases} q (t) \\ μ \end{cases}} = {\begin{array}{c} {[Φ^{p}]}^{T} f {(t)}_{ext}^{p} \\ 0 \end{array}}

(2)

where

f {(t)}_{ext}^{p} = {[f {(t)}_{ext}^{(1)}, . . ., f {(t)}_{ext}^{(i)}, . . ., f {(t)}_{ext}^{(N_{s})}]}^{T}

(2a)

Λ^{p} = D i a g [Λ^{(1)}, Λ^{(2)}, . . ., Λ^{(i)}, . . ., Λ^{(N_{s})}], Λ^{(i)} = D i a g [λ_{1}^{(i)}, λ_{2}^{(i)}, . . ., λ_{j}^{(i)}, . . ., λ_{n_{i}}^{(i)}]

(2b)

a₁ and a₂ represent Rayleigh damping coefficients of mass and stiffness matrices, respectively.

λ_{j}^{(i)}

is the jth eigenvalue of the ith substructure, and

f {(t)}_{ext}^{(i)}

represents the external excitation on the ith substructure. D represents the connection matrix indicating the connected relation between interface dofs. μ stands for the intensity of interface force.

This motion equation in dual form is rank deficient and cannot be solved directly. Besides, calculating all eigenpairs (eigenvectors and eigenvalues) of each substructure would consume too much time. Accordingly, the substructural eigenpairs are divided into master and slave ones herein, and the master eigenpairs are calculated and the contribution of slave eigenpairs is compensated by master eigenpair related term (Zhu et al., 2021). The master and eigenpairs are represented as

Λ_{s}^{p} = D i a g [Λ_{s}^{(1)}, . . ., Λ_{s}^{(i)}, . . ., Λ_{s}^{(N_{s})}], Λ_{s}^{(i)} = D i a g [λ_{m_{i} + 1}^{(i)}, λ_{m_{i} + 2}^{(i)}, . . ., λ_{n_{i}}^{(i)}]

(3)

Φ_{s}^{p} = D i a g [Φ_{s}^{(1)}, . . ., Φ_{s}^{(i)}, . . ., Φ_{s}^{(N_{s})}], Φ_{s}^{(i)} = [ϕ_{m_{i} + 1}^{(i)}, ϕ_{m_{i} + 2}^{(i)}, . . ., ϕ_{n_{i}}^{(i)}]

(4)

Λ_{m}^{p} = D i a g [Λ_{m}^{(1)}, . . ., Λ_{m}^{(i)}, . . ., Λ_{m}^{(N_{s})}], Λ_{m}^{(i)} = D i a g [λ_{1}^{(i)}, λ_{2}^{(i)}, . . ., λ_{m_{i}}^{(i)}]

(5)

Φ_{m}^{p} = D i a g [Φ_{m}^{(1)}, . . ., Φ_{m}^{(i)}, . . ., Φ_{m}^{(N_{s})}], Φ_{m}^{(i)} = [ϕ_{1}^{(i)}, ϕ_{2}^{(i)}, . . ., ϕ_{m_{i}}^{(i)}]

(6)

The master eigenpairs are the first lower order ones, denoted as subscript m, and the slave eigenpairs are the remaining higher order ones, denoted as subscript s. m_i is the number of master eigenpairs. n_i stands for the number of dofs of ith substructure.

Expanding the motion equation of global structure (equation (2)) with respect to the master and slave eigenpairs has

[\begin{array}{l} I_{m}^{p} & 0 & 0 \\ 0 & I_{s}^{p} & 0 \\ 0 & 0 & 0 \end{array}] {\begin{cases} \ddot{q} {(t)}_{m} \\ \ddot{q} {(t)}_{s} \\ \ddot{μ} \end{cases}} + [\begin{array}{c} a_{1} I_{m}^{p} + a_{2} Λ_{m}^{p} & 0 & 0 \\ 0 & a_{1} I_{s}^{p} + a_{2} Λ_{s}^{p} & 0 \\ 0 & 0 & 0 \end{array}] {\begin{cases} \dot{q} {(t)}_{m} \\ \dot{q} {(t)}_{s} \\ \dot{μ} \end{cases}} + [\begin{array}{l} Λ_{m}^{p} & 0 & - {(D Φ_{m}^{p})}^{T} \\ 0 & Λ_{s}^{p} & - {(D Φ_{s}^{p})}^{T} \\ - D Φ_{m}^{p} & - D Φ_{s}^{p} & 0 \end{array}] {\begin{cases} q {(t)}_{m} \\ q {(t)}_{s} \\ μ \end{cases}} = {\begin{array}{c} {[Φ_{m}^{p}]}^{T} f {(t)}_{ext}^{p} \\ {[Φ_{s}^{p}]}^{T} f {(t)}_{ext}^{p} \\ 0 \end{array}}

(7)

where q ( t )_m and q ( t )_s are the master and slave modal coordinates, respectively.

According to the second line of equation (7), the motion equation with respect to slave eigenvectors is

\ddot{q} {(t)}_{s} + (a_{1} I_{s}^{p} + a_{2} Λ_{s}^{p}) \dot{q} {(t)}_{s} + Λ_{s}^{p} q {(t)}_{s} = {[Φ_{s}^{p}]}^{T} D^{T} μ + {[Φ_{s}^{p}]}^{T} f {(t)}_{ext}^{p}

(8)

In traditional substructuring method, the inertial and damping terms of slave eigenvectors are neglected, that is, the terms (

\ddot{q} {(t)}_{s}

and

(a_{1} I_{s}^{p} + a_{2} Λ_{s}^{p}) \dot{q} {(t)}_{s}

) are assumed to be zeroes (Zhu et al., 2019). In this paper, the inertial and damping terms of slave eigenvectors are retained, and μ is solved from equation (8) as

μ = {(D F^{p} D^{T})}^{- 1} {D F^{p} M^{p} Φ_{s}^{p} \ddot{q} {(t)}_{s} + (a_{1} D F^{p} M^{p} + a_{2} D) Φ_{s}^{p} \dot{q} {(t)}_{s} + D Φ_{s}^{p} q {(t)}_{s} - D F^{p} f {(t)}_{ext}^{p}}

(9)

Assuming the acceleration, velocity, and external force in the motion equation to be zeroes, the static response is solved from the following static motion equation as

[\begin{array}{l} Λ_{m}^{p} & 0 & - {(D Φ_{m}^{p})}^{T} \\ 0 & Λ_{s}^{p} & - {(D Φ_{s}^{p})}^{T} \\ - D Φ_{m}^{p} & - D Φ_{s}^{p} & 0 \end{array}] {\begin{cases} q {(t)}_{m} \\ q {(t)}_{s} \\ μ \end{cases}} = {\begin{cases} 0 \\ 0 \\ 0 \end{cases}}

(10)

According to the second and third lines of equation (10), the master and slave static responses satisfy the relation of

Φ_{s}^{p} q {(t)}_{s} = - Θ_{s} q {(t)}_{m}

(11)

where

Θ_{s} = F^{p} D^{T} {(D F^{p} D^{T})}^{- 1} D Φ_{m}^{p}, F^{p} = Φ_{s}^{p} {[Λ_{s}^{p}]}^{- 1} {[Φ_{s}^{p}]}^{T} = {(K^{p})}^{- 1} - Φ_{m}^{p} {[Λ_{m}^{p}]}^{- 1} {[Φ_{m}^{p}]}^{T}

(11a)

K^{p} = D i a g [K^{(1)}, . . ., K^{(i)}, . . ., K^{(N_{s})}]

(11b)

F^p is the first order residual flexibility matrix. Θ_s relates the master and slave eigenvectors, which is named as the equivalent eigenvector herein. K⁽ⁱ⁾ represents stiffness matrix of the ith substructure.

Substituting equation (11) into equation (9) gets

μ = - {(D F^{p} D^{T})}^{- 1} {D F^{p} M^{p} Θ_{s} \ddot{q} {(t)}_{m} + a_{1} D F^{p} M^{p} Θ_{s} \dot{q} {(t)}_{m} + a_{2} D Φ_{m}^{p} \dot{q} {(t)}_{m} + D Φ_{m}^{p} q {(t)}_{m} + D F^{p} f {(t)}_{ext}^{p}}

(12)

in which μ is only related to master eigenvectors. Substituting equation (12) into the motion equation with respect to master eigenvectors, i.e., the first line of equation (7), leads to the reduced motion equation

M_{d} \ddot{q} {(t)}_{m} + C_{d} \dot{q} {(t)}_{m} + K_{d} q {(t)}_{m} = f_{d} (t)

(13)

where

M_{d} = I_{m}^{p} + Θ_{s}^{T} M^{p} Θ_{s}, C_{d} = a_{1} (I_{m}^{p} + Θ_{s}^{T} M^{p} Θ_{s}) + a_{2} (Λ_{m}^{p} + Z_{m}), K_{d} = Λ_{m}^{p} + Z_{m}

(13a)

f_{d} (t) = ({[Φ_{m}^{p}]}^{T} - Θ_{s}^{T}) f {(t)}_{ext}^{p}, Z_{m} = {(D Φ_{m}^{p})}^{T} {(D F^{p} D^{T})}^{- 1} D Φ_{m}^{p}

(13b)

M_d, C_d, K_d, and f _d( t ) represent the reduced mass matrix, damping matrix, stiffness matrix, and the external excitation. q ( t )_m,

\dot{q} {(t)}_{m}

and

\ddot{q} {(t)}_{m}

are then calculated from the reduced motion equation (Equation (13)) using Newmark-β method (Newmark, 1959). Afterwards, combining Equation (11) and Equation (1) gets the structural response as

x {(t)}^{p} = [Φ_{m}^{p}] q {(t)}_{m} + [Φ_{s}^{p}] q {(t)}_{s} = (Φ_{m}^{p} - Θ_{s}) q {(t)}_{m}

(14)

\dot{x} {(t)}^{p} = (Φ_{m}^{p} - Θ_{s}) \dot{q} {(t)}_{m}, \ddot{x} {(t)}^{p} = (Φ_{m}^{p} - Θ_{s}) \ddot{q} {(t)}_{m}

(15)

The reduced motion equation is only related to the master eigenvectors, and the slave eigenvectors are not required to be calculated. In addition, the inertial effect of discarded slave eigenvectors are considered through adding the term ( $Θ_{s}^{T} M^{p} Θ_{s}$ ) in the reduced mass matrix, which ensures the accuracy of the proposed substructuring method with a few master eigenvectors retained. The size of the reduced motion equation is equal to the number of master eigenvectors, which means the size of motion equation is reduced a lot. Accordingly, the proposed substructuring method can calculate the time history response accurately and efficiently.

First order time history response derivatives

The design variable is chosen as the elemental stiffness factor herein. The reduced first order sensitivity equation is derived by differentiating the reduced motion equation (Equation (13)) with respect to the design variable r_j as

M_{d} \frac{\partial \ddot{q} {(t)}_{m}}{\partial r_{j}} + C_{d} \frac{\partial \dot{q} {(t)}_{m}}{\partial r_{j}} + K_{d} \frac{\partial q {(t)}_{m}}{\partial r_{j}} = \frac{\partial f_{d} (t)}{\partial r_{j}} - \frac{\partial M_{d}}{\partial r_{j}} \ddot{q} {(t)}_{m} - \frac{\partial C_{d}}{\partial r_{j}} \dot{q} {(t)}_{m} - \frac{\partial K_{d}}{\partial r_{j}} q {(t)}_{m}

(16)

According to equations (13a) and (13b), the reduced system matrices (M_d, C_d, K_d and f _d( t )) are expressed in terms of master eigenpairs, mass matrix, and stiffness matrix. Because the design variable is not related to mass matrix, the first order derivative of mass matrix with respect to the design variable is zero. Therefore, the first order derivative of reduced system matrices (∂M_d/∂r_j, ∂C_d/∂r_j, ∂K_d/∂r_j and ∂ f _d( t )/∂r_j) with respect to the design variable is calculated from those of master eigenpairs ( $\partial Φ_{m}^{p} / \partial r_{j}$ ; $\partial Λ_{m}^{p} / \partial r_{j}$ ) and stiffness matrix. Those sensitivity matrices of one substructure (for example, the Qth substructure) that includes the design variable are calculated and those of other substructures are zeroes. Accordingly, the first order assembled eigenpair derivatives ( $\partial Φ_{m}^{p} / \partial r_{j}$ ; $\partial Λ_{m}^{p} / \partial r_{j}$ ) are formed as

\frac{\partial Λ_{m}^{p}}{\partial r_{j}} = [\begin{array}{c} 0 & 0 & 0 \\ 0 & \frac{\partial Λ_{m}^{(Q)}}{\partial r_{j}} & 0 \\ 0 & 0 & 0 \end{array}], \frac{\partial [Φ_{m}^{p}]}{\partial r_{j}} = [\begin{array}{c} 0 & 0 & 0 \\ 0 & \frac{\partial [Φ_{m}^{(Q)}]}{\partial r_{j}} & 0 \\ 0 & 0 & 0 \end{array}]

(17)

The first order eigenpair derivatives of Qth substructure ( $\partial Φ_{m}^{(Q)} / \partial r_{j}$ ; $\partial Λ_{m}^{(Q)} / \partial r_{j}$ ) are calculated by Nelson’s method (Nelson, 1976).

Afterwards, the first order derivative of master modal coordinates (∂ q ( t )_m/∂r_j, $\partial \dot{q} {(t)}_{m} / \partial r_{j}$ , and $\partial \ddot{q} {(t)}_{m} / \partial r_{j}$ ) are calculated from the reduced first order sensitivity equation (Equation (16)) by using Newmark-β method (Newmark, 1959). Differentiating equations (14) and (15) with respect to the design variable gets the first order time history response derivatives of global structure as

\frac{\partial x {(t)}^{p}}{\partial r_{j}} = \frac{\partial (Φ_{m}^{p} - Θ_{s})}{\partial r_{j}} q {(t)}_{m} + (Φ_{m}^{p} - Θ_{s}) \frac{\partial q {(t)}_{m}}{\partial r_{j}}

(18)

\frac{\partial \dot{x} {(t)}^{p}}{\partial r_{j}} = \frac{\partial (Φ_{m}^{p} - Θ_{s})}{\partial r_{j}} \dot{q} {(t)}_{m} + (Φ_{m}^{p} - Θ_{s}) \frac{\partial \dot{q} {(t)}_{m}}{\partial r_{j}}

(19)

\frac{\partial \ddot{x} {(t)}^{p}}{\partial r_{j}} = \frac{\partial (Φ_{m}^{p} - Θ_{s})}{\partial r_{j}} \ddot{q} {(t)}_{m} + (Φ_{m}^{p} - Θ_{s}) \frac{\partial \ddot{q} {(t)}_{m}}{\partial r_{j}}

(20)

The size of reduced first order sensitivity equation is greatly reduced by using a small number of master eigenvectors by the proposed substructuring method. The calculation of first order time history response derivatives only requires calculation of reduced first order sensitivity equation and the first order eigenpair derivatives of one substructure, which thus alleviates the computational burden.

Higher order time history response derivatives

The second and higher order time history response derivatives with respect to design variables are derived herein. Differentiating the reduced first order sensitivity equation with respect to the design variable r_k leads to the reduced second order sensitivity equation as

\begin{array}{l} M_{d} \frac{\partial^{2} \ddot{q} {(t)}_{m}}{\partial r_{j} \partial r_{k}} + C_{d} \frac{\partial^{2} \dot{q} {(t)}_{m}}{\partial r_{j} \partial r_{k}} + K_{d} \frac{\partial^{2} q {(t)}_{m}}{\partial r_{j} \partial r_{k}} = \frac{\partial^{2} f_{d} (t)}{\partial r_{j} \partial r_{k}} - \frac{\partial^{2} M_{d}}{\partial r_{j} \partial r_{k}} \ddot{q} {(t)}_{m} - \frac{\partial^{2} C_{d}}{\partial r_{j} \partial r_{k}} \dot{q} {(t)}_{m} - \frac{\partial^{2} K_{d}}{\partial r_{j} \partial r_{k}} q {(t)}_{m} \\ - \frac{\partial M_{d}}{\partial r_{j}} \frac{\partial \ddot{q} {(t)}_{m}}{\partial r_{k}} - \frac{\partial M_{d}}{\partial r_{k}} \frac{\partial \ddot{q} {(t)}_{m}}{\partial r_{j}} - \frac{\partial C_{d}}{\partial r_{j}} \frac{\partial \dot{q} {(t)}_{m}}{\partial r_{k}} - \frac{\partial C_{d}}{\partial r_{k}} \frac{\partial \dot{q} {(t)}_{m}}{\partial r_{j}} - \frac{\partial K_{d}}{\partial r_{j}} \frac{\partial q {(t)}_{m}}{\partial r_{k}} - \frac{\partial K_{d}}{\partial r_{k}} \frac{\partial q {(t)}_{m}}{\partial r_{j}} \end{array}

(21)

where

\frac{\partial^{2} M_{d}}{\partial r_{j} \partial r_{k}} = \frac{\partial^{2} Θ_{s}^{T}}{\partial r_{j} \partial r_{k}} M^{p} Θ_{s} + \frac{\partial Θ_{s}^{T}}{\partial r_{j}} M^{p} \frac{\partial Θ_{s}}{\partial r_{k}} + \frac{\partial Θ_{s}^{T}}{\partial r_{k}} M^{p} \frac{\partial Θ_{s}}{\partial r_{j}} + Θ_{s}^{T} M^{p} \frac{\partial^{2} Θ_{s}}{\partial r_{j} \partial r_{k}}

(21a)

\frac{\partial^{2} K_{d}}{\partial r_{j} \partial r_{k}} = \frac{\partial^{2} Λ_{m}^{p}}{\partial r_{j} \partial r_{k}} + \frac{\partial^{2} Z_{m}}{\partial r_{j} \partial r_{k}}, \frac{\partial^{2} C_{d}}{\partial r_{j} \partial r_{k}} = a_{1} \frac{\partial^{2} M_{d}}{\partial r_{j} \partial r_{k}} + a_{2} \frac{\partial^{2} K_{d}}{\partial r_{j} \partial r_{k}}, \frac{\partial^{2} f_{d} (t)}{\partial r_{j} \partial r_{k}} = (\frac{\partial^{2} {[Φ_{m}^{p}]}^{T}}{\partial r_{j} \partial r_{k}} - \frac{\partial^{2} Θ_{s}^{T}}{\partial r_{j} \partial r_{k}}) f {(t)}_{ext}^{p}

(21b)

In equations (21a)–(21b), the first order derivative of master eigenpairs, residual flexibility, and master modal coordinates is obtained when the first order time response derivatives are calculated. The second order derivatives (

\frac{\partial^{2} Θ_{s}}{\partial r_{j} \partial r_{k}}

、

\frac{\partial^{2} Z_{m}}{\partial r_{j} \partial r_{k}}

) are calculated as

\begin{array}{l} \frac{\partial^{2} Θ_{s}}{\partial r_{j} \partial r_{k}} = \frac{\partial^{2} F^{p}}{\partial r_{j} \partial r_{k}} D^{T} {(D F^{p} D^{T})}^{- 1} D Φ_{m}^{p} + \frac{\partial F^{p}}{\partial r_{j}} D^{T} \frac{\partial {(D F^{p} D^{T})}^{- 1}}{\partial r_{k}} D Φ_{m}^{p} + \frac{\partial F^{p}}{\partial r_{j}} D^{T} {(D F^{p} D^{T})}^{- 1} D \frac{\partial Φ_{m}^{p}}{\partial r_{k}} \\ + \frac{\partial F^{p}}{\partial r_{k}} D^{T} \frac{\partial {(D F^{p} D^{T})}^{- 1}}{\partial r_{j}} D Φ_{m}^{p} + F^{p} D^{T} \frac{\partial^{2} {(D F^{p} D^{T})}^{- 1}}{\partial r_{j} \partial r_{k}} D Φ_{m}^{p} + F^{p} D^{T} \frac{\partial {(D F^{p} D^{T})}^{- 1}}{\partial r_{j}} D \frac{\partial Φ_{m}^{p}}{\partial r_{k}} \\ + \frac{\partial F^{p}}{\partial r_{k}} D^{T} {(D F^{p} D^{T})}^{- 1} D \frac{\partial Φ_{m}^{p}}{\partial r_{j}} + F^{p} D^{T} \frac{\partial {(D F^{p} D^{T})}^{- 1}}{\partial r_{k}} D \frac{\partial Φ_{m}^{p}}{\partial r_{j}} + F^{p} D^{T} {(D F^{p} D^{T})}^{- 1} D \frac{\partial^{2} Φ_{m}^{p}}{\partial r_{j} \partial r_{k}} \end{array}

(22)

\begin{array}{l} \frac{\partial^{2} Z_{m}}{\partial r_{j} \partial r_{k}} = \frac{\partial^{2} {(Φ_{m}^{p})}^{T}}{\partial r_{j} \partial r_{k}} D^{T} {(D F^{p} D^{T})}^{- 1} D Φ_{m}^{p} + \frac{\partial {(Φ_{m}^{p})}^{T}}{\partial r_{j}} D^{T} \frac{\partial {(D F^{p} D^{T})}^{- 1}}{\partial r_{k}} D Φ_{m}^{p} + \frac{\partial {(Φ_{m}^{p})}^{T}}{\partial r_{j}} D^{T} {(D F^{p} D^{T})}^{- 1} D \frac{\partial Φ_{m}^{p}}{\partial r_{k}} \\ + \frac{\partial {(Φ_{m}^{p})}^{T}}{\partial r_{k}} D^{T} \frac{\partial {(D F^{p} D^{T})}^{- 1}}{\partial r_{j}} D Φ_{m}^{p} + {(Φ_{m}^{p})}^{T} D^{T} \frac{\partial^{2} {(D F^{p} D^{T})}^{- 1}}{\partial r_{j} \partial r_{k}} D Φ_{m}^{p} + {(Φ_{m}^{p})}^{T} D^{T} \frac{\partial {(D F^{p} D^{T})}^{- 1}}{\partial r_{j}} D \frac{\partial Φ_{m}^{p}}{\partial r_{k}} \\ + \frac{\partial {(Φ_{m}^{p})}^{T}}{\partial r_{k}} D^{T} {(D F^{p} D^{T})}^{- 1} D \frac{\partial Φ_{m}^{p}}{\partial r_{j}} + {(Φ_{m}^{p})}^{T} D^{T} \frac{\partial {(D F^{p} D^{T})}^{- 1}}{\partial r_{k}} D \frac{\partial Φ_{m}^{p}}{\partial r_{j}} + {(Φ_{m}^{p})}^{T} D^{T} {(D F^{p} D^{T})}^{- 1} D \frac{\partial^{2} Φ_{m}^{p}}{\partial r_{j} \partial r_{k}} \end{array}

(23)

and the second order derivative (

\frac{\partial^{2} {(D F^{p} D^{T})}^{- 1}}{\partial r_{j} \partial r_{k}}

) is expressed as

\begin{array}{l} \frac{\partial^{2} {(D F^{p} D^{T})}^{- 1}}{\partial r_{j} \partial r_{k}} = {(D F^{p} D^{T})}^{- 1} D \frac{\partial F^{p}}{\partial r_{k}} D^{T} {(D F^{p} D^{T})}^{- 1} D \frac{\partial F^{p}}{\partial r_{j}} D^{T} {(D F^{p} D^{T})}^{- 1} - {(D F^{p} D^{T})}^{- 1} D \frac{\partial^{2} F^{p}}{\partial r_{j} \partial r_{k}} D^{T} {(D F^{p} D^{T})}^{- 1} \\ + {(D F^{p} D^{T})}^{- 1} D \frac{\partial F^{p}}{\partial r_{j}} D^{T} {(D F^{p} D^{T})}^{- 1} D \frac{\partial F^{p}}{\partial r_{k}} D^{T} {(D F^{p} D^{T})}^{- 1} \end{array}

(24)

It can be seen from equations (22)–(24) that, the second order derivatives ( $\frac{\partial^{2} Θ_{s}}{\partial r_{j} \partial r_{k}}$ 、 $\frac{\partial^{2} Z_{m}}{\partial r_{j} \partial r_{k}}$ ) are calculated from the master eigenpairs, residual flexibility matrix, and their first and second order derivatives.

When the two design variables r_j and r_k belong to the Qth substructure, the second order derivative of master eigenpairs and residual flexibility matrix ( $\frac{\partial^{2} Λ_{m}^{p}}{\partial r_{j} \partial r_{k}}$ ; $\frac{\partial^{2} Φ_{m}^{p}}{\partial r_{j} \partial r_{k}}$ ; $\frac{\partial^{2} F^{p}}{\partial r_{j} \partial r_{k}}$ ) is calculated from those of Qth substructure only, which has the form of

\frac{\partial^{2} Λ_{m}^{p}}{\partial r_{j} \partial r_{k}} = [\begin{array}{c} 0 & 0 & 0 \\ 0 & \frac{\partial^{2} Λ_{m}^{(Q)}}{\partial r_{j} \partial r_{k}} & 0 \\ 0 & 0 & 0 \end{array}], \frac{\partial^{2} Φ_{m}^{p}}{\partial r_{j} \partial r_{k}} = [\begin{array}{c} 0 & 0 & 0 \\ 0 & \frac{\partial^{2} Φ_{m}^{(Q)}}{\partial r_{j} \partial r_{k}} & 0 \\ 0 & 0 & 0 \end{array}]

(25)

\frac{\partial^{2} F_{m}^{p}}{\partial r_{j} \partial r_{k}} = [\begin{array}{c} 0 & 0 & 0 \\ 0 & \frac{\partial^{2} ({(K^{(Q)})}^{- 1} - Φ_{m}^{(Q)} {(Λ_{m}^{(Q)})}^{- 1} {[Φ_{m}^{(Q)}]}^{T})}{\partial r_{j} \partial r_{k}} & 0 \\ 0 & 0 & 0 \end{array}]

(26)

The first and second order derivatives of master eigenpairs of Qth substructure are required to form the coefficients of the reduced second order sensitivity matrix (equation (21)). The second order derivative of master eigenpairs of Qth substructure can be derived by differentiating the eigenequation a second time (Nelson, 1976).

When the the two design variables r_j and r_k belong to the two different substructures, $\frac{\partial^{2} Λ_{m}^{p}}{\partial r_{j} \partial r_{k}}$ , $\frac{\partial^{2} Φ_{m}^{p}}{\partial r_{j} \partial r_{k}}$ , and $\frac{\partial^{2} (Φ_{m}^{p} {(Λ_{m}^{p})}^{- 1} {[Φ_{m}^{p}]}^{T})}{\partial r_{j} \partial r_{k}}$ in equations (25) and (26) are zeroes. Only the first order derivative of master eigenpairs of two substructures that contain the design variables is required for calculating the coefficients of reduced second order sensitivity equation, and second order derivative of master eigenpairs is zero. In addition, the first order derivative of assembled master eigenpairs ( $\frac{\partial Λ_{m}^{p}}{\partial r_{j}}$ 、 $\frac{\partial Φ_{m}^{p}}{\partial r_{j}}$ 、 $\frac{\partial Λ_{m}^{p}}{\partial r_{k}}$ 、 $\frac{\partial Φ_{m}^{p}}{\partial r_{k}}$ ) is diagonally assembled by those of substructures (equation (17)). Since the two design variables belong to two different substructures, the multiplication of the first order derivative of assembled master eigenpairs with respect to different variable is zero.

In equation (21), the first order derivative of master modal coordinates (∂ q ( t )_m/∂r_j、 $\partial \dot{q} {(t)}_{m} / \partial r_{j}$ 、 $\partial \ddot{q} {(t)}_{m} / \partial r_{j}$ 、 ∂ q ( t )_m/∂r_k 、 $\partial \dot{q} {(t)}_{m} / \partial r_{k}$ 、 $\partial \ddot{q} {(t)}_{m} / \partial r_{k}$ ) has been calculated in Section 3. The second order derivative of master modal coordinates ( $\frac{\partial^{2} \ddot{q} {(t)}_{m}}{\partial r_{j} \partial r_{k}}$ 、 $\frac{\partial^{2} \dot{q} {(t)}_{m}}{\partial r_{j} \partial r_{k}}$ 、 $\frac{\partial^{2} q {(t)}_{m}}{\partial r_{j} \partial r_{k}}$ ) can be calculated from the reduced second order sensitivity equation (Equation (21)) by Newmark-β method. Differentiating the equations (18)–(20) with respect to the second design variable gets the second order time history response derivatives as

\frac{\partial^{2} x {(t)}^{p}}{\partial r_{j} \partial r_{k}} = (\frac{\partial^{2} Φ_{m}^{p}}{\partial r_{j} \partial r_{k}} - \frac{\partial^{2} Θ_{s}}{\partial r_{j} \partial r_{k}}) q {(t)}_{m} + (\frac{\partial Φ_{m}^{p}}{\partial r_{j}} - \frac{\partial Θ_{s}}{\partial r_{j}}) \frac{\partial q {(t)}_{m}}{\partial r_{k}} + (\frac{\partial Φ_{m}^{p}}{\partial r_{k}} - \frac{\partial Θ_{s}}{\partial r_{k}}) \frac{\partial q {(t)}_{m}}{\partial r_{j}} + (Φ_{m}^{p} - Θ_{s}) \frac{\partial^{2} q {(t)}_{m}}{\partial r_{j} \partial r_{k}}

(27)

\frac{\partial^{2} \dot{x} {(t)}^{p}}{\partial r_{j} \partial r_{k}} = (\frac{\partial^{2} Φ_{m}^{p}}{\partial r_{j} \partial r_{k}} - \frac{\partial^{2} Θ_{s}}{\partial r_{j} \partial r_{k}}) \dot{q} {(t)}_{m} + (\frac{\partial Φ_{m}^{p}}{\partial r_{j}} - \frac{\partial Θ_{s}}{\partial r_{j}}) \frac{\partial \dot{q} {(t)}_{m}}{\partial r_{k}} + (\frac{\partial Φ_{m}^{p}}{\partial r_{k}} - \frac{\partial Θ_{s}}{\partial r_{k}}) \frac{\partial \dot{q} {(t)}_{m}}{\partial r_{j}} + (Φ_{m}^{p} - Θ_{s}) \frac{\partial^{2} \dot{q} {(t)}_{m}}{\partial r_{j} \partial r_{k}}

(28)

\frac{\partial^{2} \ddot{x} {(t)}^{p}}{\partial r_{j} \partial r_{k}} = (\frac{\partial^{2} Φ_{m}^{p}}{\partial r_{j} \partial r_{k}} - \frac{\partial^{2} Θ_{s}}{\partial r_{j} \partial r_{k}}) \ddot{q} {(t)}_{m} + (\frac{\partial Φ_{m}^{p}}{\partial r_{j}} - \frac{\partial Θ_{s}}{\partial r_{j}}) \frac{\partial \ddot{q} {(t)}_{m}}{\partial r_{k}} + (\frac{\partial Φ_{m}^{p}}{\partial r_{k}} - \frac{\partial Θ_{s}}{\partial r_{k}}) \frac{\partial \ddot{q} {(t)}_{m}}{\partial r_{j}} + (Φ_{m}^{p} - Θ_{s}) \frac{\partial^{2} \ddot{q} {(t)}_{m}}{\partial r_{j} \partial r_{k}}

(29)

The reduced higher order sensitivity equation is derived by differentiating the reduced motion equation with respect to k design variables as

\begin{array}{l} M_{d} \frac{\partial^{k} \ddot{q} {(t)}_{m}}{\partial r_{1} \partial r_{2} \dots \partial r_{k}} + C_{d} \frac{\partial^{k} \dot{q} {(t)}_{m}}{\partial r_{1} \partial r_{2} \dots \partial r_{k}} + K_{d} \frac{\partial^{k} q {(t)}_{m}}{\partial r_{1} \partial r_{2} \dots \partial r_{k}} = \frac{\partial^{k} f_{d} (t)}{\partial r_{1} \partial r_{2} \dots \partial r_{k}} \\ - (\frac{\partial^{k} M_{d}}{\partial r_{1} \partial r_{2} \dots \partial r_{k}} \ddot{q} {(t)}_{m} + \frac{\partial^{k - 1} M_{d}}{\partial r_{1} \partial r_{2} \dots \partial r_{k - 1}} \frac{\partial \ddot{q} {(t)}_{m}}{\partial r_{k}} + \dots + \frac{\partial M_{d}}{\partial r_{1}} \frac{\partial^{k - 1} \ddot{q} {(t)}_{m}}{\partial r_{2} \dots \partial r_{k - 1} \partial r_{k}}) \\ - (\frac{\partial^{k} C_{d}}{\partial r_{1} \partial r_{2} \dots \partial r_{k}} \dot{q} {(t)}_{m} + \frac{\partial^{k - 1} C_{d}}{\partial r_{1} \partial r_{2} \dots \partial r_{k - 1}} \frac{\partial \dot{q} {(t)}_{m}}{\partial r_{k}} + \dots + \frac{\partial C_{d}}{\partial r_{1}} \frac{\partial^{k - 1} \dot{q} {(t)}_{m}}{\partial r_{2} \dots \partial r_{k - 1} \partial r_{k}}) \\ - (\frac{\partial^{k} K_{d}}{\partial r_{1} \partial r_{2} \dots \partial r_{k}} q {(t)}_{m} + \frac{\partial^{k - 1} K_{d}}{\partial r_{1} \partial r_{2} \dots \partial r_{k - 1}} \frac{\partial q {(t)}_{m}}{\partial r_{k}} + \dots + \frac{\partial K_{d}}{\partial r_{1}} \frac{\partial^{k - 1} q {(t)}_{m}}{\partial r_{2} \dots \partial r_{k - 1} \partial r_{k}}) \end{array}

(30)

When the k design variables belong to the same substructure, only kth and lower order derivatives of master eigenpairs of one substructure are required to form the reduced kth order sensitivity equation and derivative matrices of other substructures are zeroes. When the k design variables belong to the different substructures, kth order derivatives are zeroes and only (k-1)th and lower order derivatives of substructural master eigenpairs of k substructures are required for calculation of kth order time history response derivatives. Differentiating the equations (14) and (15) with respect to the k design variable gets the kth order time history response derivatives as

\frac{\partial^{k} x {(t)}^{p}}{\partial r_{1} \partial r_{2} \dots \partial r_{k}} = \frac{\partial^{k} (Φ_{m}^{p} - Θ_{s})}{\partial r_{1} \partial r_{2} \dots \partial r_{k}} q {(t)}_{m} + \frac{\partial^{k - 1} (Φ_{m}^{p} - Θ_{s})}{\partial r_{1} \partial r_{2} \dots \partial r_{k - 1}} \frac{\partial q {(t)}_{m}}{\partial r_{k}} + \dots + \frac{\partial (Φ_{m}^{p} - Θ_{s})}{\partial r_{1}} \frac{\partial^{k - 1} q {(t)}_{m}}{\partial r_{2} \dots \partial r_{k - 1} \partial r_{k}}

(31)

\frac{\partial^{k} \dot{x} {(t)}^{p}}{\partial r_{1} \partial r_{2} \dots \partial r_{k}} = \frac{\partial^{k} (Φ_{m}^{p} - Θ_{s})}{\partial r_{1} \partial r_{2} \dots \partial r_{k}} \dot{q} {(t)}_{m} + \frac{\partial^{k - 1} (Φ_{m}^{p} - Θ_{s})}{\partial r_{1} \partial r_{2} \dots \partial r_{k - 1}} \frac{\partial \dot{q} {(t)}_{m}}{\partial r_{k}} + \dots + \frac{\partial (Φ_{m}^{p} - Θ_{s})}{\partial r_{1}} \frac{\partial^{k - 1} \dot{q} {(t)}_{m}}{\partial r_{2} \dots \partial r_{k - 1} \partial r_{k}}

(32)

\frac{\partial^{k} \ddot{x} {(t)}^{p}}{\partial r_{1} \partial r_{2} \dots \partial r_{k}} = \frac{\partial^{k} (Φ_{m}^{p} - Θ_{s})}{\partial r_{1} \partial r_{2} \dots \partial r_{k}} \ddot{q} {(t)}_{m} + \frac{\partial^{k - 1} (Φ_{m}^{p} - Θ_{s})}{\partial r_{1} \partial r_{2} \dots \partial r_{k - 1}} \frac{\partial \ddot{q} {(t)}_{m}}{\partial r_{k}} + \dots + \frac{\partial (Φ_{m}^{p} - Θ_{s})}{\partial r_{1}} \frac{\partial^{k - 1} \ddot{q} {(t)}_{m}}{\partial r_{2} \dots \partial r_{k - 1} \partial r_{k}}

(33)

In a word, the proposed substructuring method calculates the first, second, and higher order time history response derivatives from the reduced sensitivity equation and eigensensitivity of a few substructures that contain the design variables. Compared with the traditional global method, the proposed substructuring method avoids calculating higher order derivative matrices of the global structure and reduces the size of higher order sensitivity equation by using a few master eigenvectors, which would accelerate the computational efficiency. Compared with the traditional substructuring method considering the elastic effect of slave eigenvectors only, the proposed substructuring method considers the inertial effect of slave eigenvectors, which also ensures the accuracy of time history response derivatives with a small number of master eigenvectors meanwhile.

Numerical simulation 1: a one-bay plane frame

A one-bay steel plane frame is used to test the accuracy of the proposed substructuring method to calculate the time history response derivatives. As Figure 1 shows, the finite element model of frame has 30 elements, 31 nodes, and 87 dofs. The length of each element is 100 mm. The cross sectional area of the column and beam is 50.50 × 6.0 mm² and 40.5 × 6.0 mm², respectively. The mass density is 7.76 × 10³ kg/m³. The damping constants for Rayleigh damping are a₁ = 0.2932 s⁻¹, and a₂ = 0.0055 s. The Young’s modulus of the steel frame is 2.10 × 10¹¹ N·m². The frame is subject to the earthquake excitation in horizontal direction, as shown in Figure 2. The global structure is divided into two same substructures with 15 elements, 16 nodes, and 45 dofs, respectively.

Figure 1.

The finite element model of frame.

Figure 2.

Earthquake excitation.

The global method and traditional substructuring method are used for comparison. The global method solves the time history response derivatives by treating the structure as a whole. The results calculated by global method are regarded as the exact ones. The traditional substructuring method retains two or 10 master eigenvectors in each substructure to calculate the time history response derivatives, where only the elastic effect of slave eigenvectors is considered [39]. The proposed substructuring method retains 2 master eigenvectors in each substructure and both the elastic and inertial effects are considered. For brevity, only the second order time history response derivatives are demonstrated herein. r₁₂ and r₁₆ represent the bending rigidity factor of Element 12 and Element 16, respectively. r₁₂ and r₁₆ belong to the Substructure one and Substructure 2, respectively.

Figure 3 compares the global and substructuring methods in calculation of the second order time history response derivatives of Node 15 in the horizontal direction with respect to the design variable pairs (r₁₂, r₁₂) located in the same substructure, where only the first and second order master eigenvector derivatives of Substructure one are required to recover the second order time history response derivatives. Figure 4 compares the global and substructuring methods in calculation of the second order time history response derivatives of Node 15 in the horizontal direction with respect to the design variable pairs (r₁₂, r₁₆) located in the different substructures, where only the first order master eigenvector derivatives of the two substructures are required. It can be seen from Figures 3 and 4 that, when merely 2 master eigenvectors are retained in each substructure, the curves of the proposed substructuring method overlap with those of the global method while the curves of the traditional substructuring method shows a noticeable shift from the reference curves. When the number of master eigenvectors reaches up to 10, the curves of the traditional substructuring method are in good agreement with the reference curves.

Figure 3.

The second order response derivatives with respect to the two design variables located in the same substructure: (a) Second order displacement derivative; (b) Second order velocity derivative; (c) Second order acceleration derivative.

Figure 4.

The second order response derivatives with respect to the two design variables located in the different substructures: (a) Second order displacement derivative; (b) Second order velocity derivative; (c) Second order acceleration derivative.

To quantify the difference between the substructuring and global methods to calculate the second order time history response derivatives, the relative error is defined as

Relative error (\frac{\partial^{2} x_{G}}{\partial r_{j} \partial r_{k}}, \frac{\partial^{2} x_{Sub}}{\partial r_{j} \partial r_{k}}) = \frac{{‖ \frac{\partial^{2} x_{G}}{\partial r_{j} \partial r_{k}} - \frac{\partial^{2} x_{Sub}}{\partial r_{j} \partial r_{k}} ‖}_{2}}{{‖ \frac{\partial^{2} x_{G}}{\partial r_{j} \partial r_{k}} ‖}_{2}}

(34)

where

\frac{\partial^{2} x_{G}}{\partial r_{k} \partial r_{j}}

and

\frac{\partial^{2} x_{Sub}}{\partial r_{k} \partial r_{j}}

represent the second order displacement derivative by global and substructuring methods, respectively. Figure 5 shows the relative error values of second order displacement derivatives of all nodes in the horizontal direction with respect to the design variable pairs (r₁₂, r₁₂) and (r₁₂, r₁₆) by the proposed and traditional substructuring methods, respectively. Figure 5 shows that the accuracy of the proposed substructuring method is obviously higher than the traditional substructuring method when each substructure retains 2 master eigenvectors. When the number of master eigenvectors increases up to 10, the traditional substructuring method reaches the comparable accuracy in comparison with the proposed substructuring method with 2 master eigenvectors in each substructure. It means that the traditional substructuring method requires more master eigenvectors to get accurate results if the inertial effect of slave eigenvectors is not considered.

Figure 5.

The relative error for second order displacement derivatives by proposed and traditional substructuring methods: (a) the two design variables located in the same substructure; (b) the two design variables located in the different substructures.

To investigate the influence of master eigenvectors on the accuracy of the proposed substructuring method, the number of master eigenvectors in each substructure is chosen as 2, 5, and 10 for comparison. Figure 6 shows the relative error values of proposed substructuring method to calculate the second order displacement derivatives of all nodes in horizontal direction with respect to the design variable pairs (r₁₂, r₁₂) and (r₁₂, r₁₆), respectively. Table 1 lists the maximum and mean of relative errors shown in Figure 6. As shown in Figure 6 and Table 1, the maximum and mean of relative errors of second order displacement derivatives with respect to the two design variables located in one substructure are 6.80 × 10⁻³ and 1.50 × 10⁻³, respectively with each substructure retaining 2 master eigenvectors, which are reduced to 1.37 × 10⁻⁴ and 6.88 × 10⁻⁶ when 10 master eigenvectors in each substructure are retained. The maximum of relative error values of the second order displacement derivatives with respect to the two design variables located in the different substructures are reduced from 4.70 × 10⁻³ to 9.52 × 10⁻⁵ when the number of master eigenvectors increases from 2 to 10. It can be concluded that the proposed substructuring method is accurate to calculate the second order time history response derivatives since the contribution of discarded slave eigenvectors are accurately considered and the accuracy would increase with the increase of the number of master eigenvectors.

Figure 6.

The relative error for second order displacement derivatives by proposed substructuring method with different number of master eigenvectors: (a) the two design variables located in the same substructure; (b) the two design variables located in the different substructures.

Table 1.

The relative error of proposed substructuring method to second order displacement derivatives.

	In the same substructure			In the different substructures
No. of master modes	2	5	10	2	5	10
Maximum	6.80 × 10⁻³	2.89 × 10⁻⁴	1.37 × 10⁻⁴	4.70 × 10⁻³	3.17 × 10⁻⁴	9.52 × 10⁻⁵
Mean	1.50 × 10⁻³	3.29 × 10⁻⁵	6.88 × 10⁻⁶	1.60 × 10⁻³	4.42 × 10⁻⁵	8.02 × 10⁻⁶

Numerical simulation 2: A highway bridge

Application to a highway bridge is used to verify the time efficiency of the proposed substructuring method. The highway bridge is made up of a nine-span box girder and 10 circular columns. The sectional inertial moment of girder is 1.86 m⁴ in X direction and 34.04 m⁴ in Y direction. The sectional inertial moment of column is 0.64 m⁴ in both X and Y directions. The finite element model of highway bridge is shown in Figure 7, which has 750 elements, 751 nodes, and 2223 dofs. The damping constants for Rayleigh damping are a₁ = 0.0063 s⁻¹, and a₂ = 0.3940 s. The bridge is excited by the earthquake force in horizontal direction and the earthquake acceleration is shown in Figure 2. The bridge is divided into three substructures as Figure 7. As before, the time history response derivatives of the global method are regarded as the exact results. The second order displacement derivatives are presented for brevity in this numerical case. A total of 30 master eigenvectors is retained in each substructure for the calculation of the second order displacement derivatives by using the proposed substructuring method. The traditional substructuring method retains 30 or 60 master eigenvectors in each substructure.

Figure 7.

The finite element model of a highway bridge.

The elemental bending rigidity factor is chosen as design variable. Figure 8 shows the second order displacement derivative of a randomly selected node in horizontal direction with respect to the design variable pairs (r₁₂, r₁₃) and (r₁₂, r₁₅₆) using the substructuring and global methods. r₁₂, r₁₃, and r₁₅₆ stand for the bending rigidity factor of Elements 12, 13, and 156, respectively. Elements 12 and 13 belong to the first span box girder of Substructure 1. Element 156 belongs to the sixth span box girder, which belongs to Substructure 2. The design variable pairs (r₁₂, r₁₃) belong to the same substructure and the design variable pairs (r₁₂, r₁₅₆) belong to the different substructures. It can be seen from Figure 8 that, the results calculated by the proposed substructuring method overlap with those of the global method, which proves that the second order displacement derivative calculated by proposed substructuring method has high accuracy. This is because the proposed substructuring method considers the elastic and inertial effects of slave eigenvectors, which ensures the accuracy of the proposed substructuring method with a small number of master eigenvectors. However, the curves of the traditional substructuring method overlap with the exact curves only when the number of master eigenvectors increases up to 60. When the inertial effect of slave eigenvectors is not considered, the traditional substructuring method has to retain more master eigenvectors to get accurate results.

Figure 8.

The second order displacement derivatives with respect to the two design variables: (a) located in the same substructure; (b) located in the different substructures.

Table 2 compares the computational time in calculation of the second order displacement derivatives with respect to the design variable pairs (r_j, r_k) by using the substructuring and global methods. r_j stands for the bending rigidity factor of a total of 90 box girder elements located in Substructure 1. r_k stands for the bending rigidity factor of a total of 270 box girder elements of the global structure. The computational time of calculating the second order displacement derivatives involves calculating the displacement and the first order displacement derivatives by both the aforementioned methods. The global method consumes 41 626.20 s to calculate the second order displacement derivatives with respect to 90 design variables of r_j and 270 design variables of r_k. When each substructure retains 30 master eigenvectors, the traditional substructuring method takes up 7.39% of the computational time by the global method to calculate the second order displacement derivatives, which is slightly less than the proposed substructuring method. This is because the proposed substructuring method requires a little extra time to consider the inertial effect of slave eigenvectors. To reach the comparable accuracy as the proposed substructuring method, the traditional substructuring method has to retain as much as 60 master eigenvectors and the time consumed by the traditional substructuring method accounts for 11.80% of that of the global method.

Table 2.

Comparison of computational time by the global and substructuring methods to second order displacement derivatives.

	Global method	Traditional substructuring method (30 master eigenvectors)			Traditional substructuring method (60 master eigenvectors)			Proposed substructuring method (30 master eigenvectors)
		Location of design parameter r_k			Location of design parameter r_k			Location of design parameter r_k
		Sub 1	Sub 2	Sub 3	Sub 1	Sub 2	Sub 3	Sub 1	Sub 2	Sub 3
Number of design parameter	270	90	90	90	90	90	90	90	90	90
Time (second)		1190.41	346.93	346.93	1870.08	585.95	585.95	1486.67	602.44	602.44
Total time (second)	41 626.20		3074.68			4912.06			4178.22
Ratio	100.00%		7.39%			11.80%			10.04%

The computational time consumed by the proposed substructuring method accounts for only 10.04% of that of global method. When the two design variables r_j and r_k belong to Substructure 1, the proposed substructuring method consumes 1486.67 s in calculation of the second order displacement derivatives with respect to 90 design variables of r_j and 90 design variables of r_k. When the two design variables r_j and r_k belong to Substructure 1 and 2, respectively, the proposed substructuring method takes 602.44 s to calculate the second order displacement derivatives with respect to the 90 design variables of r_j and 90 design variables of r_k. For the proposed substructuring method, calculating the second order response derivatives with respect to the two design variables located in the same substructure consumes more time than the two design variables located in the different substructures. This is because the former requires calculating both the first and second order derivatives of master eigenvectors whereas the latter only requires calculating the first order derivative of master eigenvectors.

The proposed substructuring method assembles a few master eigenvectors and the derivative of master eigenvectors of each substructure to recover the times history response and response derivatives of the global structure. The size of reduced motion and sensitivity equations is equal to the number of master eigenvectors, which is reduced from 2223 × 2223 to 90 × 90. Only the derivative matrices of substructure containing the design variables are calculated and those of other substructures are not required. In addition, the proposed substructuring method increases the accuracy by considering the inertial effect of slave eigenvectors rather than by retaining more master eigenvectors compared with the traditional substructuring method. Accordingly, the proposed substructuring method can calculate the time history response derivatives of global structure efficiently.

Conclusions

This paper proposes a substructuring method for efficient calculation of the first, second, and higher order time history response derivatives of large-scale structures. The large-scale global structure is disassembled into small independent substructures. The derivative matrices of a few master eigenvectors with respect to the design variables of each independent substructure are then assembled to form the reduced first, second, and higher order sensitivity equations. The first, second, and higher order time history response derivatives of the global structure are efficiently solved from the reduced sensitivity equations by using Newmark-β method. As the derivative matrices are zeroes for substructures that do not include the design variables, the time history response derivatives of global structure are determined from the derivative matrices of a few substructures that contain the design variables. The size of sensitivity equations is reduced to the number of master eigenvectors. The equivalent eigenvector is derived to relate the substructural master and slave eigenvectors. The inertial effect of substructural slave eigenvectors is compensated by the equivalent eigenvector related term. Accordingly, the proposed substructuring method is efficient and accurate to calculate the time history response derivatives.

The proposed substructuring method is applied to calculate the second order time history response derivatives of a numerical one-bay plane frame and a numerical highway bridge. The results by proposed substructuring method show high accuracy with a small number of master eigenvectors and the computational time is far less than the global method. While the traditional substructuring method, which only considers the elastic effect of slave eigenvectors, has to retain more master eigenvectors to achieve the comparable accuracy. The proposed substructuring method allows for rapid and accurate calculation of the time history response derivatives of large-scale structure. Besides, the proposed substructuring method is promising to alleviate the computational load of the time history response based finite element model updating since the time history response derivatives with respect to all updating variables need to be calculated iteratively.

Footnotes

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported by the National Natural Science Foundation of China (NSFC, contract number: 52078395, 51922046, 51778258), Wuhan Institute of Technology Research Start-Up Fund (23QD73), Natural Science Foundation of Hubei province for Distinguished Young Scholars (2023AFA103), Young Topnotch Talent Cultivation Program of Hubei Province, and Hubei Provincial Engineering Research Centre for Green Civil Engineering Materials and Structures.

ORCID iDs

Jiajing Li

Qiaoyun Wu

References

Sheikhi Azqandi

Hassanzadeh

(2021) First- and second-order sensitivity analysis of finite element models using extended complex variables method. Archive of Applied Mechanics 91: 4263–4277.

Chen

Sun

(2021) A component mode synthesis method for reduced-order modeling of cable networks in cable-stayed bridges. Journal of Sound and Vibration 491: 115769.

Chung

Kwon

Park

, et al. (2009) Sensitivity analysis by the adjoint variable method for optimization of the die compaction process in particulate materials processing. Finite Elements in Analysis and Design 45(11): 836–844.

Chung

Kim

Chae

, et al. (2021) An iterative scheme of flexibility‐based component mode synthesis with higher‐order residual modal compensation. International Journal for Numerical Methods in Engineering 122(13): 3171–3190.

Craig

Bampton

MCC

(1968) Coupling of substructures for dynamic analyses. AIAA Journal 6(7): 1313–1319.

De Klerk

Rixen

Voormeeren

(2008) General framework for dynamic substructuring: history, review, and classification of techniques. AIAA Journal 46(5): 1169–1181.

Lim

Kim

, et al. (2020) A family of Craig–Bampton methods considering residual mode compensation. Applied Mathematics and Computation 369: 124822.

Gomez-Farias

Montoya

Millwater

(2015) Complex finite element sensitivity method for creep analysis. International Journal of Pressure Vessels and Piping 133: 27–42.

Gruber

Gille

Rixen

(2019) A strategy to stabilize the transient analysis and increase the approximation accuracy of dual Craig-Bampton reduced systems. Finite Elements in Analysis and Design 160: 32–45.

10.

Hou

Xia

(2021) Review on the new development of vibration-based damage identification for civil engineering structures: 2010–2019. Journal of Sound and Vibration 491: 115741.

11.

Hou

Jankowski

(2015) Frequency-domain substructure isolation for local damage identification. Advances in Structural Engineering 18(1): 137–153.

12.

Hurty

(1965) Dynamic analysis of structural systems using component modes. AIAA Journal 3(5): 678–685.

13.

Ilinca

Pelletier

Borggaard

(2007) A continuous second-order sensitivity equation method for time-dependent incompressible laminar flows. International Journal for Numerical Methods in Fluids 55(6): 565–587.

14.

Jalali

Rideout

(2022) Substructural damage detection using frequency response function based inverse dynamic substructuring. Mechanical Systems and Signal Processing 163: 108166.

15.

Kim

Lee

(2015) An enhanced Craig-Bampton method. International Journal for Numerical Methods in Engineering 103: 79–93.

16.

Kim

Lee

(2017) Improving the accuracy of the dual Craig-Bampton method. Computers & Structures 191: 22–32.

17.

Koh

See

Balendra

(1991) Estimation of structural parameters in time domain: a substructure approach. Earthquake Engineering & Structural Dynamics 20(8): 787–801.

18.

Kron

(1963) Diakoptics-Piecewise Solution of Large-Scale Systems. London: Macdonald and Co.

19.

Lei

Lai

, et al. (2015) Identification of flexible buildings with bending deformation and the unmeasured earthquake ground motion. Science China Technological Sciences 58(3): 454–461.

20.

Lei

Tian

, et al. (2019) Direct method for second-order sensitivity analysis of modal strain energy. Journal of Sound and Vibration 462: 114926.

21.

Law

(2012) Substructural damage detection with incomplete information of the structure. Journal of Applied Mechanics 79(4): 041003.

22.

Liew

Leek

Liu

(1996) Mixed-interface substructures for dynamic analysis of flexible multibody systems. Engineering Structures 18(7): 495–503.

23.

Lin

Mottershead

(2020) A state-of-the-art review on theory and engineering applications of eigenvalue and eigenvector derivatives. Mechanical Systems and Signal Processing 138: 106536.

24.

MacNeal

(1971) A hybrid method of component mode synthesis. Computers & Structures 1(4): 581–601.

25.

Mottershead

Link

Friswell

(2011) The sensitivity method infinite element model updating: a tutorial. Mechanical Systems and Signal Processing 25(7): 2275–2296.

26.

Nelson

(1976) Simplified calculation of eigenvector derivatives. AIAA Journal 14(14): 1201–1205.

27.

Newmark

(1959) A method of computation for structural dynamics. Journal of the Engineering Mechanics Division 85(1): 67–94.

28.

Park

(2004) Partitioned component mode synthesis via a flexibility approach. AIAA Journal 42(6): 1236–1245.

29.

Pradeepkumar

Nagaraj

(2022) Dynamic substructuring method for vibration analysis of complex structures. Journal of Vibration Engineering & Technologies 10: 313–333.

30.

Rubin

(1975) Improved component-mode representation for structural dynamic analysis. AIAA Journal 13(8): 995–1006.

31.

Tian

Weng

Xia

(2021) Kron’s substructuring method to the calculation of structural responses and response sensitivities of nonlinear systems. Journal of Sound and Vibration 502: 116101.

32.

Voormeeren

Rixen

(2012) A family of substructure decoupling techniques based on a dual assembly approach. Mechanical Systems and Signal Processing 27: 379–396.

33.

Wang

Lim

, et al. (2020) A strategy for decoupling of nonlinear systems using the inverse sub-structuring method and the parametric modal identification technique. Mechanical Systems and Signal Processing 140: 106695.

34.

Weng

Xia

, et al. (2011) Substructure based approach to finite element model updating. Computers & Structures 89(10): 772–782.

35.

Weng

Zhu

Xia

, et al. (2013) Substructuring approach to the calculation of higher-order eigensensitivity. Computers & Structures 117(2): 23–33.

36.

Weng

Zhu

Xia

, et al. (2013) Damage detection using the eigenparameter decomposition of substructural flexibility matrix. Mechanical Systems and Signal Processing 34(2): 19–38.

37.

Weng

Tian

Zhu

, et al. (2017) Dynamic condensation approach to calculation of structural responses and response sensitivities. Mechanical Systems and Signal Processing 88: 302–317.

38.

Weng

Zhu

Xia

, et al. (2020) A review on dynamic substructuring methods for model updating and damage detection of large-scale structures. Advances in Structural Engineering 23(3): 584–600.

39.

Yuen

Huang

(2018) Real-time substructural identification by boundary force modeling. Structural Control and Health Monitoring 25(5): e2151.

40.

Zhang

(2020) Taylor approximation of second degree for complex mode in asymmetrical damped system. AIAA Journal 58(2): 929–938.

41.

Zhang

Yang

Wang

, et al. (2018) Improving substructure identification accuracy of shear structures using virtual control system. Smart Materials and Structures 27(2): 025013.

42.

Zhang

(2023) High-order sensitivity analysis of complex modal parameters and their comparison. Archive of Applied Mechanics 93: 1123–1138.

43.

Zhu

Mao

Weng

(2013) Calculation of dynamic response sensitivity to substructural damage identification under moving load. Advances in Structural Engineering 16(9): 1621–1631.

44.

Zhu

Weng

, et al. (2019) Calculation of structural response and response sensitivity with improved substructuring method. Journal of Aerospace Engineering 32(3): 1–17.

45.

Zhu

Tian

, et al. (2021) An enhanced substructure-based response sensitivity method for finite element model updating of large-scale structures. Mechanical Systems and Signal Processing 154: 107359.