Observation vector reconstruction-based nonparametric nonlinear restoring force identification for granules-structures coupled vibrating system

Abstract

The nonlinear restoring force (NRF) generated by the collision and friction between particles and structures is the leading cause of the complex dynamic response of the granules-structures coupled vibrating system (GSCVS). Identification of NRF can provide critical information for post-event damage diagnosis and structural design of immersed structures. However, the spatial distribution and dynamic response of the particles near the structures are diverse and complex, making it difficult to describe the NRF with an accurate mathematical model. This paper proposed a data-based nonparametric method to estimate the NRF in the GSCVS. A nonparametric model of NRF that considered the additional effects of particles on both sides of the structures and consisted of system response and undetermined coefficients was developed. The observation vector of the conventional Extended Kalman Filter (EKF) was reconstructed by the sparse measurement of the strain response. The reconstructed observation vector contains three response components: translational displacement, translational acceleration, and rotational acceleration, in which the rotational acceleration response is difficult to measure in engineering applications. The proposed EKF based on observation vector reconstruction (EKF-OVR) can identify the undetermined coefficients in the nonparametric model, and then the NRF can be calculated. Numerical studies showed that EKF-OVR achieved higher accuracy and noise robustness than the conventional EKF and the data fusion based EKF. A dynamic experimental study on granules-beam coupled vibrating system (GBCVS) was conducted, and the proposed algorithm was employed to identify the NRF of the GBCVS. The effects of excitation amplitude, particle size, and immersed depth on NRF are analyzed, and it is found that higher harmonic components in the NRF led to period doubling and chaos of the beam response.

1. Introduction

The granules-structures coupled vibrating system (GSCVS) is a typical nonlinear system in industrial and agricultural production, such as coal vibration screening, concrete mixing preparation, crops processing, storage and transportation, particle damper, excavator operation, ballasted railway track bed design, mine detection, and desert vehicles (Chandravanshi and Mukhopadhyay, 2017; Chung and Wu, 2019; Jahani et al., 2015; Kuriyama and Saeki, 2022; Matthews et al., 1998; Xiong et al., 2017). At present, some important progress has been made in studying the dynamic characteristics of GSCVS. Lee (1998) studied the motion of strongly inelastic particles in a narrow vibrating tube using molecular-dynamics simulation, and found the center of mass of the particle cloud can be described by a superposition of modes of different frequencies. Donskoy et al. (2005) presented a theoretical model and detection algorithm considering quadratic and cubic nonlinearity of soil-mine interface. It was found that the cubic nonlinearity could be a significant contributor to the nonlinear response. Generally, such cubic nonlinear systems, called Duffing-type systems, will lead to a chaotic response under certain conditions (Li et al., 2023). Kang et al. (2007) conducted experiments to examine how the resonance frequency of a vibrating plate varied with different particle sizes of dry sand, and discovered a point of minimum frequency. Rattanadit et al. (2009) developed a coupled discrete element method-finite element method (DEM-FEM) computational model to study the interaction between a granular layer and an elastic foundation. The model was used to simulate the quasi-static bending of the granular layer and observed the changes in the force chain under two scenarios: with or without rolling resistance. An experiment involving a sinusoidally excited beam immersed in glass bead particles contained in a box was conducted, and the period-doubling bifurcation and the jump phenomenon were found, which were related to the particle size, immersed depth, and excitation amplitude (Liu et al., 2018). Similarly, Hong et al. (2018) performed an experiment on the particle-plate coupling system and discovered that the plate is more prone to jump under the action of small particle size. Further studies showed that the jump phenomenon was caused by the mutation in equivalent mass, stiffness, and damping, and the negative mass effect was first discovered in the particle-plate coupling system and explained by a simplified spring-mass model (Hong et al., 2021; Pan and Li, 2023). A mapping method was presented to study the dynamic behaviors of the vibrating plate-particle coupling system by partitioning the continuum into a mesh of mass points and simulated by the DEM, where the mass, stiffness, and damping were mapped from the global mass, stiffness, and Rayleigh damping matrix in the FEM (Pan et al., 2021). Liu et al. (2021) proposed a dilated-polyhedron-based DEM for the analysis of ice resistance on ship hulls during escort operations in level ice. This method can effectively simulate the irregular geometric shape, breaking process, and contact behavior of sea ice with the ship hulls. The nonlinear restoring force (NRF) is mainly generated by collision and friction between particles and vibrating structures. NRF affects the stability and diversity of GSCVS and is one of the main factors leading to complex dynamic behaviors of the system, such as the arrangement, fluidization, and accumulation of particle materials, as well as the bifurcation, period doubling, and chaos of structural responses (Dong et al., 2022; Lin et al., 2022). Moreover, NRF can serve as an indicator of system damage and failure. By monitoring and identifying the changes in NRF, the abnormal state of the system can be detected and prevented in time (He et al., 2012). Therefore, investigating the NRF between granular media and structures can provide a new perspective for the nonlinear dynamic behavior analysis and optimal control of GSCVS. However, no general method has been reported to identify the NRF in such systems. A major challenge is that the high discreteness and multiphase nature of granular materials make the NRF difficult to parameterize.

The analysis of traditional nonlinear dynamical system problems relies on the dynamical model describing the physical mechanism and its solution. However, as the research on complex dynamical systems deepens, the model dimension and computational cost increase significantly, and it becomes difficult to obtain the critical characteristic model of complex dynamical systems. In this context, the nonparametric model method based on data science is expected to be an effective way to simplify the complex modeling process and reduce the computational cost. One of the first developed nonparametric identification methods is the restoring force surface method, proposed more than 40 years ago, to identify and express the nonlinear system characteristics in terms of orthogonal functions, which was first developed for a broad class of single-degree-of-freedom dynamic systems, and then was extended for multi-degree-of-freedom (MDOF) systems (Masri and Caughey, 1979; Masri et al., 1982). In the following years, Yang and Ibrahim (1985) developed a relatively simple nonparametric technique for identifying nonlinearities in various discrete nonlinear vibrating systems, using power series expansions. Masri et al. (2004, 2006) extended their previous work to propose a general data-based approach for developing reduced-order, nonparametric models of nonlinear MDOF systems, using power series fitting techniques. Using Chebyshev polynomials, the nonparametric restoring force mapping method was applied to identify the nonlinear behavior of a passive elasto-magnetic suspension system (Bonisoli and Vigliani, 2007). Yun et al. (2009) applied model-free identification techniques based on nonparametric system identification methods to detect and interpret the changes of the nonlinear system, and to quantify the uncertainty associated with them. Recently, various nonparametric models have been used to identify the NRF of different nonlinear components, such as the power series and Legendre polynomials for magnetorheological (MR) dampers (Xu et al., 2012; Zhao et al., 2022), and the Chebyshev polynomials for shape memory alloys (Xu et al., 2015). Lei Y et al. (2021) proposed a two-stage nonparametric technique for identifying the NRF of MR dampers under unknown excitations. The NRF are proposed to be treated as “unknown fictitious forces” to the corresponding linear bare structure. Further, by considering the wavelet energy distribution in the nonlinear element become different from those of other linear elements, they employed the relative wavelet entropy to localize structural nonlinear elements. And the generalized Kalman filter with unknown input was used to reconstruct the characteristics of the NRF without its models (Lei et al., 2023). The actual hysteresis behavior of structures under granular dynamic load is hard to determine and model with parametric methods in engineering. Thus, inspired by the previous research, a general data-based and nonparametric approach can be devised to identify the NRF in the granules-structures coupled system. Moreover, the granular medium will generate various additional effects and impact with the structures, and its NRF differs from the nonlinear components previously studied, such as MR dampers and shape memory alloys. A new form of nonparametric model is required to estimate the NRF in the GSCVS.

The classical Kalman filter (KF) is a well-known effective algorithm for system identification, which can recursively estimate the full states from sparse measurements in real time (Kalman, 1960). However, it is restricted to linear systems with known structure and input, limiting its applicability. To overcome this, some researchers developed extended Kalman filter (EKF) methods that treat the structure parameters as part of the augmented state (Lei et al., 2012; Qi et al., 2024; Yao et al., 2015). In the following years, many alternative techniques have been proposed to address various engineering challenges that EKF is unable to cope with. Julier and Uhlmann (1997) proposed the unscented Kalman filter (UKF) to overcome the issues of linearization in the EKF, which employs a sampling technique to estimate the state distribution of the system. Mariani and Ghisi (2007) demonstrated that the UKF performs better than the EKF in nonlinear parameter identification problems, albeit with a higher computational burden. To cope with the high computational costs, Azam et al. (2012) suggested a parallel implementation of the UKF. Mu and Yuen (2015) proposed a novel robust outlier-resistant EKF for online parameter identification of linear time-varying structural systems. This solves the issue of parameter instability and enables online and continuous monitoring. In recent years, some KF-based system identification methods under unknown inputs have been developed. Lourens et al. (2012) proposed an augmented Kalman filter (AKF) to identify forces, where the unknown forces were included in the augmented state vector, and a technique similar to the EKF was used to estimate the state and the unknown forces. Wei et al. (2022) also presented an AKF based on sparse constraint theory, which assumed a random walk model for the unknown inputs. In order to reduce the negative effects of linearization in traditional Kalman-based methods, Al-Hussein and Haldar (2015) combined UKF with iterative least squares technology and proposed an unscented Kalman filter with unknown input (UKF-UI). Lei et al. (2019) developed a novel UKF for recursively identifying the state-input-system of nonlinear systems. The method is real-time and requires only partial observations of the responses. In subsequent studies, they proposed a method called smoothing EKF with unknown input without direct feedthrough, which can recursively identify the structural states, parameters, and unknown seismic excitations without direct feedthrough simultaneously (Lei et al., 2023). Huang et al. (2024) proposed an adaptive generalized EKF with unknown input based on the adaptive discrete state equation and minimum variance unbiased estimation principle. The algorithm eliminates the limitation that the real-time performance is restricted by whether the measurement equation has a direct feedthrough of the input, and realizes the optimization of the state and input estimates in the sense of minimum variance. Currently, the system identification methods based on KF and EKF, including the studies mentioned above, mainly start from the measured translational acceleration or displacement responses to estimate the system states, identify the system parameters and unknown inputs. However, some engineering structures may also generate rotational responses under external loads, which are often excluded from the observation vector due to measurement difficulties in practical engineering. Lei et al. (2012) also highlighted this problem. Some prior studies have demonstrated that the absence of response components in the observation vector of EKF may induce drift in the identification results (Liu et al., 2016; Xu et al., 2022). Therefore, if the observation vector containing the complete response components can be obtained by the existing measurement methods, it will help to improve the accuracy of the identification results, especially for civil structures with high noise.

