A General Shock Waveform and Characterisation Method

TL;DR

This paper introduces a comprehensive shock waveform model based on the mechanical origins of shock, along with a signal processing method that effectively represents pyroshock and ballistic shocks in both time and frequency domains.

Contribution

It presents a novel shock waveform with parameters that have clear mechanical meanings and a processing method that captures shock signals efficiently in multiple domains.

Findings

Shock complexity correlates with the number of waveform components ($ta_{90\%}$).

Shock dominance is quantified by a normalized parameter ($ppa$).

The method effectively analyzes different categories of pyroshock and mechanical shocks.

Abstract

A shock waveform is proposed based on the mechanical mechanism of shock generation in a structure. The parameters in the shock waveform have clear mechanical meanings about the generation and development of the shock. A shock signal processing method is proposed and applied to represent a pyroshock or ballistic shock signal in both temporal and frequency domain using finite terms of the shock waveform components. It is found that complexity and dominant shock distance of a shock can be described quantitatively by the number of waveform components ($\eta_{90\%}$) and a normalised parameter ($\kappa$), respectively. Pyroshock and mechanical shocks in different categories are analysed using the proposed shock waveform to demonstrate its importance and generality in shock representation.

Figures

Tables

Table 1. Table 1: Shock environment classification based on κ 𝜅 \kappa

Category	$\kappa$ Range
Near-field shock
Mid-field shock
Far-field shock	$\kappa ⩾10$

Table 2. Table 2: Special Case of Shock Wave

Waveforms	Case in Shock Waveform	Shape
Proposed Waveform	No restriction on parameters	Including any following shape
Harmonic Wave $e^{i\omega t}$	$\zeta =0$
Prony Mode $e^{i\omega t-\zeta \omega t}$	$\zeta \ne 0$, $\tau =0$
Kern and Hayes’ Function $t{e}^{i\omega t-\zeta \omega t}$	$n=\tau \omega \zeta =1$
Advanced Prony Mode	$\zeta \ne 0,\tau \ne 0$
Asymmetric Wavelet	$\zeta \ne 0,\kappa \gg 0$
Symmetric Wavelet	$\zeta \ne 0,\kappa \gg 0,\phi =-\left(\tau \omega -\pi \lfloor \frac{\tau \omega }{\pi }\rfloor \right)$

Table 3. Table 3: Shock waveform components of NFS

Component	A($m/s^2$)	$\frac{\omega }{2\pi }$(Hz)	$\mathring{t}\left(ms\right)$	$\tau \left(ms\right)$	$\zeta$	$\phi$	$\kappa$	$ϵ(\%)$
1	3939.41	2077	0.53	0.15	0.066	5.93	0.32	52.48
2	1375.47	2534	1.52	0.08	0.027	3.31	0.21	9.22
3	2703.28	1543	0.21	0.02	0.201	0.56	0.04	7.56
4	2435.72	5469	0.59	0.22	0.212	3.03	1.23	6.59
5	968.94	3987	-10.89	13.03	1.001	0.01	51.99	4.16
6	1367.09	3289	1.16	0.00	0.042	0.37	0.01	3.80
7	2737.31	6056	-3.15	3.66	18.987	3.18	22.19	3.28
8	1939.25	6512	1.62	0.17	1.364	3.96	1.14	1.29
9	766.91	5896	2.13	0.18	0.050	1.89	1.07	1.27
10	781.65	9444	0.66	0.12	0.033	5.41	1.13	1.08
11	128.05	1016	-10.89	12.09	1.845	3.30	12.29	0.09
12	62.27	0.08	0.49	0.04	3333	6.22	0.00	0.01

Table 4. Table 4: Shock waveform components of MFS

Component	A($m/s^2$)	$\frac{\omega }{2\pi }$(Hz)	$\mathring{t}\left(ms\right)$	$\tau \left(ms\right)$	$\zeta$	$\phi$	$\kappa$	$ϵ(\%)$
1	3231.82	4683	0.99	1.3	0.056	6.01	6.11	30.02
2	3489.69	4315	-21.08	21.99	5.272	0.11	94.89	15.20
3	1206.31	22068	-3.93	5.46	0.015	1.90	120.69	7.14
4	1222.74	4086	1.84	1.07	0.021	5.87	4.39	7.13
5	1786.27	16446	-3.88	6.66	0.385	1.59	109.63	4.15
6	615.19	4863	4.25	0.53	0.005	3.44	2.60	3.21
7	2138.04	16626	0.09	0.04	0.021	5.46	0.71	2.89
8	762.57	5554	-18.25	19.51	0.239	6.11	108.40	2.52
9	1059.31	20996	-3.64	4.37	0.099	0.04	91.91	2.01
10	-693.35	3249	0.65	0.18	0.014	5.75	0.59	1.99

Table 5. Table 5: Shock waveform components of FFS

Component	A($m/s^2$)	$\frac{\omega }{2\pi }$(Hz)	$\mathring{t}\left(ms\right)$	$\tau \left(ms\right)$	$\zeta$	$\phi$	$\kappa$	$ϵ(\%)$
1	1804.18	4437	0.22	5.73	0.0079	3.63	25.45	72.73
2	592.20	2829	0.20	8.36	0.0071	6.11	23.66	12.53
3	324.35	702	0.67	1.57	0.0061	0.64	1.10	5.14
4	274.58	127	-39.99	67.21	0.2375	1.36	8.57	4.61
5	-32.78	90	1.44	1.00	0.0307	6.22	0.09	0.06

