Application of wave propagation to pyroshock analysis

Pyroshocks are high frequency transients due to pyrotechnic devices used in aerospace engineering in order to deploy solar arrays and antennae, separate subsystems from the spacecraft or separate the spacecraft itself from the base stage booster; their prediction is usually complex and very time consuming. The aim of this article is to investigate on the application of the analysis of transmission of waves inelastic media in order to predict the dynamic response to pyroshocks. The work is completed by numerical examples, related to components of common use in the aerospace engineering field, showing the comparison between results obtained by using both MSC-NASTRAN and this novel application of wavepropagation analysis.


Introduction
Launch vehicles make use of explosive charges to separate structural subsystems (such as a satellite) from the launcher, deploy appendices, activate or deactivate subsystems (such as valves of the propulsion subsystem).These devices have the drawback of giving rise to high accelerations (called pyroshocks) and extremely short transients with a very broad-band frequency spectrum that may be highly dangerous for flight equipments such as electronic devices.
The evaluation of the dynamic response to pyroshocks is a hard task: in this case classical techniques based on the Finite Element Method (FEM) are not efficient, since at high frequency a huge number of elements is required to match the smallest characteristic wavelength of the problem.
In [1] some methods are described permitting to deal with pyroshocks, nevertheless at the moment none of them is well established.A technique that usually is employed is that one suggested by NASA in [2]: it allows to determine the acceleration at the mounting points of most critical components by using empirical Transmissibilities that depend basically on the type of structure under analysis (e.g., honeycomb, skin-frame, monocoque, etc.), on the distance between source and the receiver locations and on the presence of junctions.These data have been collected mainly during an experimental activity performed several years ago [3].Anyway, it must be stressed that the procedure does not takes into account a number of parameters, e.g., material and thickness of the panels, so that it often leads to inaccurate results.
It follows that alternative techniques are needed.Lee and Henricks in [4] proposed an analytical method for predicting the upper bound for an envelope of mechanical responses due to a shock.In [5] Ruotolo suggested to use the Spectral Element Method (SEM) proposed by Doyle in [6] and described also in [7] and [8], taking advantage of the efficiency of SEM to deal with high frequency transients.Nevertheless, even though the computational effort required by SEM is much lower than that required by FEM, the usage of this method during the initial stage of the design process is not advisable.Indeed, during this phase sensitivity analysis is performed with the aim of determining the best choice for materials, thickness etc., so that a very fast computational method is needed while a reduced accuracy of results can be accepted.It follows that SEM is the best method to be used during the following phases, when structure properties have been already fixed, while during the preliminary design a faster tool is required.
Waves propagation in elastic media is a subject extensively studied, with direct application on the prediction of structure borne sound transmission [9].This approach allows to determine the amplitude of waves transmitted from one component to the next one by passing through a junction with a given geometry.In a following step, the ratio among transmitted, reflected and incident power can be evaluated, providing the basic data (i.e., coupling loss factors) to be used in Statistical Energy Analysis calculations in the high frequency range.
The basic idea behind this article is to take advantage of the capability of wave propagation analysis of determining the amplitude of transmitted waves in order to define Transmissibility Functions to be used in the framework of pyroshock analysis; it follows that the same basic philosophy that is behind the procedure proposed by NASA in [2] is applied.
As a result, the aim of this article is to investigate on the application of wave propagation analysis in order to predict the dynamic response to pyroshocks.In the following sections, firstly the analysis of waves propagating through a structure is reviewed and, subsequently, comparisons between results related to a simplified structure of a satellite and obtained by using MSC-NASTRAN and this novel application of wave propagation analysis are described and discussed.