The NRF in the GSCVS is difficult to parameterize due to the discreteness and multiphase nature of granular materials, which impedes the advancement of GSCVS research. This paper proposes a data-based nonparametric method to identify the NRF in the GSCVS, using a nonparametric model that accounts for the additional effects of the granules on both sides of the structures. This method reconstructs the observation vector of the conventional EKF by using the sparse measurement of the strain response. The reconstructed observation vector contains three response components: translational displacement, translational acceleration, and rotational acceleration. A numerical study was conducted on a simply supported beam under forced vibration. A pulse generator (PG) and a nonlinear damper (ND) were attached to the beam to represent the impact and additional effects of discrete media, respectively. The results of the data-based nonparametric method are compared with those of existing methods. An experimental study was performed on a simply supported beam immersed in granular media. The state estimation and NRF identification of the system were carried out by applying the proposed data-based nonparametric method. The dynamic characteristics of NRF and the effects of excitation amplitude, particle size, and immersed depth on NRF were analyzed.

2. Identification algorithm

2.1. Nonparametric modeling for NRF in the GSCVS

In the study of granular-structures coupled systems, researchers often simplify the complex modeling process by treating the particulate medium in an equivalent manner. For instance, Michon et al. (2013) considered the soft particles inside a honeycomb cantilever beam as oscillators with viscous damping to study the damping efficiency of the particles. Li et al. (2023) utilized several inertia-viscous dampers to equivalently represent the additional inertia and viscous effects of granules on the beam. Although these equivalent models cannot fully reveal the dynamic characteristics of the particulate medium, they provide an effective approach to study the interaction between particles and structures, which helps to reveal the nonlinear dynamic behaviors of the structures.

Inspired by the aforementioned research, we consider the effects of granular materials act on a vibrating structure as a series of concentrated NRF, and study the dynamic characteristics of the structure by identifying these NRFs. Figure 1(a) shows the coupled vibrating system composed of the granular medium and the embedded simply supported beam. To equivalently represent the granules-beam coupled system, multiple NRFs are incorporated into the finite element model of the simply supported beam, as depicted in Figure 1(b).

Figure 1.

(a) Particles-beam coupled vibrating system. (b) The equivalent finite element model of the particles-beam coupled vibrating system.

Then the dynamic equation of the system can be expressed as equation (1)

$\mathit{M}\ddot{\mathit{x}}+\mathit{C}\dot{\mathit{x}}+\mathit{K}\mathit{x}+\mathit{\psi }\mathit{R}\left[\dot{\mathit{x}},\mathit{x},\mathit{\eta }\right]=\mathit{\varphi }\mathit{f}$

(1)

In which $\ddot{\mathit{x}}$ , $\dot{\mathit{x}}$ , and $\mathit{x}$ are vectors of acceleration, velocity, and displacement of the beam, respectively, $\mathit{M}$ is the mass matrix, $\mathit{C}$ is the damping matrix, and $\mathit{K}$ is the stiffness matrix, respectively, $\mathit{R}\left[\dot{\mathit{x}},\mathit{x},\mathit{\eta }\right]$ represents the NRF vector, which reflects the effect of particulate matter on the vibrating beam, $\mathit{\eta }$ is the coefficients vector to be identified, $\mathit{f}$ is the excitation vector, $\mathit{\psi }$ and $\mathit{\varphi }$ are the influence matrix associated with $\mathit{R}\left[\dot{\mathit{x}},\mathit{x},\mathit{\eta }\right]$ and $\mathit{f}$ respectively.

The effects of the upper and lower granular medium on the immersed vibrating structures are generally different. The lower particles usually produce additional stiffness and damping effects, while the upper particulate matter produce additional mass effects and will collide with the structures in a wide range of excitation frequency. Considering the aforementioned characteristics in the GSCVS, the NRF can be represented in two parts, one of which reflects the additional stiffness, damping, mass and nonlinear effects, and the other represents the impact of the upper particles. The additional effects in GSCVS are similar to the nonlinear characteristics of an improved Dahl model (Zhou and Qu, 2002). Xu et al. (2012) have estimated the NRF generated by this model in the form of a power series polynomial, and numerical simulations have verified the accuracy of this method. On the other hand, in combination with our previous research, the impact effect can be fitted with a double exponential function (Dong et al., 2022; Wu, 2014). Taking into account the above two aspects, we propose a new nonparametric model for estimating the NRF in GSCVS, as shown in equation (2)

${\mathit{R}}_{\mathbf{i}}\left[\dot{\mathit{x}},\mathit{x},\mathit{\eta }\right]=\sum _{m=0}^{\text{U}}\sum _{n=0}^{\text{U}}{\eta }_{i,m,n}{\dot{\mathit{x}}}_i^m{\mathit{x}}_i^n{+\eta }_i^a{\text{e}}^{{(\eta }_i^b{\dot{\mathit{x}}}_i)}{+\eta }_i^c{\text{e}}^{{(\eta }_i^d{\mathit{x}}_i)}$

(2)

where ${\mathit{R}}_{\mathbf{i}}\left[\dot{\mathit{x}},\mathit{x},\mathit{\eta }\right]$ is the NRF acting on the ith node of the beam, ${\dot{\mathit{x}}}_i$ and ${\mathit{x}}_i$ are velocity and displacement vectors of the ith node, respectively. In this paper, ${\mathit{x}}_i^n$ denotes taking each element in ${\mathit{x}}_i$ with the power of n. ${\eta }_{i,m,n}$ , ${\eta }_i^a$ , ${\eta }_i^b$ , ${\eta }_i^c$ , and ${\eta }_i^d$ are the coefficients of the terms that constitute the NRF, the integers m and n satisfy the condition that $m+n\ne 0$ , the integer U depends on the degree of nonlinearity of the system. In the numerical validations and experimental investigations of this paper, the excitation applied to the beam induces only micro-vibrations. Therefore, the effects of higher-order nonlinear terms are not considered, and we set U = 2. Equation (2) can be specifically expressed as

${\mathit{R}}_{\mathbf{i}}\left[\dot{\mathit{x}},\mathit{x},\mathit{\eta }\right]={\eta }_{i,0,0}\mathit{j}+{\eta }_{i,0,1}{\mathit{x}}_i+{\eta }_{i,0,2}{\mathit{x}}_i^2+{\eta }_{i,1,0}{\dot{\mathit{x}}}_i{+\eta }_{i,1,1}{\dot{\mathit{x}}}_i{\mathit{x}}_i{+\eta }_{i,1,2}{\dot{\mathit{x}}}_i{\mathit{x}}_i^2{+\eta }_{i,2,0}{\dot{\mathit{x}}}_i^2{+\eta }_{i,2,1}{\dot{\mathit{x}}}_i^2{\mathit{x}}_i+{\eta }_{i,2,2}{\dot{\mathit{x}}}_i^2{\mathit{x}}_i^2{+\eta }_i^a{\text{e}}^{{(\eta }_i^b{\dot{\mathit{x}}}_i)}{+\eta }_i^c{\text{e}}^{{(\eta }_i^d{\mathit{x}}_i)}$

(3)

where j is a vector with all elements being 1.

2.2. EKF-OVR based on measured strain response

EKF is an effective method for dealing with inverse problems in engineering structures. It identifies system parameters or unknown inputs through observation vectors containing measured responses. Related studies have shown that for classical shear structures and truss structures (the response vector $x$ only contains translational responses), the identification results of unknown forces and displacements may show significant drifts when the observation vector only has translational acceleration responses. Adding translational displacement responses to the observation vector or using dummy-measurements can eliminate the drifts (Naets et al., 2015; Xu et al., 2022). This means that translational displacement responses or dummy-measurements can provide more information for inverse problems than translational acceleration responses alone. The response components in the observation vector significantly influence the identification results. However, in the existing studies that use EKF or its variants to deal with the inverse problem of continuum structures, the observation vector generally contains only the translational response measured by accelerometers but not the rotational response (Lei et al., 2012; Lei et al., 2012; Mu and Yuen, 2015; Xu et al., 2022), which is mainly because it is still a problem to have appropriate instruments for measuring dynamic responses of rotational angles. The inverse problem research of systems with strong nonlinearity, such as immersed continuous, is easily disturbed by noise. If the observation vector lacks rotational response components, it is more likely to make the identification results unstable.