Equations

$$W(t)=t^n{e}^{-\omega \zeta t+i(\omega t+\phi )}H(t)$$

$$H(t)=\{\begin{array}{l}0,t<0\\ 1,t\ge 0\end{array}$$

$$M(t)=t^n{e}^{-\omega \zeta t}H(t)$$

$$M^{\prime }(\tau )$$

$$\tau$$

$$\bar{W}(t)$$

$$=\frac{M(t)}{M(\tau )}{e}^{i(\omega t+\phi )}$$

$$=t^{\zeta \omega \tau }{\tau }^{-\zeta \omega \tau }{e}^{\zeta \omega (\tau -t)+i(\omega t+\phi )}H(t)$$

$$w(t)=A{t}^{\zeta \omega \tau }{\tau }^{-\zeta \omega \tau }{e}^{\zeta \omega (\tau -t)+i(\omega t+\phi )}H(t)$$

$$r(t)=\sum _{i=1}{T}_i(w_i(t))$$

$$T_i(w_i(t,\tau ,A,\omega ,\zeta ,\phi ))=w_i(t-{\mathring{t}}_i,\tau -{\mathring{t}}_i,A,\omega ,\zeta ,\phi )$$

$$\phi =-(\tau \omega -\pi \lfloor \frac{\tau \omega }{\pi }\rfloor )$$

$$\kappa =\frac{\tau \omega }{2\pi }=\frac{\tau }{T}=\frac{\tau c_0}{T{c}_0}=\frac{d}{\lambda }$$

$$\hat{f}(\xi )$$

$$={\int }_{-\infty }^{\infty }A{t}^{\zeta \omega \tau }{\tau }^{-\zeta \omega \tau }{e}^{\zeta \omega (\tau -t)+i(\omega t+\phi )}H(t)e^{-i\xi t}dt$$

$$={\int }_0^{\infty }A{t}^{\zeta \omega \tau }{\tau }^{-\zeta \omega \tau }{e}^{\zeta \omega (\tau -t)+i(\omega t+\phi )}{e}^{-i\xi t}dt$$

$$=A{\tau }^{-\zeta \tau \omega }(\zeta \omega +i(\xi -\omega ){)}^{-\zeta \tau \omega -1}\Gamma (\zeta \tau \omega +1)e^{\zeta \tau \omega +i\phi }$$

$$|\hat{f}(\xi )|=A{e}^{\zeta \tau \omega }{\tau }^{-\zeta \tau \omega }({\zeta }^2{\omega }^2+(\xi -\omega {)}^2{)}^{\frac{1}{2}(-\zeta \tau \omega -1)}\Gamma (\zeta \tau \omega +1)$$

$$|\hat{f}(\omega +\xi )|=|\hat{f}(\omega -\xi )|$$

$$w_i(t)$$

$${\omega }_i$$

$$w(t)=A{e}^{\zeta \omega \tau (\ln t-\ln \tau )+\zeta \omega (\tau -t)+i(\omega t+\phi )}H(t)$$

$${ϵ}_{wi}=\frac{E_{wi}}{E_r}$$

$$E_{wi}=⟨w_i(t),w_i(t)⟩={\int }_{-\infty }^{\infty }|w_i(t){|}^{2}dt$$

$$\underset{{\mathbf{x}}_1}{\min }(E_{r2})=\underset{{\mathbf{x}}_1}{\min }{\int }_{-\infty }^{\infty }|r_1(t)-T_1(w_1(t)){|}^{2}dt$$

$${\mathbf{x}}_1=[A_1,{\omega }_1,{\mathring{t}}_1,{\tau }_1,{\zeta }_1,{\phi }_1{]}^{⊺}$$

$$\begin{array}{rl}r2& =r1-T_1(w_1)\\ r3& =r2-T_2(w_2)\\ & ⋮\\ r_{i+1}& =r_i-T_i(w_i)\end{array}$$

$$r_{i+1}$$

$${ϵ}_{wi}>{ϵ}_{w(i+1)}$$

$$\{w_i{\}}_E=\{w_i:i=1,2,3,\dots ,{\eta }_{90\%}\}$$

$$\{w_i{\}}_L=\{w_i:{\omega }_i<{\omega }_j,j=1,2,3,\dots ,i-1\},\text{for}i⩾({\eta }_{90\%}+1)$$

$$\{w_i{\}}_S=\{w_i{\}}_E\cup \{w_i{\}}_L$$

$$\hat{r}(t)=\sum _{i=1}{T}_i(w_i(t)),{\omega }_i(t)\in \{w_i{\}}_S$$

$$\left\{\begin{array}[]{lr}\frac{1}{n}\sum_{j=1}^{n}\omega_{j}=\bar{\omega}\\ \frac{|\omega_{j}-\bar{\omega}|}{\bar{\omega}}\ll 1,\quad j\in[1,n]\end{array}\right.$$

Full text

K****eywords Shock Waveform $\cdot$ Shock Decomposition $\cdot$ Pyroshock $\cdot$ Dominant Shock Distance $\cdot$ Shock Complexity

1 Introduction

Shock environment has to be considered for the design of equipment and components used in aerospace, transportation, defence and package engineering. Shock signals are normally recorded in the form of acceleration-time-history data, which is not straightforward to be linked to its severity and its effect on equipment components[1, p.385]. It is often necessary to get more specific mechanical information from shock signal using signal processing tools to understand the features of the shock and the associated phenomena.

Shock signal can be decomposed into a combination of basis signals. The most well-known decomposition method is Fourier transform, with which the frequency content of a measured signal can be identified by representing the signal with a series of harmonic waveforms. NASA proposed a decomposition method and an associated algorithm for shock signal using wavelet[2], with intention to replace the shock response spectrum (SRS) for shock specification. Prony decomposition method is usually adopted to express the shock as the superposition of exponentially damped harmonic waves in recent research study[3, 4, 5]. The Prony decomposition has advantage to describe near-field shocks characterised by the simultaneous peak location and sharp initial rise of its components. ESA provided an advanced scheme for Prony decomposition[1, 6], where the advanced Prony mode is given by the standard Prony mode convoluted with a basic Gaussian pulse. This scheme has advantage for the far-field shock environment characterised by the different peak time and the gradual initial rise of each advanced Prony mode. However, ESA did not give explicit mathematical expression for this scheme and its decomposition algorithm is not available. Similarly, the decomposed components of shock signals by empirical mode decomposition[7], discrete wavelet decomposition or wavelet package decomposition methods[8] are discrete signals without mathematical expressions.

Other basis waveforms can be used to synthesise mechanical shock environment for design and testing purposes. The WAVSYN waveform was initially proposed for the simulation of earthquakes, and was also found effective in shock synthesise[9, 10]. A waveform function similar to the response of a single-degree-of-freedom system to a Dirac impulse excitation was proposed in Ref.[11]. The ZERD waveform function was proposed to simulate shock signals with zero residual displacements, which has been well adopted with a shaker[12, 13]. A combination of simple impulses, sine waves, damped sine waves and modified Morlet wavelet has been used to synthesise shock time histories[14]. These methods are mainly for the synthesis of shock signal corresponding to a given SRS, rather than for the decomposition of a measured shock time history signal.

This study analyses the features of the generation mechanism of pyroshock and ballistic shock, and proposes a characteristic shock waveform to represent shock signals. This new shock waveform is determined by a group of parameters, i.e., amplitude, frequency, damping ratio, initial time, peak time and phase, which have explicit mechanics meaning. By decomposing the shock signal into the sum of limited terms of shock waveform components, mathematical expression of a shock signal can be explicitly given, which can support the understanding of shock effects.

2 Shock Waveform Expression from Shock Response

2.1 Shock Waveform

The shock waveform is proposed in the form of Eq.(1) based on the mechanical mechanism described in A.

[TABLE]

where $n$ is the order, $\omega$ is the angular frequency, $\zeta$ is the damping ratio, $i$ is the imaginary unit, $\varphi$ is the phase, and the $H(t)$ is the Heaviside step function defined by

[TABLE]

Step function is included because the magnitude of shock signal is zero before the start of shock event.

If $W(t)$ is considered as a modulated harmonic signal, the envelope of the shock waveform $W(t)$ is

[TABLE]

By calculating the stationary point of $M(t)$ in $t>0$, the time $\tau$ associated with the maximum of $M(t)$ is determined by

[TABLE]

The shock waveform proposed in Eq.(1) can be normalized in terms of the maximum value of its envelope, i.e.

[TABLE]

which can be generalised as

[TABLE]

where $A$ is the amplitude of the waveform.

With the relationship shown in Eq.(5), the parameter $n$, which has no direct mechanical meaning in Eq.(1), is now replaced by the parameter $\tau$ representing the peak time of a single shock waveform. $\tau$ is an important parameter when analysing the temporal structure of the shock waveform.

Considering a practical shock signal, the initial time for each shock waveform component may be different. Therefore, in the superposition process, the initial time for each shock waveform should be taken into consideration. Then, a shock signal $r(t)$ can be generally expressed as the superposition of a series of shock waveform components, i.e.

[TABLE]

where $T_{i}$ is a temporal translation transform for the consideration of the initial time for each waveform component, i.e.,

[TABLE]

and $\mathring{t}_{i}$ is the initial time of the $i$th shock waveform.

2.2 Mechanical Significance of Waveform Parameters

To better illustrate the shock waveform, the mechanical significances of waveform parameters in Eq.(9) are discussed below.

Frequency $\omega$

The shock waveform can be considered as an amplitude modulated harmonic wave. Here $\omega$ is actually the frequency of carrier wave, as depicted in Fig.1. Such a shock waveform will cause resonance for those modes whose natural frequencies are closed to $\omega$. Therefore, $\omega$ is one of the key parameters of the shock waveform.

Amplitude $A$

This is another key parameter of the shock waveform. The amplitude $A$ is the maximum of envelope curve, rather than the absolute maximum of a shock waveform $max|w(t)|$. As shown in Fig.2, usually the difference between maximum value of envelope and shock waveform are small. The condition of $A=max|w(t)|$ is that absolute maximum of carrier wave meets the maximum of its envelope, which is equivalent to

[TABLE]

where $\lfloor x\rfloor$ is the floor function with input $x$.

Damping ratio $\zeta$

This parameter mainly influences the duration of shock waveform. Smaller damping ratio leads to longer shock waveform duration, while larger damping ratio results in shorter duration, as illustrated in Fig.3. In the extreme case when $\zeta=0$, shock waveform is a harmonic wave.

Peak Time $\tau$

As shown in Fig.4, $\tau$ defines the temporal structure of each waveform component, i.e., the duration between the initial time and the peak time. If peak time aligns with initial time, i.e., $\tau=0$, shock waveform represent a Prony mode (ordinary damped harmonic wave). With the increase of $\tau$, the temporal structure of the waveform changes from Prony mode to the advanced Prony mode. When $\tau$ is large enough, the shape of shock waveform approaches symmetric wavelet. Their relationship will be further discussed in Section 2.3.

The physical meaning of peak time $\tau$ can be used to evaluate dominant shock distance, i.e., near-field, mid-field and far-field. The smaller the peak time, the closer the shock measurement point to the shock source. Considering the effect of frequency of carrier wave in time domain, it is proposed here that the dominant shock distance can be indicated by a parameter normalised by the period of the carrier harmonic wave, i.e.

[TABLE]

where $T$ is the period of the carrier harmonic wave; $d$ is the distance between the excitation source and the measurement point; $c_{0}$ and $\lambda$ are the phase velocity and the wavelength of the wave, respectively. Therefore, the parameter $\kappa$ is equivalent to the distance between the excitation source and the measurement point normalised by the wavelength of the wave.

In this study, a quantitative classification of the dominant shock distance is given in Table 1, based on the order of the magnitude of $\kappa$. With this classification, the waveforms of the near-field shock are similar to the damped sine waves, while the waveforms of the far-field shocks are close to the advanced Prony mode (or wavelet if $\zeta$ is high). Similar to the definition in Ref.[1], the transitional shocks with characters between near- and far-field shocks are classified as mid-field shock. Further analysis and discussion about $\kappa$ are given in Section 4.

Initial Time $\mathring{t}$

This parameter is necessary for the reconstruction of shock signal by the superposition of shock waveform component. The proposed shock waveform $w(t)$ has $t=0$ as the initial time, while it is possible that a shock signal is composed of shock waveform components having different non-zero initial times. Therefore, the simple translation transform in time domain can move the waveform component from [math] to $\mathring{t}$. The parameter $\tau$ denotes the duration between the initial time $\mathring{t}$ and the peak time of a translated shock waveform $T(w(t))$, so the translated peak time is $\mathring{t}+\tau$.

Phase $\varphi$

This is a trivial parameter since it does not offer any valuable mechanical information, but it is still necessary when a shock signal is decomposed numerically. This can be considered as a fine tuning of the temporal structure.

2.3 Connection between the Proposed Shock Waveform and Other Basis Signals

The shock waveform proposed in this study can cover a wide range of existing waveforms. Table 2 shows six special cases of the proposed shock waveform, and their corresponding conditions. If damping does not exist ($\zeta=0$), the proposed shock waveform is a harmonic wave with frequency $\omega$ and phase $\varphi$. The Prony mode (damped harmonic wave) is the same as the proposed shock waveform, if damping influence is included and the peak time is the initial time (i.e. $\tau=0$). Kern and Hayes proposed their waveform by multiplying Prony mode by time variable $t$[11], which is also a special case of shock waveform when $n=\tau\omega\zeta=1$. In far-field condition ($\kappa\gg 0$), the shock waveform looks like an asymmetric wavelet (e.g. B-spline wavelet), whose envelop is symmetric but not its phase. By further restricting the phase parameter of a far-field shock waveform with Eq.(12), a symmetric wavelet can be obtained, which is very close to the frequently used Morlet wavelet.

The shock waveform is generally similar to the advanced Prony mode proposed by ESA. A detailed description for the relationship between the proposed shock waveform and the advanced Prony mode is shown in B. It shall be noted that, the advance Prony mode denotes those waveforms similar to the far-field shocks in this study, although it can generate a wide range of existing waveforms.

2.4 Fourier Transform of Shock Waveform

The spectrum of the shock waveform in frequency domain can be determined by Fourier transform.

[TABLE]

where $\xi$ represents angular frequency and $\Gamma$ is the Gamma function. The amplitude of each frequency component is

[TABLE]

which satisfies

[TABLE]

Hence, the amplitude of frequency component is symmetric about frequency $\omega$, and reaches its peak when $\xi=\omega$. Shock waveform has a bell-like shape in frequency domain, as given by Eq.(18) and depicted in Fig.5.

According to the mechanical mechanism described in A, this study assumes that each bell spectrum in frequency domain of a shock signal ($r$) can be represented by a shock waveform $w_{i}$ component in time domain. If there exists a frequency band $[\xi_{i},\xi^{i}]$ for the $i$th bell, then the following relationship between the band-pass filtered signal, i.e., $r(t)|_{\xi_{i}}^{\xi^{i}}$, and a shock waveform is expected.

[TABLE]

and $\omega_{i},A_{i},\zeta_{i},\mathring{t}_{i}$ and $\tau_{i}$ are determined by both temporal and frequency structures of $r(t)|_{\xi_{i}}^{\xi^{i}}$; $\varphi_{i}$ can only fine tune the temporal structure in this case. Examples of Eq.(20) will be shown in Section 4.3.

3 Shock Decomposition Method

3.1 Program Implementation

A signal processing tool is proposed to decompose any shock signal into a sum of shock waveform components. To prevent potential arithmetic overflow during exponentiation calculation, Eq.(9) is rewritten in the following form,

[TABLE]

Generally, the main idea of shock decomposition algorithm as shown in Fig.6, is to fit a series of shock waveform components described by Eq.(22) to a shock signal. The energy ratio of a waveform component ($w_{i}$) is calculated in Eq.(23) to evaluate its weight of contribution,

[TABLE]

where $E_{wi}$ is the signal energy of $w_{i}$ calculated by

[TABLE]

and $E_{r}$ is the signal energy of shock signal $r$.

The first shock waveform component $T_{1}(w_{1})$ can be obtained by fitting the original shock signal $(r)$, i.e., $r_{1}=r$. The residue, $r_{2}$, comes from separating the first shock waveform component from the original shock signal by $r_{2}=r_{1}-T_{1}(w_{1})$, whose energy ratio $\epsilon_{r2}$ is still not negligible. This non-linear curve fitting problem

[TABLE]

is solved in least-squares sense with ‘lsqcurvefit’ algorithm in Matlab, where $\bm{x}_{1}$ is the solution vector of all parameters to be calculated.

[TABLE]

There may be many local minima of this non-linear problem, in contrast to linear curve fitting where the local minima is also the global minima. The global solution can be found by using various starting points ($\mathring{\bm{x}}_{1}$) for the fitting, which is achieved by ‘MultiStart’ strategy in Matlab. The lower bound (LB) and upper bound (UB) of fitting parameters can be set according to the mechanical condition of each shock signal, e.g., peak amplitude, interesting frequency range, etc. These starting points are discrete combination of parameters in the feasible region defined by LB and UB.

Normally $\epsilon_{w1}$ can account for the main proportion of the shock signal $r$, but the energy ratio $\epsilon_{r2}$ still accounts for non-negligible energy. Hence $r_{2}$ is treated as the new curve to be fitted. This non-linear fitting procedure can be repeated for all the subsequent residues $r_{i}$ ($i\geqslant 2$), with $r_{i+1}=r_{i}-T_{i}(w_{i})$.

[TABLE]

[TABLE]

According to Eqs.(27) and (28), $r_{i+1}$ is the global error between the shock signal and the reconstructed signal. The decomposition procedure finally stops when the energy ratio $\epsilon_{r(i+1)}$ of the global error $r_{i+1}$ is within the acceptable tolerance, which is 10% in this study.

3.2 Goodness of Fitting and Decomposition

The goodness of shock waveform decomposition is mainly depended on the set-up of ‘MultiStart’ strategy. The more starting points, the closer to the global solution (best fit). But having too many starting points can dramatically increase the computational time. Therefore, attention should be paid to the choice of feasible starting points.

Basically, starting points shall be decided based on the goodness required. In this study, starting points are chosen as follow:

• •

$A$: [$max(r)$];

• •

$\omega$: Every octave from $\omega_{LB}$ to $\omega_{UB}$ [$2^{0}\cdot\omega_{LB},2^{1}\cdot\omega_{LB},2^{2}\cdot\omega_{LB},\cdots,\omega_{UB}$];

• •

$\mathring{t}:[-\tau_{r1},\,0,\,\tau_{r1}]$;

• •

$\tau:[\tau_{r1}]$;

• •

$\zeta:[10^{-2},10^{-1},10^{0},10^{1}]$;

• •

$\varphi:[0,\,\frac{\pi}{2},\,\pi,\,\frac{3\pi}{2}]$.

Once the starting points are determined, the decomposition process and its results are repeatable.

The goodness of fitting and decomposition can be evaluated by the energy ratio ($\epsilon_{wi}$) of each shock waveform component $w_{i}$. If the solution in each fitting process is global solution, then following relation exists between the values of the energy ratio of consecutive shock waveform components according to Eq.(25).

[TABLE]

The starting points can be set-up based on practical need and engineering experience. However, it is recommended to increase starting points if $E_{wi}$ check is not satisfied.

3.3 Selection of Characteristic Shock Waveform Components

The decomposition process is implemented in least-squares sense, which usually concentrates on the waveform component with higher energy (usually at high frequency). However, only picking up high energy waveform components may hide low frequency components (usually with less energy), which are also important for the description of shock environment and the analysis of shock response. Such a problem can be solved below by the union of two waveform component sets.

Signal energy ratio of the residual signal $\epsilon_{ri}$ is very effective to evaluate the difference of two signals in time domain. In this study, the tolerance of energy ratio of the residual signal ($r_{i+1}$) is set to 10%, which is discarded for its limited influence. This means that the selected shock waveform components shall account for at least 90% energy of the shock signal. The least number of waveform components required to compose 90% energy of shock signal is defined as $\eta_{90\%}$, which is a parameter to show the complexity of a shock environment. A set $\{w_{i}\}_{E}$ can be used to include these waveform components, i.e.,

[TABLE]

Usually the waveform components set $\{w_{i}\}_{E}$ does not cover all important low frequency components. The method to select important low frequency component is given by

[TABLE]

where all qualified waveform components become the set $\{w_{i}\}_{L}$. Only waveform components with outstanding energy and lower frequency are included in this set.

The united set $\{w_{i}\}_{S}$ of $\{w_{i}\}_{E}$ and $\{w_{i}\}_{L}$ are a set of limited shock waveform components which can represent the time-frequency structure of a shock signal.

[TABLE]

The simplified shock signal $\hat{r}(t)$ can be reconstructed by waveform components in set $\{w_{i}\}_{S}$.

[TABLE]

4 Shock Signal Analyses

Shock waveform can analyse shocks originated from similar mechanical mechanisms, e.g., pyroshock, navy shock and ballistic shock, which are mainly determined by the free vibration of metallic structure (e.g., spacecraft, warship and armoured vehicle).

A near-field pyroshock, measured in the separation stage of an unnamed re-entry rocket vehicle[15] is taken as an example in this section. As a comparison, mid- and far-field shocks are produced by mechanical impact (metal-metal impact as shown in Fig. 7 (a)) in laboratory environment where a resonant plate and a sub-structure shown in Fig.7 (b) and (c), respectively, are used to produce acceleration shock signals for mid- and far-field shocks. The shocks are actually produced by the impact of a steel projectile (mass is 58g) and a fully fixed steel resonant plate with a sub-structure, at speed of 13.8m/s. The accelerometer for data collection is a KISTLER model K-shear 8742A20, whose acceleration and frequency measurement range are within $\pm$20,000g and 100kHz, respectively.

In the subsequent analysis, the near-field pyroshock from re-entry vehicle is termed ‘NFS’; the mid-field shock from resonant plate and the far-field shock from the sub-structure are referred to ‘MFS’ and ‘FFS’, respectively. It will be noted that there is no definite distinction between near- and mid-field shocks and between mid- and far-field shocks in common standards. Here MFS and FFS are considered a ‘farther’ shock than MFS and NFS, respectively, according to the characters described in ECSS Shock Handbook[1, p.33-35]. All the three shocks are shown in Fig.8.

4.1 Near-Field Pyroshock from Re-entry Vehicle (NFS)

Fig.9 shows the first ten waveform components of NFS from shock waveform decomposition, where each decomposed component of the shock signal is presented in the same scales of acceleration and time. The $\eta_{90\%}$ of NFS is 10, which means that 90% of NFS’s energy can be represented by the first 10 decomposed components according to the procedure in Section 3.1. The goodness of decomposition is acceptable by evaluating the relationship of the signal energy rate of each shock waveform component with Eq.(29). The parameters of twelve waveform components in set $\{w_{i}\}_{S}$ are listed in Table 3. The first ten components belong to set $\{w_{i}\}_{E}$ in energy sense, while the 11th and 12th components are low-frequency components selected by Eq.(31).

The first shock waveform component $w_{1}$ accounts for more than half of total signal energy, which is obviously a dominant component compared with the rest of waveform components. The characterisation of temporal structure of the shock signal $r$ can be outlined by $w_{1}$, whose maximum amplitude $A_{1}$ is 3939 $m/s^{2}$ at $\tau_{1}+\mathring{t}_{1}=0.68$ ms. The frequency $\omega_{1}$ is 2077 Hz, matching the knee frequency of the SRS or the dominant frequency component of the Fourier spectrum of the shock signal $r(t)$. Although the damping coefficient of the unnamed re-entry rocket vehicle is not given, the damping ratio $\zeta$ is 6.7% for $w_{1}$, close to the general estimation of damping ratio 5% of spacecraft structural[1].

In terms of the dominant shock distance, the overall shape of this near-field pyroshock ($r$) in Fig.10 is visually like the far-field shock rather than the near-field shock shown in Shock Handbook[1, p35]. Parameter $\kappa$ and $\epsilon$ for each waveform component are given on the right of the vertical line in Table 3. The value $\kappa_{1}=0.32$ corresponds to a near-field shock scenario and the shape of $w_{1}$ also belongs to an ordinary damped harmonic wave, which looks like measured near-field shock signal[1, p.403]. Because $\epsilon_{1}$ overpowers the rest of $\epsilon_{i}$, the dominant shock distance can be indicated by $\kappa_{1}$.

The shock decomposition method can simplify the description of a shock environment while keeping its main temporal structure and frequency content. The reconstructed waveform $\hat{r}=\sum_{i}\{w_{i}\}_{S}$ shall still represent the shock signal $r$ in both time and frequency domain. The time-history comparison of $r$, $w_{1}$ and reconstructed waveform $\hat{r}$ is illustrated in Fig.10. The 10% energy ratio of the residual signal $(r-\hat{r})$ guarantees the similarity between $r$ and $\hat{r}$ in time domain.

Comparisons of $\hat{r}$ and $r$ in frequency domain are presented in two forms, i.e., frequency spectrum in Fig.11(a) and SRS in Fig.11(b). The spectrum of $\hat{r}$ can cover several most significant bells of the spectrum of $r$ in frequency domain, depending on how many shock waveform components are involved. In this case, $\eta_{90\%}$ of NFS is ten, which means that ten bells are covered. The frequencies where maximum errors take place are the centre frequencies of components $w_{i}$ outside set $\{w_{i}\}_{E}$, whose energy ratios are less than 1%. These components are discarded since the cumulative error in the whole frequency region is less than the acceptable tolerance, i.e., 10%. The SRS of $\hat{r}$ also matches the SRS of $r$ very well. There exists certain error in low-frequency domain if shock waveform components in set $\{w_{i}\}_{L}$ are not included, which shows the importance of low-frequency waveform components for mechanical response.

4.2 Mid-Field Shock from Resonant Plate (MFS)

The first ten shock waveform components are shown in Fig.12. A more complicated mechanical scenario is indicated by $\eta_{90\%}$ of MFS having 25 and 3 components in $\{w_{i}\}_{E}$ and $\{w_{i}\}_{L}$, respectively. Although 28 waveform components are required in total to characterize this complex shock environment, only parameters for the first ten components in $\{w_{i}\}_{E}$ are listed in Table 4 for illustration purpose.

In the case of MFS, even though the first shock waveform component is not as outstanding as the first component of NFS, it can still provide some important features of MFS. The parameter $\kappa_{1}=6.12$ identifies MFS as a ‘farther’ shock compared to NFS, although they have similar overall shape visually. This coincides with the description of resonant plate technology[16], and hence, the MFS signal is regarded as a mid-field shock in this study. The damping ratio of the resonant plate apparatus is indicated by its $\zeta_{1}$ (5.7%), which coincides with the 5% general estimation of damping ratio as well.

The comparison of $r$ and $\hat{r}$ in both time and frequency domain are shown in Figs.13 and 14, respectively. The frequency spectrum may also explain the complexity of frequency component of MFS.

4.3 Far-Field Shock from Sub-Structure (FFS)

This far-field shock scenario is a strong support of the proposed shock waveform as a characteristic waveform. With only three waveform components, 90% signal energy of the shock signal can be represented. Time-history and parameters of $\{w_{i}\}_{S}$ are shown in Fig.15 and Table 5, respectively, where waveform components 4 and 5 belong to $\{w_{i}\}_{L}$.

The first shock waveform component $w_{1}$ of FFS is particularly important. It constitutes 72.7% of total energy of the shock signal. This shock environment is a common far-field shock, which can be validated both from its temporal shape and the parameter $\kappa_{1}=25.46$. The damping ratio of the sub-structure is implied by $\zeta_{1}=0.79\%$, which is much smaller than the damping of the resonant plate. This difference may be mainly resulted from different mounting and boundary condition.

Validation of representativeness in time and frequency domain between simplified shock $\hat{r}$ and original shock signal $r$ are depicted in Figs.16 and 17, respectively. The signal $\hat{r}$ reconstructed from only 5 shock waveform components can keep both important temporal and frequency features of the original shock signal.

FFS is an appropriate example of Eq.(20) for its frequency structure. The band of each bell $[\xi_{i},\xi^{i}]$ can be easily identified in the FFT spectrum because the bells are clearly separated from each other ($\xi^{i}\leqslant\xi_{i+1}$), which are not satisfied in the case of NFS and MFS in this study. Fig.18 compares the decomposed shock waveform components and their corresponding filtered signals, which are extracted by the ideal band-pass filter in Matlab. The pass-bands are shown in their respective legends. It is evident that both Eqs.(20) and (21) are satisfied in this FFS case. Band-pass filtered signals are almost the same as their corresponding shock waveform components. Lower bound $\xi_{i}$ and upper bound $\xi^{i}$ of filtering band are symmetry about the centre frequency $w_{i}$.

5 Conclusive Remarks

According to the shock generation mechanism, this study proposes an explicit mathematical expression of the basic shock waveform, with which the time-frequency structure of a shock environment can be characterised, and the links between time and frequency domains can be identified. Several index parameters with clear mechanical meaning are introduced, e.g., $\kappa$ and $\eta_{90\%}$, which can quantitatively describe the dominant shock distance and shock complexity. By retaining only important waveform components, a shock environment can be simplified while still keeping its main mechanical properties in both time and frequency domain. By comparing specific parameters that are not easily available from other signal processing tools, the proposed method can be used to compare different shock environment and track the change of characteristics of the shock signal at different location when it propagates in a structure. With the mathematical expression of shock signal, more convenient and reliable tool can be developed to characterise shock environment and its effects on equipment and devices.

Acknowledgements

The authors would like to thank William Fitzpatrick, Dr. Stephen Burley and William Storey for the conduction and support of the experiment of mechanical impact.

Appendix A Mechanical Basis for the Proposed Shock Waveform

Most of the shocks can be characterised by the free vibration of a structure subjected to a local impulsive loading. The loading duration in most shock scenarios is short, comparing to the periods of the dominant vibration modes of the structure, which transmit the shock from the loading location to the concerned location. This implies that the shock signal at the concerned location carries the contributions from the modes with high order natural frequencies of the structure. If a structure is considered as a MDOF system, the shock response of such system is very complex, although the response of each modal is standard damped harmonic oscillation. This is mainly because the modal density in high frequency domain (normally one order higher than the fundamental natural frequency) of target structure is very high, and the modal responses interact with each other.

It is assumed that a shock signal can be decomposed into a series of waveforms, and each of which can be represented by the response of a MDOF substructure subjected to the impulse loading. The frequency band of each MDOF substructure is continuous and narrow. Fig.1 shows a diagram of $n$-mass system representing a substructure. As shown in Eq.(34), the nature frequencies of the $j$th and $k$th modes are expressed as $\omega_{j}$ and $\omega_{k}$, respectively, which are assumed to be different but close to each other with a mean value $\bar{\omega}$. The damping ratio of these $n$ modes are all assumed to be $\zeta$.

[TABLE]

The dynamic motion equation of free vibration is

[TABLE]

where $\bm{m}$ is the mass matrix, $\bm{c}$ is the damping matrix, $\bm{k}$ is the stiffness matrix, and $\bm{v}$ is the displacement vector.

A solution in vector form for this motion equation is

[TABLE]

where $\bm{A}$ is the solution matrix and $\bm{u}_{0}$ is the initial displacement vector.

The solution can be expressed with its eigenvector matrix $\Phi$.

[TABLE]

where $\bm{\Lambda}_{A}$ is the eigenvalue matrix of matrix $A$. According to the mechanical assumption in Eq.(34), the diagonal elements of $\bm{\Lambda}_{A}$ are $\lambda_{j}=-\omega_{j}\zeta+i\omega_{j}$, whose mean value is $\bar{\lambda}=-\bar{\omega}\zeta+i\bar{\omega}$.

Note that the eigenvalues are centralized around the mean $\bar{\lambda}$, the solution vector $\bm{v}$ can be written as

[TABLE]

where

[TABLE]

According to the property of eigenvalue method, the eigenvalue matrix of $\bm{B}$ is

[TABLE]

Since the eigenvalue $\lambda_{j}$ in matrix $\bm{\Lambda}_{A}$ is close to $\bar{\lambda}$, the eigenvalue $\lambda^{B}_{j}$ in matrix $\bm{\Lambda}_{B}$ can be considered as small value.

The matrix $e^{\bm{B}t}$ can be expanded by Taylor series, i.e.

[TABLE]

The complete solution is

[TABLE]

The $t$ is relevant to the dominant frequency, i.e., $\bar{\omega}$ in this case. In pyroshock, ballistic shock and navy shock scenarios, $t$ within shock duration is very small. Combined with the small value of $\lambda^{B}_{j}$ in matrix $\bm{\Lambda}_{B}$, the matrix series in Eq.(A) can converge fast. By this way, the equations to describe such a shock response can be reduced form $n$ equations to much fewer equations.

For instance, if only the $n$th term in Eq.(41) is outstanding, and only the response of the $m$th DOF is interested, then

[TABLE]

where constant $C$ equals to the $m$th term in vector $\frac{1}{n!}\bm{B}^{n}\bm{u}_{0}$.

Recall that

[TABLE]

so the general waveform expression is in the form of

[TABLE]

Appendix B Discussion of Advanced Prony Mode

The advanced Prony mode is obtained by convoluting the standard Prony mode $P(t)$ with a basic Gaussian pulse $G(t)$ according to Refs.[1, 6]. This section attempts to derive the explicit mathematical expression of advance Prony mode, which was not given in Refs.[1, 6].

The following form of $P(t)$ and $G(t)$ are used

[TABLE]

where $A$ is the amplitude, $\zeta$ is the damping ratio, $\omega$ is the angular frequency, $\varphi$ is the phase, $\tau$ is the peak time of Gaussian pulse, and $\sigma$ is a parameter relative to the width of Gaussian pulse.

The advanced Prony mode $P_{A}(t)$ can be obtain with Mathematica

[TABLE]

where $\text{erf}(x)$ is the Gauss error function.

The number of parameters in the explicit mathematical expression Eq.(55) are the same as the proposed shock waveform in Eq.(9), if initial time $\mathring{t}$ is taken into account. Although the expression of Eq.(55) is very different from the proposed shock waveform in Eq.(9), the advanced Prony mode can also cover a wide range of existing waveform, e.g., the waveforms listed in Table 2.

However, the advanced Prony mode is not as convenient as the proposed shock waveform in Eq.(9) in terms of practical application. It is difficult to get an intuitive understanding of mechanical meaning of each parameter. The normalisation process of Eq.(55) cannot be implemented, since the absolute maximum of this expression is difficult to be obtained. The complexity mainly comes from the existence of time variable $t$ in the Gauss error function, which is a special (non-elementary) function. Advanced Prony mode may even include some waveforms which are uncommon for mechanical shock. Two such examples of waveforms produced by Eq.(55) are shown in Fig.(2) where the ‘zero velocity change’ validation criteria[1, p.417] is not satisfied.