Wave propagation analysis
Structural vibrations can be regarded as the superposition of elastic waves propagating through the structure.In presence of structural discontinuities the waves are partially reflected and partially transmitted, this process has been referred to as "attenuation of structure borne-sound" by Cremer et al in their book [9].
Numerical procedures capable of predicting this process of reflection and transmission are based on the determination of the properties of joints, e.g., beam and plate junctions, in terms of ratios between transmitted and reflected power with respect to the incident power, and on the definition of a relation between the damping of the component under analysis and the dissipated power of waves propagating through this component.
Cremer et al. in [9] analyzed joints between a limited number of either beams or plates (the latter under the assumption of cylindrical bending), up to 4, providing the reader with the basic knowledge to analyze more complex junctions.Subsequently, Craven and Gibbs in [10] proposed a versatile method of analysis of sound waves generated at a junction of plates as a result of a wave incident on one of the plates.Langley and Heron in [11] substantially extended this method and proposed a generic plate/beam junction which consists of an arbitrary number of plates which are either coupled through a beam or directly coupled along a line.They solved the problem of determining the amplitude of both reflected and transmitted waves by taking advantage of a matrix approach, arriving to the formulation of a dynamic stiffness matrix of the junction; as a result, their work allows to deal with real structures such as those currently used in the aerospace engineering field.
The previous technique can be applied rigorously only to semi-infinite plates joined along one edge.In contrast, space structures subjected to pyroshocks are characterized by finite dimensions; the corresponding rigorous approach, leading to exact results, is described in [5].Nevertheless, during the preliminary design phase approximate results can be accepted if provided by very fast techniques permitting to perform in short times sensitivity analysis.
As a consequence, the approach described subsequently follows from investigations performed by Langley and Heron in [11] with the main, and important, difference that the effect of reflected waves is neglected.
This approch has the main advantage that the evaluation of Transmissibilities is extemely fast and is based on two considerations: first, during the explosion of pyrotechnic charges usually every equipment is subjected to high acceleration levels during the first instants of the response, i.e., before the effect of waves reflected by subsequent junctions can be felt.Second, complex engineering structures, such as in-service space structures, at high frequencies are known to behave like heavily damped mechanical systems, with the result that the effect of reflected waves can be neglected since it is small with respect to the direct one.

N y M y T y N yx
x y z Fig. 2. A panel with its local reference system and corresponding forces (as acting on the beam).

Wave transmission at junctions
As already stated, the analysis of wave transmission at junctions follows strictly that proposed by Langley and Heron in [11], so that only the basics are given here.
The geometry of the junction considered in this article is shown in Fig. 1, where an arbitrary number of semiinfinite panels are joined along a common line by means of a beam; the motion of the panels are governed by the Kirchoff theory, while, according to [11], the equation of motion of the beam is derived by taking into account Timoshenko theory, moreover all the structural elements are assumed to be isotropic.
It is assumed that the motion of panels is described by variables (u, v, w) along the local reference system (x, y, z), as shown in Fig. 2. The equations of motion for the j-th plate, neglecting rotary inertia, are: (1) represent the bending and in-plane stiffness of the j-th panel, respectively, while and m j = ρ j h j is the density per unit area of the material.Stress resultants are related to displacements according to: As underlined in [11], under the assumption that the incident wave on one panel has space/time dependency given by exp −ikx+µy+iωt , compatibility of displacements at the junction requires that all panels must have the dependency exp −ikx+iωt , leaving the dependency on y to be determined by the equation of motion of the panels, i.e., by their properties and by wave type.
By focusing the attention on the out-of-plane motion, the dynamic response of the j-th panel can be written as the superposition of flexural waves moving toward the free edge (at infinity), i.e., characterized by a negative (real or imaginary) wavenumber µ B : As a result, it follows: This expression allows to determine displacement, w e j , and rotation, θ e j = w ,y , of the panel edge at the junction, i.e., where y = 0.By inverting the relation between w e j , θ e j and W 1 , W 2 it is possible to express the out-of-plane motion of the panel as a function of the constraint edge degrees of freedom, w e j and θ e j .This last expression can be subsequently introduced into Eq.( 2) to provide a direct relation between shear, bending moment, displacement and rotation at the constraint edge: where As pointed out in [11] a similar procedure is applied in order to determine a relation between in-plane shear, N yx , traction, N y , and in-plane displacements at the constraint edge.
By recalling that wavenumbers for longitudinal and shear waves in the panels are given by: the y dependency of waves moving toward the free edge of the panel is: and in-plane motion is given by the superposition of two types of waves according to: As well as in the case of out-of-plane motion, this last equation evaluated at y = 0 allows to determine the constraint edge displacements u e j and v e j as a function of amplitudes P L and P S .By inverting this relation, in-plane