This paper proposes an observation vector reconstruction method for EKF, which can obtain an observation vector that contains both translational displacement response, translational acceleration response, and rotational acceleration response when only strain response is measured. One takes a slight deflection simply supported beam system as an example and derives the observation vector reconstruction process for EKF. By introducing an extended state vector $\mathit{Z}=\{{\mathit{x}}^{\text{T}},{{\dot{\mathit{x}}}^{\text{T}},{\mathit{\eta }}^{\text{T}}\}}^{\text{T}}$ , and considering modeling error, equation (1) can be re-written in the state space form as equation (4)

$\dot{\mathit{Z}}=\left\{\begin{array}{l} \dot{\mathit{x}}\\ {\mathit{M}}^{-1}\left(\mathit{\varphi }\mathit{f}-\mathit{C}\dot{\mathit{x}}-\mathit{K}\mathit{x}-\mathit{\psi }\mathit{R}\mathbf{[}\dot{\mathit{x}}\mathbf{,}\mathit{x}\mathbf{,}\mathit{\eta }\mathbf{]}\right)\\ \mathbf{ }\mathbf{ }\mathbf{ }\mathbf{ }\mathbf{ }\mathbf{ }\mathbf{ }\mathbf{ }\mathbf{ }\mathbf{ }\mathbf{ }\mathbf{ }\mathbf{ }\mathbf{ }\mathbf{ }\mathbf{ }\mathbf{ }\mathbf{ }\mathbf{ }\mathbf{ }\mathbf{ }\mathbf{ }\mathbf{ }\mathbf{ }\mathbf{ }\mathbf{0}\end{array}\right\}=\mathit{g}\left(\mathit{Z},\mathit{f}\right)+\mathit{w}$

(4)

where $\mathit{g}\left(·\right)$ is a nonlinear function and $\mathit{w}$ is the model noise with zero mean and a covariance matrix $\mathit{Q}$ . In the traditional application of EKF, the nonlinear discrete observation equation can generally be expressed by equation (5)

${\mathit{y}}_{k+1}={\ddot{\mathit{u}}}_{k+1}^m=\mathit{\Pi }\mathit{h}\left({\mathit{Z}}_{k+1},{\mathit{f}}_{k+1}\right)+{\mathit{v}}_{k+1}$

(5)

where ${\mathit{y}}_{k+1}$ is an observation vector at time $t=\left(k+1\right)\mathrm{\Delta }t$ , $\mathrm{\Delta }t$ is the sampling time step, ${\ddot{\mathit{u}}}_{k+1}^m$ is the translational acceleration response in the deflection direction of a few measurement points, $\mathit{\Pi }$ is the position influence matrix, $\mathit{h}\left(·\right)$ is a nonlinear function, ${\mathit{v}}_{k+1}$ is a measurement noise vector assumed as a Gaussian noise with a covariance matrix ${\mathit{R}}_{k+1}$ . However, as just mentioned, when the observation vector contains only translational response, the error of identification results will be caused by the absence of response components. Therefore, we reconstruct the observation vector to achieve that more response components can be included in the observation vector when only a partial strain response is measured.

A simply supported transverse forced vibration beam with a slight deflection at time $t=\left(k+1\right)\mathrm{\Delta }t$ is shown in Figure 2(a), N+2 strain measurement points are established on different beam positions. Taking O as the origin, the s-axis is established along the axial of the undeformed beam, and the u-axis is found along the deflection direction of the beam.

Figure 2.

(a) A forced-vibration simply supported beam with a slight deflection. (b) Reconstruction of the response of the simply supported beam.

From the mechanics of materials, a slightly deformed beam satisfies the relationship that $\mu =\left(2/\delta \right){\epsilon }_s$ , where ε _s is the surface strain of the beam, μ is the corresponding curvature of the neutral layer, and $\delta$ represents the thickness of the beam. The Region I in Figure 2(a) is a micro segment of the simply supported beam, θ is the angle between the tangent of the micro segment and the s-axis, l represents the arc length parameter of the beam. The micro segment of a small deflection beam exhibit the following geometric relationships: ${\left(\text{d}l\right)}^2={\left(\text{d}u\right)}^2+{\left(\text{d}s\right)}^2$ , $\text{d}l=\sqrt{{\left(\text{d}u/\text{d}s\right)}^2+1}\text{d}s$ , and $\theta =\text{d}u/\mathrm{d}s$ . Considering that the θ of a small deflection beam is extremely small, ${\left(\text{d}u/\text{d}s\right)}^2$ is negligible compared to 1, and therefore $\text{d}l=\text{d}s$ . Moreover, the curvature μ and the bending moment M satisfy $\mu =M/EI$ , E and I are elasticity modulus and inertia moment of the beam. When E and I remain unchanged along the length of the beam, the change of curvature along the arc length can be expressed as

$\frac{\text{d}\mu }{\text{d}l}=\frac{\text{d}\mu }{\text{d}s}=\frac{1}{EI}\frac{\text{d}M}{\text{d}s}$

(6)

Since the bending moment between adjacent concentrated loads varies linearly along the length of the beam, it can be known from equation (6) that when all the concentrated loads have strain measurement points set up nearby, the curvature between adjacent strain measurement points can be approximated as linearly varying along the arc length. Therefore, the curvature function of the deflection curve between strain measurement points i and i+1 can be expressed as $\mu \left(l\right)={\alpha }_{i,i+1}l+{\beta }_{i,i+1}$ , ${\alpha }_{i,i+1}$ and ${\beta }_{i,i+1}$ represent the linear interpolation coefficients. Then the curvature at time $t=\left(k+1\right)\mathrm{\Delta }t$ can be discrete expressed as

${\mathit{\mu }}_{k+1}=\text{diag}\left({\mathit{\alpha }}_{k+1}\right)\hspace{0.17em}\mathit{l}+{\mathit{\beta }}_{k+1}$

(7)

where ${\mathit{\mu }}_{k+1}={\left({\mu }_{k+1}\left(0\right),{\mu }_{k+1}\left(1\right),{\mu }_{k+1}\left(2\right),{\mu }_{k+1}\left(3\right),···,{\mu }_{k+1}\left(i\right),···{\mu }_{k+1}\left(N+1\right)\right)}^{\text{T}}$ , ${\mu }_{k+1}\left(i\right)$ is the neutral layer curvature at the ith measurement point, the numbers 0 and N+1 represent the left and right endpoints of the beam, respectively, $\mathit{l}={\left(l\left(0\right),l\left(1\right),l\left(2\right),l\left(3\right),···,l\left(i\right),···,l\left(N+1\right)\right)}^{\text{T}}$ , $l\left(i\right)$ is the arc length between the left endpoint and the ith measurement point, ${\mathit{\alpha }}_{k+1}={\left({\alpha }_{k+1}^{0,1},{\alpha }_{k+1}^{0,1},{\alpha }_{k+1}^{1,2},{\alpha }_{k+1}^{2,3},{\alpha }_{k+1}^{3,4},···,{\alpha }_{k+1}^{i,i+1},···,{\alpha }_{k+1}^{N,N+1}\right)}^{\text{T}}$ , ${\mathit{\beta }}_{k+1}={\left({\beta }_{k+1}^{0,1},{\beta }_{k+1}^{0,1},{\beta }_{k+1}^{1,2},{\beta }_{k+1}^{2,3},{\beta }_{k+1}^{3,4},···,{\beta }_{k+1}^{i,i+1},···,{\beta }_{k+1}^{N,N+1}\right)}^{\text{T}}$ , ${\alpha }_{k+1}^{i,i+1}$ and ${\beta }_{k+1}^{i,i+1}$ are the linear interpolation coefficients of the arc segment between the measurement points i and i+1, which can be calculated by the measured strain and expressed as follows

${\alpha }_{k+1}^{i,i+1}=\frac{2\left({\epsilon }_{k+1}^s\left(i+1\right)-{\epsilon }_{k+1}^s\left(i\right)\right)}{h\left(l\left(i+1\right)-l\left(i\right)\right)}$

(8a)

${\beta }_{k+1}^{i,i+1}=\frac{2\left({\epsilon }_{k+1}^s\left(i\right)l\left(i+1\right)-{\epsilon }_{k+1}^s\left(i+1\right)l\left(i\right)\right)}{h\left(l\left(i+1\right)-l\left(i\right)\right)}$

(8b)

where ${\epsilon }_{k+1}^s\left(i\right)$ is the surface strain in the s direction at the position of the ith measurement point. The L in Figure 2(b) is a deflection curve of the beam at time $t=\left(k+1\right)\mathrm{\Delta }t$ , the angle between the tangent at the left endpoint of L and s-axis is ${\theta }_{k+1}^0$ . Because ${\theta }_{k+1}^0$ is unknown, the deflection curve equation of L cannot be reconstructed directly. Therefore, the deflection curve L is rotated counterclockwise around point o until it becomes tangent to the s-axis at point o, and the $L^{\prime }$ is obtained, the angle of rotation is ${\theta }_{k+1}^0$ . $L^{\prime }$ is reconstructed at first, and then the deflection curve equation of L can be obtained by applying the boundary conditions at the right side of the simply supported beam.