MSCNASTRAN results
Fig. 6.Comparison between accelerations predicted by MSC-NASTRAN (continuous line) and by wave propagation analysis (dashed line) at node 226.
variables u and v are expressed as a function of the constraint edge degrees of freedom; this last equation introduced into Eq.( 2) allows to determine the dynamic relation between in-plane forces and displacements at the constrained edge: (9) where L ) According to [11], Eqs ( 5) and ( 9) can be used to produce a relation, expressed in the local reference system, between four generalized forces and four generalized displacements at the constrained end of every panel: with In order to switch to the global reference system, the following rotation matrix must be taken into account: with t 1j = y j sin θ j − z j cos θ j t 2j = y j cos θ j + z j sin θ j allowing to take into account that the junction of plates with the beam is not located on the beam neutral axis as well as that the j-th panel has orientation given the angle θ j .As a result, calling {F } the external forces acting on the junctions, to be determined as described subsequently, and {a} = u v w θ T the displacements of the junction, for m panels the following equation holds: where [ K P ] is the dynamic stiffness matrix of the connected, semi-infinite panels.
The beam connecting the various panels can be taken into account in a straightforward way.Indeed by writing its equation of motion, neglecting any warping effect of the cross section as well as rotary inertia, it is possible to derive its dynamic stiffness matrix [ C ] relating forces {F } and degrees of freedom {a} already introduced.Elements of matrix [ C ] can be found in [11].
The presence of the beam is now considered by introducing its dynamic stiffness matrix into Eq.( 12) and recalling both equilibrium and compatibility conditions: External forces {F } acting on the junction are evaluated by considering the effect of the incident wave on a given panel, e.g., panel e.By focusing the attention on bending incident waves, caused by the explosion of pyrotechnic charges, with amplitude A, inclination φ on the junction, wavenumbers k = k B cos φ, µ = ik B sin φ along x and y directions, respectively, induced displacements, collected in {b e }, at the constrained edge (supposed to be free) are given by: and corresponding forces, collected in {F e }, acting on the junction are: In order to satisfy equilibrium requirements of the whole assembly of panels: As a result, for an incident bending wave of amplitude A = 1, after defining the frequency ω, wavenumbers k B , k and µ can be evaluated, allowing to evaluate the right hand side of Eq. ( 16).The dynamic matrix of the whole junction depends only on geometric and material properties, so that the dynamic response of the junction, {a}, can be quickly evaluated.Finally, the amplitude of reflected or transmitted bending waves of the j-th panel is obtained by extracting the w degree of freedom after projecting {a} on the local reference system through The ratio between the amplitude of the wave transmitted t(= b, l) to panel no.r and the wave incident i(= b, l) acting on panel no.s gives the Transmissibility of the junction, T ti r−s (ω).A Transmissibility function for transmitted bending waves is shown in Fig. 3 for a junctions with 3 beams of equal properties, with incident bending wave on the lower beam (beam no.1).

Damping effect
The decaying effect of the dynamic response in both space and time is associated with the presence of damping in the structure under investigation.It is well known that one way of describing the presence of damping is by using a complex modulus of elasticity, i.e., E = E(1 + iη), where η is the so-called loss factor [9].
As a result of using a complex modulus of elasticity, both wave propagation velocities and wavenumbers become complex.By considering a wave propagating in the x direction, i.e., exp iωt−ikx , and writing the wavenumber k = k r − ik i , the previous relation becomes: The first term represents a reduction of the wave amplitude and allows one to determine a Transmissibility Function taking into account damping effects.According to [9], imaginary parts of the wavenumber are given by: so that the amplitude attenuation due to damping at frequency ω and over a distance ∆x is given by: for longitudinal and bending waves respectively (the effect of transverse waves is neglected in the following section, so that it is not considered in Eq. ( 18)). Figure 4 shows the Transmissibility function T d for bending waves between nodes 30 and 73 of the structure shown in Fig. 5.

Numerical results
In order to validate the procedure described in the previous section, a finite element model of the beam structure shown in Fig. 5 has been developed.The exciting force has been located at the lower end of the structure; it is described by the following expression: Moreover, it has a duration of 1 msec, so that its overall shape is that of a pulse-sine.The FE model has been developed in order to assure that the dynamic response could be predicted accurately up to a frequency of 5 kHz, so that every beam, with the length of 0.5 m, of the structure has been divided into 50 elements.the modal transient solution has been run.Results have been collected at a number of 'control' nodes located in the centers of the beams, as shown in Fig. 5.
Wave propagation analysis has been applied as described in the previous section, by considering both the effect of damping and the effect of joints on the propagation of bending, longitudinal and shear waves.In particular, the out-of-plane acceleration at node no.30 has been used as a reference, permitting to evaluate the dynamic responses at nodes 73, 124, 175 and 226 through Transmissibilities due to both joint and damping effects.
In order to show the details of the procedure, the dynamic response at node 175 is evaluated once the out-of and in-plane acceleration at node 124 is known, denoted as A op 124 (ω) and A ip 124 (ω), respectively.The out-of-plane acceleration is related to a bending wave propagating between joints A and B while the in-plane acceleration is related to a longitudinal wave propagating between these nodes.As a result, the amplitude of bending and longitudinal waves at node B are obtained from those related at node 124 by applying a transmissibility taking into account the damping effect: ω) By considering the transmissibility due to the presence of joints, the out-of-plane acceleration (i.e., that important for pyroshock prediction) at node 175 can be derived: where the last term, T b d (ω), takes into account the dissipation of bending waves from joint B to node 175 due to damping effect.
Figure 6 compares results obtained by using MSC-NASTRAN and wave propagation analysis at control node 226.This figure shows a good agreement between reference results, i.e., those obtained by using MSC-NASTRAN, and wave propagation results, demonstrating that the proposed procedure can be used in order to determine dynamic responses at control nodes that are sufficiently far from the reference node, i.e., node 30, as well as from ends and discontinuities of the beams.This last aspect is extremely important and can be considered as the expected consequence of using a theory capable of describing the propagation of waves in semi-infinite structural media.Out-of-plane accelerations have been used to evaluate the Shock Response Spectrum1 in correspondence of the various control nodes.Results are summarized in Fig. 7 where curves obtained by considering accelerations predicted by MSC-NASTRAN and wave propagation are superimposed.The upper curve is related to the reference node, node 30, so that only MSC-NASTRAN results are available in this case, while the two lower couple of curves are related to nodes 73-124 and 175-226, respectively.Figure 7 shows that the wave propagation approach permits to predict with very good accuracy the SRS at nodes 175 and 226, i.e., very far from the source of the excitation, and with sufficient accuracy at nodes 73 and 124.
Finally, Table 1 summarizes the prevailing results in terms of residual SRS at every control node evaluated at about 1 kHz, i.e., where the SRS reaches a maximum value.
A similar analysis has been performed on the structure made of 5 connected plates shown in Fig. 8. Properties of the plates are the same of the properties of the beams used in the previous example.Moreover, plates are connected by using a beam whose cross-section is shown in Fig. 8 (dimensions in mm).The excitation has been applied on the lower edge of plate no.1, with time-history given by Eq. ( 19), and the dynamic response has been evaluated by using MSC-NASTRAN as described previously.'Control nodes' have been set at the middle of every plate so that, by considering the node on plate no.1 as the reference one, the dynamic response on the other four plates has been derived.Figure 9 compares dynamic responses provided by MSC-NASTRAN with those predicted by using wave propagation analysis on plate no.5, showing a good agreement.Finally, Fig. 10 shows Shock Response Spectra evaluated by using MSC-NASTRAN and wave propagation results.As well as in the previous example, the agreement is excellent for plates 4 and 5 while it is only sufficient for the two plates nearer to the source.

Conclusions
In this article a new approach aimed at pyroshock prediction has been proposed.It is based on the evaluation of Transmissibility functions permitting to consider the presence of both junctions and damping in structures, so that it makes rather straightforward the evaluation of pyroshock effects, given a reference signal, far from the exciting location.Numerical investigations illustrated in this study show that the proposed technique is capable of providing satisfactory results in general, and accurate results when the location under analysis is sufficiently far from the exciting source.Finally, it must be underlined that MSC-NASTRAN took several hours to evaluate the dynamic response of the structure under investigation, while approximated results can be determined by using the approach proposed in this article in few seconds.It follows that this new approach opens new possibilities mainly during the initial phase of the design process.

Fig. 1 .
Fig.1.The junction under analysis with the global reference system; position and orientation of a panel.

Fig. 3 .
Fig.3.Transmissibilities for bending waves, both incident and transmitted, between beams 2, 3 and 1 of a junction -beams of equal properties.

Fig. 7 .
Fig. 7. SRS predicted by MSC-NASTRAN (continuous line) and by wave propagation analysis (dashed line) at several nodes.

Table 1
SRS residual levels at high frequency for the various control nodes