Since the definition of curvature is $\mu =\mathrm{d}\theta /\text{d}l$ , θ can be calculated by integrating the curvature along the arc: $\theta \left(l\right)=\int \mu \left(l\right)\mathrm{d}l=\left(1/2\right){\alpha }_{i,i+1}{l}^2+{\beta }_{i,i+1}l+{\gamma }_{i,i+1},l\in \left(l\left(i\right),l\left(i+1\right)\right)$ , ${\gamma }_{i,i+1}$ is an integral constant. Then the angle of the tangent line of $L^{\prime }$ at different measurement points concerning the s-axis can be expressed as equation (9)

${{\mathit{\theta }}_{k+1}^{\epsilon }}^{\prime }=\left(1/2\right)\text{diag}\left({\mathit{\alpha }}_{k+1}\right)\text{diag}\left(\mathit{l}\right)\mathit{l}+\text{diag}\left({\mathit{\beta }}_{k+1}\right)\mathit{l}+{\mathit{\gamma }}_{k+1}$

(9)

where ${{\mathit{\theta }}_{k+1}^{\epsilon }}^{\prime }=({\theta }_{k+1}^{\epsilon }{\left(0\right)}^{\prime },{\theta }_{k+1}^{\epsilon }{\left(1\right)}^{\prime },{\theta }_{k+1}^{\epsilon }{\left(2\right)}^{\prime },{\theta }_{k+1}^{\epsilon }{\left(3\right)}^{\prime },···,{\theta }_{k+1}^{\epsilon }{\left(i\right)}^{\prime },···{\theta }_{k+1}^{\epsilon }{\left(N+1\right)}^{\prime }{)}^{\text{T}}$ , ${\theta }_{k+1}^{\epsilon }{\left(i\right)}^{\prime }$ is the angle of the tangent line of l₂ at the ith measurement point concerning the s-axis, ${\mathit{\gamma }}_{k+1}={\left({\gamma }_{k+1}^{0,1},{\gamma }_{k+1}^{0,1},{\gamma }_{k+1}^{1,2},{\gamma }_{k+1}^{2,3},{\gamma }_{k+1}^{3,4},···,{\gamma }_{k+1}^{i,i+1},···,{\gamma }_{k+1}^{N,N+1}\right)}^{\text{T}}$ , ${\gamma }_{k+1}^{i,i+1}$ is an integral constant between measurement points i and i+1. Since the left end of $L^{\prime }$ is tangent to the s-axis, therefore ${\theta }_{k+1}^{\epsilon }{\left(0\right)}^{\prime }$ = 0, and then ${\gamma }_{k+1}^{0,1}$ = 0. Based on the continuity condition of the rotation angle at the measurement point, the calculated rotation angles by the rotation angle equations of the arc segments on either side of the measurement point should be equal. Consequently, the ${\gamma }_{k+1}^{i,i+1}$ can be sequentially determined using systems of simultaneous equations. For example, by solving equation (10)

$\{\begin{array}{l}{\theta }_{k+1}^{\epsilon }{\left(1\right)}^{\prime }=\left(1/2\right){\mathrm{\alpha }}_{k+1}^{0,1}{l}^2\left(1\right)+{\beta }_{k+1}^{0,1}l\left(1\right)+{\gamma }_{k+1}^{0,1}\\ {\theta }_{k+1}^{\epsilon }{\left(1\right)}^{\prime }=\left(1/2\right){\mathrm{\alpha }}_{k+1}^{1,2}{l}^2\left(1\right)+{\beta }_{k+1}^{1,2}l\left(1\right)+{\gamma }_{k+1}^{1,2}\end{array}$

(10)

It can be obtained that ${\gamma }_{k+1}^{1,2}=\left(1/2\right)\left({\mathrm{\alpha }}_{k+1}^{0,1}-{\alpha }_{k+1}^{1,2}\right)l^2\left(1\right)+\left({\beta }_{k+1}^{0,1}-{\beta }_{k+1}^{1,2}\right)l\left(1\right)+{\gamma }_{k+1}^{0,1}$ , and similarly, ${\gamma }_{k+1}^{2,3}=\left(1/2\right)\left({\mathrm{\alpha }}_{k+1}^{1,2}-{\alpha }_{k+1}^{2,3}\right)l^2\left(2\right)+\left({\beta }_{k+1}^{1,2}-{\beta }_{k+1}^{2,3}\right)l\left(2\right)+{\gamma }_{k+1}^{1,2}$ . Therefore, the general form of ${\gamma }_{k+1}^{i,i+1}$ can be expressed as equation (11)

${\gamma }_{k+1}^{i,i+1}={\displaystyle \sum }_{j=0}^{i-1}\left(\left(1/2\right)\left({\alpha }_{k+1}^{j,j+1}-{\alpha }_{k+1}^{j+1,j+2}\right)l^2\left(j+1\right)+\left({\beta }_{k+1}^{j,j+1}-{\beta }_{k+1}^{j+1,j+2}\right)l\left(j+1\right)\right)$

(11)

It is known from the geometric relationship in the Region I of Figure 2(a) that du = sinθdl, ds = cosθdl. Therefore, $u_{k+1}^{\epsilon }{\left(1\right)}^{\prime }={\int }_{l\left(0\right)}^{l\left(1\right)}\hspace{0.17em}\sin \left(\left(1/2\right){\mathrm{\alpha }}_{k+1}^{0,1}{l}^2+{\beta }_{k+1}^{0,1}l+{\gamma }_{k+1}^{0,1}\right)\text{d}l$ , and $u_{k+1}^{\epsilon }{\left(2\right)}^{\prime }=u_{k+1}^{\epsilon }{\left(1\right)}^{\prime }+{\int }_{l\left(1\right)}^{l\left(2\right)}\hspace{0.17em}\sin \left(\left(1/2\right){\mathrm{\alpha }}_{k+1}^{1,2}{l}^2+{\beta }_{k+1}^{1,2}l+{\gamma }_{k+1}^{1,2}\right)\text{d}l$ . Then the u-axis coordinates of the ith measurement point on $L^{\prime }$ represented by $u_{k+1}^{\epsilon }{\left(i\right)}^{\prime }$ can be generally expressed as equation (12)

$u_{k+1}^{\epsilon }{\left(i\right)}^{\prime }={\displaystyle \sum }_{j=0}^{i-1}{\int }_{l\left(j\right)}^{l\left(j+1\right)}\hspace{0.17em}\sin \left(\left(1/2\right){\mathrm{\alpha }}_{k+1}^{j,j+1}{\lambda }^2+{\beta }_{k+1}^{j,j+1}\lambda +{\gamma }_{k+1}^{j,j+1}\right)\text{d}\lambda$

(12)

Similarly, the s-axis coordinates of the ith measurement point on $L^{\prime }$ represented by $s_{k+1}^{\epsilon }{\left(i\right)}^{\prime }$ can be expressed as equation (13)

$s_{k+1}^{\epsilon }{\left(i\right)}^{\prime }={\displaystyle \sum }_{j=0}^{i-1}{\int }_{l\left(j\right)}^{l\left(j+1\right)}\hspace{0.17em}\cos \left(\left(1/2\right){\mathrm{\alpha }}_{k+1}^{j,j+1}{\lambda }^2+{\beta }_{k+1}^{j,j+1}\lambda +{\gamma }_{k+1}^{j,j+1}\right)\text{d}\lambda$

(13)

According to the geometric relationship in Region II of Figure 2(b), $u_{k+1}^{\epsilon }\left(i\right)$ can be expressed by $u_{k+1}^{\epsilon }{\left(i\right)}^{\prime }$ , $s_{k+1}^{\epsilon }{\left(i\right)}^{\prime }$ , and ${\theta }_{k+1}^0$ , as shown in equation (14)

$u_{k+1}^{\epsilon }\left(i\right)=-s_{k+1}^{\epsilon }{\left(i\right)}^{\prime }\sin \left({\theta }_{k+1}^0\right)+u_{k+1}^{\epsilon }{\left(i\right)}^{\prime }\cos \left({\theta }_{k+1}^0\right)$

(14)

Since $L^{\prime }$ is obtained from L by rotating it through an angle of ${\theta }_{k+1}^0$ , the angle of the tangent line of L at the ith measurement point concerning the s-axis satisfy the following expressions

${\theta }_{k+1}^{\epsilon }\left(i\right)={\theta }_{k+1}^{\epsilon }{\left(i\right)}^{\prime }-{\theta }_{k+1}^0$

(15)

where ${\theta }_{k+1}^0=\arctan \left(u_{k+1}^{\epsilon }{\left(N+1\right)}^{\prime }/s_{k+1}^{\epsilon }{\left(N+1\right)}^{\prime }\right)$ . Then the rotational acceleration response ${\ddot{\mathit{\theta }}}_{k+1}^{\epsilon }$ and the translational acceleration response ${\ddot{\mathit{u}}}_{k+1}^{\epsilon }$ can be obtained by the central difference method. To include more response components in the observation vector, equation (5) can be rewritten as

${\mathit{y}}_{k+1}^{\epsilon }=\left\{\begin{array}{l}{\mathit{L}}^a{\ddot{\mathit{x}}}_{k+1}^{\epsilon }\\ {\mathit{L}}^d{\mathit{x}}_{k+1}^{\epsilon }\end{array}\right\}={\mathit{h}}^{\prime }\left({\mathit{Z}}_{k+1},{\mathit{f}}_{k+1}\right)+{\mathit{v}}_{k+1}$

(16)

where ${\mathit{h}}^{\prime }\left(·\right)$ is a nonlinear function, ${\mathit{y}}_{k+1}^{\epsilon }$ is the reconstructed observation vector containing both the translational displacement response, translational acceleration response, and the rotational acceleration response derived by the measured strain response ${\mathit{\epsilon }}_{k+1}^s$ through the above process. ${\ddot{\mathit{x}}}_{k+1}^{\epsilon }$ and ${\mathit{x}}_{k+1}^{\epsilon }$ are the reconstructed acceleration response and the displacement response of the beam. Assume that the finite element model of the simply supported beam we are studying contains M nodes. Apart from the nodes at both ends, which only consider rotational degrees of freedom, the other nodes consider both translational and rotational degrees of freedom. And then, ${\ddot{\mathit{x}}}_{k+1}^{\epsilon }=({\ddot{\theta }}_{k+1}^{\epsilon }\left(1\right),{\ddot{\mathrm{u}}}_{k+1}^{\epsilon }\left(2\right),{\ddot{\theta }}_{k+1}^{\epsilon }\left(2\right),\cdots ,{\ddot{u}}_{k+1}^{\epsilon }\left(M-1\right),{\ddot{\theta }}_{k+1}^{\epsilon }\left(M-1\right),{\ddot{\theta }}_{k+1}^{\epsilon }\left(M\right){)}^{\text{T}}$ , ${\mathit{x}}_{k+1}^{\epsilon }={\left({\theta }_{k+1}^{\epsilon }\left(1\right),u_{k+1}^{\epsilon }\left(2\right),{\theta }_{k+1}^{\epsilon }\left(2\right),\cdots ,u_{k+1}^{\epsilon }\left(M-1\right),{\theta }_{k+1}^{\epsilon }\left(M-1\right),{\theta }_{k+1}^{\epsilon }\left(M\right)\right)}^{\text{T}}$ . ${\mathit{L}}^a$ and ${\mathit{L}}^d$ are 2M-2 by 2M-2 positional influence matrices representing locations of the reconstructed acceleration and displacement, composed of 0 and 1. If ${\mathit{y}}_{k+1}^{\epsilon }$ is reconstructed through the observed strain responses at nodes i and j, then it follows that ${\mathit{L}}_{\hspace{0.17em}2i-2,2i-2}^{\hspace{0.17em}a}={\mathit{L}}_{\hspace{0.17em}2i-1,2i-1}^{\hspace{0.17em}a}={\mathit{L}}_{\hspace{0.17em}2j-2,2j-2}^{\hspace{0.17em}a}={\mathit{L}}_{\hspace{0.17em}2j-1,2j-1}^{\hspace{0.17em}a}=1$ , ${\mathit{L}}_{\hspace{0.17em}2i-2,2i-2}^{\hspace{0.17em}d}={\mathit{L}}_{\hspace{0.17em}2j-2,2j-2}^{\hspace{0.17em}d}=1$ , and all other elements in matrices ${\mathit{L}}^a$ and ${\mathit{L}}^d$ are zeros.

The above derivation process is conducted under simply supported boundary conditions. For other common boundary conditions of beams in engineering, such as cantilever beams, clamped-clamped beams, and beams with one end clamped and the other end hinged, this method is still applicable. This is because there is a rotational constraint at the clamped support boundary of the beam, that is, the rotational response at that end is zero. Therefore, there is no need to rotate the beam, and the above derivation process can be directly applied to reconstruct the response of the beam starting from the clamped support end. Since the above derivation process is based on the assumption of small deformations, and the one-dimensional deflection curve is reconstructed through the strain response in a single direction, the OVR method is currently applicable to one-dimensional continuous structures under small deformations. For the application of more complex structural forms, the method need to be improved.

Let ${\widehat{\mathit{Z}}}_{k|k}$ and ${\widehat{\mathit{Z}}}_{k+1|k}$ be the posterior estimation of ${\mathit{Z}}_k$ and the prior estimation of ${\mathit{Z}}_{k+1}$ , which represent the state vector at time k and k+1, respectively. Eq. (4) and Eq. (16) can be linearized at ${\widehat{\mathit{Z}}}_{k|k}$ and ${\widehat{\mathit{Z}}}_{k+1|k}$ by Taylor series expansion as follows

$\mathit{g}\left(\mathit{Z},\mathit{f}\right)\approx \mathit{g}({\widehat{\mathit{Z}}}_{k|k},\mathit{f})+{\mathit{G}}_{k|k}(\mathit{Z}-{\widehat{\mathit{Z}}}_{k|k});\hspace{0.17em}{\mathit{G}}_{k|k}={\frac{\partial \mathit{g}\left(\mathit{Z},\mathit{f}\right)}{\partial \mathit{Z}}|}_{\mathit{Z}={\widehat{\mathit{Z}}}_{k|k}}$

(17)

${\mathit{h}}^{\mathbf{\prime }}\left({\mathit{Z}}_{k+1},{\mathit{f}}_{k+1}\right)\approx {\mathit{h}}^{\mathbf{\prime }}({\widehat{\mathit{Z}}}_{k+1|k},{\mathit{f}}_{k+1})+{\mathit{H}}_{k+1|k}({\mathit{Z}}_{k+1}-{\widehat{\mathit{Z}}}_{k+1|k});\hspace{0.17em}{\mathit{H}}_{k+1|k}={\frac{\partial {\mathit{h}}^{\mathbf{\prime }}\left(\mathit{Z},\mathit{f}\right)}{\partial \mathit{Z}}|}_{\mathit{Z}={\widehat{\mathit{Z}}}_{k+1|k}}$

(18)

Similar to the conventional EKF in the first-time update (prediction) procedure

${\widehat{\mathit{Z}}}_{k+1|k}={\widehat{\mathit{Z}}}_{k|k}+{\int }_{k\mathrm{\Delta }t}^{\left(k+1\right)\mathrm{\Delta }t}\mathit{g}({\widehat{\mathit{Z}}}_{t|k},\mathit{f})\text{d}t$

(19)

The prediction error covariance matrix ${\widehat{\mathit{P}}}_{k+1|k}=\mathit{E}[({\mathit{Z}}_{k+1}-{\widehat{\mathit{Z}}}_{k+1|k})({\mathit{Z}}_{k+1}-{\widehat{\mathit{Z}}}_{k+1|k}{)}^{\text{T}}]$ , can be derived as

${\widehat{\mathit{P}}}_{k+1|k}={\mathit{\Phi }}_{k|k}{\widehat{\mathit{P}}}_{k|k}{\mathit{\Phi }}_{k|k}^{\text{T}}{+\mathit{Q}}_k$

(20)

where ${\mathit{\Phi }}_{k|k}\approx \mathit{I}+\mathrm{\Delta }t{\mathit{G}}_{k|k}$ , ${\widehat{\mathit{P}}}_{k|k}=\mathit{E}[({\mathit{Z}}_k-{\widehat{\mathit{Z}}}_{k|k})({\mathit{Z}}_k-{\widehat{\mathit{Z}}}_{k|k}{)}^{\text{T}}]$ , I is a unit matrix. The second process of the EKF is the measurement update (correction) procedure, in which

${\widehat{\mathit{Z}}}_{k+1|k+1}={\widehat{\mathit{Z}}}_{k+1|k}+{\mathit{K}}_{k+1}[{\mathit{y}}_{k+1}^{\epsilon }-{\mathit{h}}^{\mathbf{\prime }}({\widehat{\mathit{Z}}}_{k+1|k},{\mathit{f}}_{k+1})]$

(21)

${\mathit{K}}_{k+1}={\widehat{\mathit{P}}}_{k+1|k}{\mathit{H}}_{k+1|k}^{\text{T}}{({\mathit{H}}_{k+1|k}{\widehat{\mathit{P}}}_{k+1|k}{\mathit{H}}_{k+1|k}^{\text{T}}{+\mathit{R}}_{k+1})}^{-1}$

(22)

${\widehat{\mathit{P}}}_{k+1|k+1}=\left(\mathit{I}-{\mathit{K}}_{k+1}{H}_{k+1|k}\right){\widehat{\mathit{P}}}_{k+1|k}{\left(\mathit{I}-{\mathit{K}}_{k+1}{\mathit{H}}_{k+1|k}\right)}^{\text{T}}+{\mathit{K}}_{k+1}{\mathit{R}}_{k+1}{\mathit{K}}_{k+1}^{\text{T}}$

(23)

where ${\mathit{K}}_{k+1}$ is the Kalman gain. From the above estimation by EKF for the extended state vector Z , the unknown coefficients vector $\mathit{\eta }$ can be identified, then the NRF of the granular materials act on the vibrating beam can be calculated by equation (3).

3. Numerical validations

To demonstrate and validate the performances of the proposed nonparametric model and EKF-OVR, numerical studies on a forced-vibration simply supported beam with a pulse generator (PG) and nonlinear damper (ND) are performed. The attached PG and ND represent the impact and additional effects of discrete media, respectively. The finite element model of the beam is shown in Figure 3, in which the beam is divided into 10 elements with equal lengths. Each element has two nodes with two DOFs in the vertical and rotational directions, respectively, except nodes 1 and 11 have only rotational DOF due to the constraints at the two ends. In this example, the beam is assumed to be in bending deformation without axial and shear forces.

Figure 3.

A finite element model of a forced-vibration simply supported beam with a pulse generator (PG) and nonlinear damper (ND).

The structure parameters are selected: $l_{Num}=0.22m$ , $b_{Num}=0.028m$ , $h_{Num}=0.0033m$ , $\rho =7800kg/m^3$ , $E=2.1\times 10^{11}\mathrm{P}\mathrm{a}$ , and $v=0.3$ , where $l_{Num}$ , $b_{Num}$ , $h_{Num}$ , $\rho$ , E and v are the length, cross-section width, cross-section thickness, mass density, elasticity modulus, and Poisson’s ratio of the beam, respectively. The structural damping is assumed as Rayleigh damping $\mathit{C}={\alpha }_{Num}\mathit{M}+{\beta }_{Num}\mathit{K}$ with ${\alpha }_{Num}=50.4$ and ${\beta }_{Num}=4.9\times 10^{-5}$ . The beam is subjected to an external harmonic excitation $f\left(t\right)$ with an amplitude of 10N and frequency 120 Hz in the vertical direction at node 6.

A PG is installed at node 4 and applies pulse excitations periodically, which is used to simulate the impact of the discrete media. The impact force function of the discrete medium on a vibrating structure has been obtained in our earlier work (Dong et al., 2022)

${\mathit{F}}^{\text{PG}}=A\hspace{0.17em}\exp \left({-\left(\left(\mathit{t}-∆t/2\right)/b\right)}^2\right)$

(24)

where A is the amplitude of the impact force, $\mathrm{\Delta }t$ is the impact duration, b is a coefficient related to the granular material. In this study, the parameters take the following values: $A=1.5N$ , $\mathrm{\Delta }t=0.0008\mathrm{s}$ , and $b=3.5\times 10^{-4}$ .

In addition to the PG, an ND is installed at node 4. The mechanical behavior of the ND is described by an improved Dahl model (Xu et al., 2012), and the mathematical expression is as follows

${\mathit{F}}_i^{\text{DH}}=k_0^{\text{DH}}{\mathit{u}}_i+c_0^{\text{DH}}{\dot{\mathit{u}}}_i+i_0^{\text{DH}}{\ddot{\mathit{u}}}_i+f_d^{\text{DH}}{\mathit{z}}_{\mathbf{i}}+f_0^{\text{DH}}$

(25)

${\dot{\mathit{z}}}_i=\sigma {\mathit{u}}_i\left(1-{\mathit{z}}_i\text{sgn}\left({\dot{\mathit{u}}}_i\right)\right)$

(26)

where $k_0^{\text{DH}}$ , $c_0^{\text{DH}}$ , and $i_0^{\text{DH}}$ are the stiffness coefficient, damper coefficient, and inertia coefficient of the ND, $f_d^{\text{DH}}$ represents the adjustable coulomb frictional force, $f_0^{\text{DH}}$ is the initial force, ${\mathit{z}}_i$ is a hysteretic parameter which describes the coulomb friction, $\sigma$ is the coefficient used to control the shape of the hysteretic curve, ${\mathit{u}}_i$ , ${\dot{\mathit{u}}}_i$ , and ${\ddot{\mathit{u}}}_i$ represent the displacement, velocity, and acceleration in the vertical direction of the node i. In determining the parameter values for equations (25) and (26), we referred to the research of Xu et al. (2012). We made adjustments based on the parameter values they choose, so that the displacement response of the beam in our numerical model is close to experiments by orders of magnitude. The parameters take the following values: $k_0^{\text{DH}}=3\times 10^3\mathrm{N}/\mathrm{m}$ , $c_0^{\mathrm{DH}}=200\mathrm{Ns}/\mathrm{m}$ , $i_0^{\text{DH}}=0.15\mathrm{kg}$ , $f_d^{\text{DH}}=5\mathrm{N}$ , $f_0^{\mathrm{DH}}=0$ , and $\sigma =5\times 10^4\mathrm{s}/\mathrm{m}$ . Equation (25) contains both the stiffness, damping, inertia, and nonlinear term, which aptly simulates the additional stiffness, damping, mass, and nonlinear effects mainly generated by the discrete medium.

The Newmark method and the relation between curvature and translational displacement are used to calculate the strain response at different beam nodes. The rotational acceleration response and the translational acceleration response are calculated by the central difference method. The sampling interval h is 2 × 10⁻⁴s, and the orders of magnitude of the translational acceleration response and the rotational acceleration response are not less than 10⁻² and 10⁻⁴, respectively. Consequently, the truncation error O(h²) introduced by employing the central difference method is inconsequential. Similarly, the truncation error in the experimental study is also negligible. Combine the NRF and the response of the beam at the initial moment, the unknown coefficients in the nonparametric model can be fitted using the least squares method, and the fitted values are taken as the initial values in the identification process. To consider the influence of measurement noise on the identification results, white noises are added to the computed responses to simulate the measured responses, which could be expressed as follows (Pan et al., 2017; Zhang et al., 2022)

${\mathit{\chi }}_n=\mathit{\chi }+lev\left(1/n_{\chi }\right){\displaystyle \sum }_{i=1}^{n_{\chi }}\left|{\chi }_i\right|\mathit{r}\mathit{a}\mathit{n}\mathit{d}$

(27)

where ${\mathit{\chi }}_n$ is the noise polluted signal vector, $\mathit{\chi }$ is the initially calculated signal vector, $n_{\chi }$ represents the number of the elements in $\mathit{\chi }$ , lev is the strength level of the noise, and the rand is a standard normal distribution vector.

To assess the identification accuracy of the displacements, the velocities, and the NRF time history, a relative root-mean-square error (RRMSE) is used (Xu et al., 2022), which is defined by the following equation

$\text{RRMSE}=\sqrt{\frac{\frac{1}{N_s}{\displaystyle \sum }_{k=1}^{N_s}{\left({\mathit{\gamma }}_k^I-{\mathit{\gamma }}_k^T\right)}^2}{\frac{1}{N_s}{\displaystyle \sum }_{k=1}^{N_s}{\left({\mathit{\gamma }}_k^T\right)}^2}}\times 100\%$

(28)

where N _s is the number of sampling points, ${\mathit{\gamma }}_k^I$ and ${\mathit{\gamma }}_k^T$ are the identified results and their true value at the kth time step. The true values of displacement and velocity are obtained by solving the finite element model shown in Figure 3 numerically. Then the true value of NRF can be computed by equations (24) and (25) after calculating the displacement, velocity, and acceleration responses of the beam.

Three strain measurements at the 4th, 6th, and 8th nodes are used for the identification. The identification results of NRF of node 4 are compared when the observation vector contains different response components, that is, only the translational acceleration response (EKF), the translational displacement response, and the translational acceleration response (EKF-DF), and simultaneously included the translational displacement response, the translational acceleration response and the rotational acceleration response (EKF-OVR), the latter three responses can be calculated by the measured strain response with equations (6)–(15). In addition, EKF-DF uses power series polynomials to model the NRF, while EKF-OVR uses a new expression in the form of equation (3) proposed in this paper. The accuracy of the identification methods can be evaluated by comparing the results obtained by different methods with the true values. Substitute the parameters identified by EKF-OVR into equation (3), the following equation is obtained

${\mathit{R}}_4\left[\dot{\mathit{x}},\mathit{x},\mathit{\eta }\right]=1.09\times 10^4\times \mathit{j}+3.59\times 10^7\times {\mathit{x}}_4+5.98\times 10^{10}\times {\mathit{x}}_4^2+1.76\times 10^3\times {\dot{\mathit{x}}}_4+1.11\times 10^6\times {\dot{\mathit{x}}}_4{\mathit{x}}_4-1.42\times 10^{11}\times {\dot{\mathit{x}}}_4{\mathit{x}}_4^2+1.77\times 10^4\times {\dot{\mathit{x}}}_4^2-2.25\times 10^8\times {\dot{\mathit{x}}}_4^2{\mathit{x}}_4-54.92\times {\mathrm{e}}^{\left(23.47\times {\dot{\mathit{x}}}_4\right)}-1.09\times 10^4\times {\mathrm{e}}^{\left(3.26\times 10^3\times {\mathit{x}}_4\right)\hspace{0.17em}}$

(29)

The temporal evolution of the identified NRF can be determined by employing equation (29). Figure 4 presents a comparative analysis of the time histories for both the true value and identified NRF using the EKF, EKF-DF, and EKF-OVR. The identified NRF using the EKF shows significant drifts, while those of the EKF-DF and EKF-OVR do not. Similar results have been reported in previous studies (Liu et al., 2016; Xu et al., 2022). These drifts are caused by the low-frequency components stemming from double integration errors (Azam et al., 2015). The enlargement of the dotted frame with light blue background in Figure 4 shows that EKF-DF and EKF-OVR are approximately roughly the true value under 3% noise, the RRMSEs of EKF-DF and EKF-OVR are 4.1% and 2.7%, respectively. Still, EKF-OVR is more accurate in impact recognition, and the conventional EKF and EKF-DF can hardly detect a low-intensity impact. The phase diagram shows how the NRF varies with different beam displacements, and the impact is framed by solid black lines. Figure 5 compares identified time histories of displacement and velocity with their true values. Similar to Figure 4, the displacement identified by EKF showed deep drifts. However, no drift is observed in the velocity identification results, and the RRMSEs of EKF, EKF-DF and EKF-OVR are 2.4%, 1.4%, and 0.3%, respectively.

Figure 4.

The NRF estimated at node 4 by EKF, EKF-DF, and EKF-OVR with partial measurements. The measurements are polluted by 3% noise.

Figure 5.

The displacement and velocity estimated at node 4 by EKF, EKF-DF, and EKF-OVR with partial measurements. The measurements are polluted by 3% noise.

Due to the interference of discrete media, the measured response of GSCVS is usually accompanied by large noise. Therefore, the noise of the measured strain response is increased to 10%, and the above simulation processes are repeated. The identification result of the NRF is

${\mathit{R}}_4\left[\dot{\mathit{x}},\mathit{x},\mathit{\eta }\right]=9.19\times 10^3\times \mathit{j}+2.55\times 10^7\times {\mathit{x}}_4+3.55\times 10^{10}\times {\mathit{x}}_4^2+488.74\times {\dot{\mathit{x}}}_4+2.05\times 10^5\times {\dot{\mathit{x}}}_4{\mathit{x}}_4-7.79\times 10^{10}\times {\dot{\mathit{x}}}_4{\mathit{x}}_4^2+3.99\times 10^3\times {\dot{\mathit{x}}}_4^2-1.06\times 10^8\times {\dot{\mathit{x}}}_4^2{\mathit{x}}_4-1.93\times {\mathrm{e}}^{\left(51.53\times {\dot{\mathit{x}}}_4\right)}-9.19\times 10^3\times {\mathrm{e}}^{\left(2.76\times 10^3\times {\mathit{x}}_4\right)\hspace{0.17em}}$

(30)

Figure 6 is the time history identification results of the NRF, and two local enlargements at time 0.46–0.48s and 0.94–0.96s are drawn. Comparison between Figures 6 and 4 illustrate the different effects of increasing noise on the three methods. EKF exhibits more significant drift. EKF-DF suffers from reduced accuracy, the RRMSE increased from 4.1% to 17.2%. And EKF-OVR maintains high recognition performance, the RRMSEs of 3% and 10% noise are 2.7% and 2.9%, respectively. However, the convergence rate of EKF-OVR decreases when the noise level increases. Figure 7 compares identified time histories of displacement and velocity with their true values. The displacement identified by EKF also drifted significantly.

Figure 6.

The NRF was estimated at node 4 by EKF, EKF-DF, and EKF-OVR with partial strain measurements. The measurements are polluted by 10% noise.

Figure 7.

The displacement and velocity estimated at node 4 by EKF, EKF-DF, and EKF-OVR with partial strain measurements. The measurements are polluted by 10% noise.

The following conclusions can be drawn based on the simulation results in this section. The traditional EKF tends to cause drifts of NRF and displacement identification results, even under low noise levels. When identifying NRF with 3% noise, the RRMSE of EKF-DF is 4.1%, while that of EKF-OVR is slightly lower at 2.7%. Both EKF-DF and EKF-OVR produce identification results that are close to the true value at low noise level. However, EKF-OVR can detect the low-intensity impact force, leading to a smaller RRMSE than EKF-DF. When identifying NRF with 10% noise, the drift phenomenon of EKF becomes more evident, and the RRMSE of EKF-DF increases to 17.2%. In contrast, the RRMSE of EKF-OVR only rises slightly to 2.9%. Therefore, EKF-DF suffers from a decline in identification ability as the noise level increases, whereas EKF-OVR preserves its noise robustness.

4. Experimental study of NRF in GBCVS

To reveal the characteristics of NRF in granules-beam coupled vibrating system (GBCVS), this paper conducted an experimental study on a simply supported beam immersed in granular media and identified the NRF of granules acting on the beam using the proposed EKF-OVR. The experimental setup of the GBCVS is shown in Figure 8, and the equipment names and parameters corresponding to the numerical expressions in the figure are as follows:

1. Lenovo T440P laptop, equipped with data acquisition and analysis software.

2. Donghua Test DH-5956 dynamic signal test and analysis system, offers features such as synchronized inputs (up to 16 channels), waveform outputs, 32-bit DSP processing, 24-bit ADC/DAC resolution.

3. The Beijing Purion precision electric technology DG1000 series signal generator, features direct digital frequency synthesis (DDS) technology, dual-channel output, 5 standard and 48 built-in arbitrary waveforms, with a maximum output frequency of 25 MHz, a sampling rate of 100MSa/s, and a vertical resolution of 14 bits.

4. The Denmark B&K power amplifier Type 2718 boasts a 75 VA power output capacity at 3 Ω, features a continuously variable current limit ranging from 1 A to 5 A (RMS), and offers a 40 dB voltage gain.

5. The Denmark B&K Type 4809 vibration exciter, with capable of delivering up to 60 N of sine peak force with air cooling, operates within a frequency range of 10 Hz to 20 kHz, and allows for a continuous peak-to-peak displacement of 8 mm

6. The Donghua Test DH3810N-2 strain gauge adapter features bridge completion resistors, electrical calibration circuitry, supports both two and three wire connections, and is compatible with strain gauges having gauge factors between 1.5 and 2.5.

7. Spherical glass beads, primarily composed of SiO₂ and Na₂O, have a density of 2.8 g/cm³ and a static angle of repose of 5°. The experiment utilizes three different particle sizes: 1.0–1.5 mm, 2.0–3.0 mm, and 4.0–5.0 mm.

8. The Denmark B&K Type8230-002 force sensor, features a sensitivity of 2.2 mV/N, a full-scale range of ±2000 N, and an operating temperature range of −55 to 125°C.

9. The strain measurement point, which is essentially a half-bridge circuit made up of strain gauges. The BA350-3BB strain gauge from AVIC electronic measurement company, features a nominal resistance of 350 Ω, a gauge factor range of 1.86 to 2.20.

10. A steel beam of size $220\text{mm}\times 28\text{mm}\times 3.3\text{mm}$ . The density is 7800 kg/m³, elastic modulus is $2.1\times 10^{11}\text{Pa}$ , and Poisson’s ratio is 0.3.

Figure 8.

(a) Experimental setup of the GBCVS. A steel beam is immersed in the athermal glass beads. The vibrator obtains sinusoidal signals from the signal generator and excites the beam through the connecting rod. The signal collector collects the force signals and strain signals. (b) Schematic diagram of the experimental setup of the GBCVS.

As shown in Figure 8(a), the left end of the beam is machined with a groove, while the right end is not. The manner in which the beam and the fixture are matched is shown in Figure 8(b). The fixture at the left end matches the groove, restricting its axial translational degree of freedom, while retaining its rotational degree of freedom. The right end of the beam has both axial translational and rotational degrees of freedom. This configuration constructs a simply supported boundary. The beam is immersed in the glass beads, and the immerged depth is defined as the distance between the upper surface of the particle pile and the upper surface of the beam. The vertical force, as $F=F_0\hspace{0.17em}\sin \left(2\mathrm{\pi }ft\right)$ , where $F_0$ is the excitation force amplitude and $f$ is the excitation frequency, acts on the site, which is 10 mm away from the middle point on the center axis to avoid the second modal node, and the frequencies range from 60 to 320 Hz. The force sensor is mounted on the steel shaft connecting the exciter and the beam to provide feedback on the excitation force. The strain measurement points located about 60 mm, 120 mm, and 180 mm away from the left end, respectively.

This paper considers the effects of granular materials on the vibrating beam as a series of concentrated NRF (Li et al., 2023; Michon et al., 2013). The finite element model of the experimental GBCVS can be shown in Figure 9, the steel beam is divided into 11 equal-length elements with 12 nodes, and the action of the particles around each node on the beam is reduced to a concentrated NRF. The external excitation is applied at the 6th node, and the strain gauges are attached at the 4th, 7th, and 10th nodes, respectively. A rough estimate of the NRF acting on the beam element at the initial moment is made, and combined with the system response, the least squares method is used to fit the coefficients in the nonparametric model, taking the fitted value as the initial value of the unknown parameters.

Figure 9.

Finite element model of the experimental GBCVS.

Figure 10 shows the identification results of node 8, where the external harmonic excitation has an amplitude of 8 N and a frequency of 205 Hz. The diameter of the athermal glass beads is 1–1.5 mm and the immersed depth of the beam is 35 mm. The time history of NRF is calculated by the identified equation (31)

${\mathit{R}}_8\left[\dot{\mathit{x}},\mathit{x},\mathit{\eta }\right]=2.87\times 10^4\times \mathit{j}+1.88\times 10^7\times {\mathit{x}}_8+6.13\times 10^9\times {\mathit{x}}_8^2-112.46\times {\dot{\mathit{x}}}_8+7.22\times 10^3\times {\dot{\mathit{x}}}_8{\mathit{x}}_8+5.49\times 10^8\times {\dot{\mathit{x}}}_8{\mathit{x}}_8^2-211.31\times {\dot{\mathit{x}}}_8^2-8.49\times 10^5\times {\dot{\mathit{x}}}_8^2{\mathit{x}}_8+1.31\times 10^9\times {\dot{\mathit{x}}}_8^2{\mathit{x}}_8^2+23.22\times {\mathrm{e}}^{\left(4.22\times {\dot{\mathit{x}}}_8\right)}-2.87\times 10^4\times {\mathrm{e}}^{\left(6.53\times 10^2\times {\mathit{x}}_8\right)\hspace{0.17em}}$

(31)

Figure 10.

The displacement, velocity, and NRF are estimated at node 8 with EKF-OVR.

The reference values of the displacement and velocity are calculated by the measured strain response with equations (6)–(15). The identification results of displacement and velocity obtained by EKF-OVR have RRMSE values of 3.2% and 2.4%, respectively, compared to their reference values. Prior to 0.1 s, notable deviation are observed in both the displacement and NRF identification results. In the process of employing EKF-OVR for the estimation of system state, it is imperative to set initial values for the extended state vector Z . Due to the inability to ascertain the true initial values of displacement, velocity, and the NRF, an error arises between the presupposed and the true initial values of Z . This error is the source of the initial deviation observed in the identification results. The enlarged image of NRF identification results clearly shows that EKF-OVR recognized the impact between the particles and the beam. The identification results of NRF at different positions on the beam are shown in Figure 11. The NRF near the midpoint is the largest and decreases toward both sides. This is mainly because the vibration amplitude of the beam is larger near the midpoint, which drives the surrounding particles to have higher velocity, that is, higher energy, and thus a larger NRF acting on the beam (Bai et al., 2009). The NRF on both sides of the beam is not symmetrical about the excitation because the excitation is not applied at the midpoint of the beam.

Figure 11.

(a) Identification results of NRF at different positions on the beam. (b) The projection of (a) onto the y-z plane reflects the amplitude of the NRF at different positions.

Figure 12 shows the variations of the velocity amplitude and NRF amplitude of node 6, as well as their period-doubling modes concerning the excitation amplitude, particle size, and immersed depth. The relative changes of the amplitudes of velocity and NRF are more concerned, rather than their absolute values. Therefore, the amplitudes of velocity and NRF are reduced by their values under the external excitation of 4N, 1.5 mm particle diameter, and 35 mm immersed depth, respectively. The reduced values are defined as reduction coefficients, then the amplitudes of velocity and NRF are plotted in the same coordinate system. At the same time, the different period-doubling modes of the velocity and NRF are marked. Firstly, both the velocity and NRF raised when the excitation amplitude increased from 4N to 8N, and both exhibited a steep rise at 8N. The disruption of the uniform motion of the particle system by the vibrating structures is one of the critical factors that create the steep rise (Pan et al., 2021). The period-doubling modes of the velocity response vary as period-1, period-2, and period-3, while the NRF stayed at period-1. Similar results are found in increasing particle size but no abrupt increase in velocity or NRF. The period-4 of the velocity response is observed when the particle diameter is 3 mm, and the NRF remains free of the period-doubling phenomenon. As the immersed depth of the beam increased (by adding more particles), the beam response and NRF showed an opposite trend, that is, the beam velocity decreased, and the NRF increased. This is because increasing the immersed depth implies that the mass of the particles acting on the beam increases. Even though the reduction of the beam velocity leads to a lower velocity of the surrounding particles, the total kinetic energy of the particles still increases, resulting in a larger NRF. Based on the NRF at different positions on the beam, as shown in Figure 11, and the velocity response and the NRF under different excitation amplitudes and particle diameters, we conclude that the NRF is positively correlated with the beam velocity when the mass of the particle pile interacting with the beam remains constant.

Figure 12.

The amplitude of velocity and NRF of node 6 varied with excitation amplitude, particle diameter, and immersed depth. Meanwhile, the period-doubling modes of the velocity and the NRF are reflected.

The forces of particles’ impact on a vibrating continuum have been studied when the continuum drives the overall motion of the particle pile at the bottom. Researchers observed that the impact strength, controlled by the normalized vibration acceleration, undergoes a series of period-doubling bifurcations, typically in the sequence of period-2, period-4, chaos, period-3, period-6, chaos, period-4, period-8, and chaos (Jiang et al., 2006). However, there has yet to be a public report on the characteristics of the force exerted by the granules on the vibrating structures immersed in the granular pile. The submerged beam is excited by maintaining a constant excitation frequency and varying the excitation amplitudes. The time history, phase trajectory, Poincaré map, and frequency spectrum of the beam response data are plotted to identify the period-doubling mode, as shown in Figure 13. Each column consists of four subfigures represents one working condition. As the excitation amplitude increases successively under a specific excitation frequency, a period-doubling bifurcation process occurs: period-1, period-2, period-4, and chaos. However, the corresponding NRF showed only period-1 and period-2, as shown in Figure 14. Compared with the coupling system composed of particles and the immersed vibrating structures, the unilateral excited granular bed presents various global motion phases such as solid-like, bouncing bed, and collect-and-collide (Yin et al., 2017; Zhang et al., 2017, 2023). The completely inelastic collision between the granule bed moving as a whole and the vibrating continuum leads to period-doubling bifurcations in the motion of the particle pile (Gilet et al., 2009; Melo et al., 1995), which result in period-doubling bifurcations of the impact force. However, the immersed vibrating structures only drive the neighboring particles, and the particle pile generally exhibits only local fluidization, hardly separating from the structures to show a global movement. Therefore, NRF shows no rich period-doubling phenomena. As shown in Figure 14, the NRF under a sinusoidal external excitation can be decomposed into fundamental and higher harmonic components. The fundamental component has the same frequency as the external driving force, whereas the higher harmonic components include high-frequency components such as two and three times the external excitation frequency. Figures 13 and 14, compared to show that when the external excitation is weak, the NRF mainly consists of the fundamental component, resulting in period-1 or period-2 of the beam response. As the amplitude of external excitation increases, higher harmonic components in the NRF cause period-doubling bifurcation in the beam response. When the external excitation amplitude reaches 10N, higher harmonic components dominate the NRF, leading to a chaotic state of the beam response.

Figure 13.

A process of the period-doubling bifurcation of the beam response in 1–1.5 mm particle diameter and 45 mm immersed depth. The excitation frequency remains unchanged at 200 Hz. The excitation amplitude increases at intervals from 4N to 10N. Each column consists of four subfigures represents one working condition. The time history, phase trajectory, Poincare map, and frequency spectrum are adopted to identify the period-doubling modes.

Figure 14.

The time history and frequency spectrum of the NRF in 1–1.5 mm particle diameter and 45 mm immersed depth. The excitation amplitude increases at intervals from 4N to 10N. The excitation frequency remains unchanged at 200 Hz. Each column consists of two subfigures represents one working condition.

5. Conclusions

The study of the nonlinear restoring force (NRF) between particles and structures can enhance our understanding on the mechanical properties of granules-structures coupled vibrating system (GSCVS), and provide a scientific basis for their design, optimization and management. However, due to the complex dynamic interactions between particles and structures, it is difficult to parametrically model the NRF. This paper proposed a data-based nonparametric method to estimate the NRF in the GSCVS. The proposed data-based nonparametric method has two novel features. Firstly, it developed a nonparametric model to estimate the NRF, accounting for the additional effects of particles on both sides of the structures. Secondly, it reconstructed the observation vector in the conventional EKF using sparse measurements of the strain response. The reconstructed observation vector consists of three response components: translational acceleration, translational displacement, and rotational acceleration. The rotational acceleration response is challenging to measure in engineering applications. A numerical study was conducted on the lateral forced vibration of a simply supported beam system with a nonlinear damper and a pulse generator. The identification results obtained by the proposed extended Kalman filter based on observation vector reconstruction (EKF-OVR) were compared with those obtained by the conventional EKF and the data fusion based EKF (EKF-DF). When identifying NRF with 3% measurement noise, the relative root-mean-square error (RRMSE) of EKF-DF was 4.1%, while that of EKF-OVR was slightly lower at 2.7%. As the noise increased to 10%, the RRMSE of EKF-DF rised to 17.2%, whereas the RRMSE of EKF-OVR only increased marginally to 2.9%. An experimental study on a simply supported beam system submerged in granular media was conducted, and the NRF in the granules-beam coupled vibrating system (GBCVS) was identified by the proposed EKF-OVR. The following conclusions can be drawn from the numerical and experimental studies:

1. Numerical results demonstrated that the proposed EKF-OVR effectively addressed the issue of identification results drift in conventional EKF. Furthermore, at low noise levels, both EKF-DF and EKF-OVR yielded identification results that were close to the true values. However, EKF-OVR was able to detect the low-intensity impact force, resulting in a smaller RRMSE than EKF-DF. As the noise level increased, EKF-DF suffered from a decline in identification ability, whereas EKF-OVR maintained its noise robustness.

2. As the external excitation amplitude increased, the higher harmonic components of the NRF of the GBCVS also increased, resulting in period-doubling and chaos in the beam response. Unlike the rich periodicity of impact forces between a unilaterally excited particle bed and the vibrating continuum, the NRF in GBCVS only observed period-1 and period-2.

This paper provides a general methodology for identifying NRF of one-dimensional continuous structures under small deformations while immersed in the granular medium, which can be used to analyze the nonlinear dynamic behavior of submerged one-dimensional continuous structures under small deformations and optimize the design of structures. Future studies should include further investigation of the performance of the proposed approach on structures with more complex shapes and boundary conditions.