- Email: [email protected]
one-dimensional wave model of stress wave propagation in a shock tube
Compurrrr& srrucrurrrVol. 2. No. 4. pp. 583-587.1986 F’rinted in GreatBritain.
0045-7949/86 53.00* .oO t I986Pergamon PressL:d.
A COUPLED FINITE ELEMENT/ONE-DIMENSIONAL WAVE MODEL OF STRESS WAVE PROPAGATION IN A SHOCK TUBE W. J. T. Mechanical Engineering Department,
DANIEL
University of Queensland, Australia
(Received
1 I December
St. Lucia 4067, Queensland,
1984
finite element model coupled to a numerical representation of one-dimensional wave propagation is used to model stress wave propagation in the structure of a large shock tunnel. Modelling
Abstract-A
difftculties and the method of coupling the two types of model are described. Results are presented which indicate that stress wave damper rods, attached to the tunnel, are effective as a means of controlling stress levels.
INTRODIJCTION
At the University of Queensland, a large shock tunnel has been designed to test theoretical principles aimed at developing a high speed ramjet engine to power space craft. The shock tunnel is illustrated in Fig. I and has a launch tube, 22 m in length, down which a piston is fired to compress helium in the compression tube. A diagraphm at the end of this tube (the “burster disk”) ruptures, releasing a shock wave which travels down the shock tube, enabling very high pressure levels to be attained momentarily at the test section. A pressure pulse peaking at 200 MPa occurs as the piston nears the end of the compression tube, as shown in Fig. 2. This leads to stress wave propagation in the structure of the shock tunnel. As severe stress levels are possible in the shock tube, two long rods known as stress wave dampers are attached to the shock tube by a yoke as shown in Fig. 1. These are intended to reduce stress levels in the shock tube downstream from the yoke. The structure can be modelled assuming one-dimensional wave propagation in the long rods and tubes. However, the joint between the compression tube and shock tube and the yoke region,. both shown in Fig. 3, are not likely to be modelled successfully in one dimension. The arrows on Fig. 3 indicate alternate paths by which a stress wave could travel between the compression tube and shock tube, demonstrating that the joint is not a simple change of section. A model of stress wave propagation combining two-dimensional finite element analysis for these two regions and a one-dimensional wave propagation model elsewhere, has been successfully used to predict transient stress levels in the shock tunnel. A program “WAVE” was written to implement the analysis.
583
ANALYSISOF STRESSWAVE PROPAGATIONWITH THE FINITE ELEMENT METHOD
The finite element method has to be applied to direct integration of transient dynamics problems with caution. The mesh tends to filter high frequencies from the solution[l] and changes in mesh refinement can cause spurious reflection[2]. Lumped mass and consistent mass matrices have been used, the choice depending on the method of time integration used. To avoid smoothing effects at the wave front, elements based on linear interpolation of displacement are preferred to quadratic ones. An example of variation in the wave front resulting in response to a step load is shown in Fig. 4. A row of rod elements is integrated with time steps of from half to twice the time the wave takes to transit one element. Smoothing of the wave front occurs if the time step is too large and small oscillations ahead of the wave front can occur if the time step is too small. The numerical response continues to oscillate about the solution obtained from simple wave theory. Oscillations of similar amptitude are expected from an exact elastic analysis[3], but the numerical oscillations are not of the correct periodicity. The best representation is obtained when the time step is matched to the time to transit one element[4]. Numerical problems are alleviated when examining the response to ramp loading, as high frequencies contribute less to the response. This is the situation for the present analysis, as the pulse of Fig. 2 can be approximated as a series of ramps. Reference [5] describes and compares numerical integration algorithms. The Newmark algorithm with integration parameters set to represent constant average acceleration has been used in the present work, as in Ref. 141. It is unconditionally
W. J. T.
DANIEL
burster launch
tube,
compression
tube,
disk
I shock
test section
tube,
I I 22m
L.67m
lm’
8m
Fig. I. Plan view of the shock tunnel
Fig. 2. Pressure
pulse exciting
the stress
facility.
wave
propagation.
effect of internal damping in steel and hence not unrealistic. High frequencies were progressively filtered from the pulse as it travelled.
ube
Fig. 3. Detail of the joint between the compression tube and the shock tube and of the yoke. stable, and leads to solutions showing little numerical decay. A test of propagation of a square pulse along a row of one-dimensional rod elements, with consistent mass matrices, showed a decay of e-0.00’r, where x is position along the rod. This decay of average stress values is comparable to the
COTTSCHALL ARRAYS Gottschall[6] has introduced a simple means of analysing one-dimensional wave propagation. Two arrays are used, one for the forward wave and one for the reverse wave. At each time step, stress values in the forward arrays are moved forward one element and those in the reverse arrays moved back one element. At any instant, the resultant stress is the sum of the values in the forward and reverse arrays. Reflection due to a fixed or free end, or reflection and transmission due to a change in cross section, is easily represented. In Fig. 5, a rod with a change of section is modelled, illustrating the technique. The method can also model material damping and friction[6]. This technique is more computationally efficient than direct integration of a finite element model of a rod or tube, as it involves very few calculations. It is ideal for most of the shock tunnel. However, to model the whole system, it is necessary either to (a) determine transmission and reflection coefficients from two- or three-dimensional models of the sections that cannot be modelled in one dimension or (b) couple the two models. The latter approach is preferable for accuracy but extravagant computationally. A Gottschall array model of the yoke region il-
.A coupled finite element’one-dimensional
Iz
wave model of stress wave propagation
rays simulate the time delays that occur while stress waves traverse the yoke. The second option can be implemented as illustrated in Fig. 7. Consider a forward wave with stress ul, propagating through the linear rod elements to node C and a reverse wave of stress u2 arriving from the Gottschall arrays. If the rod is of cross-sectional area A, a force Za2A applied to node C as shown, will cause a stress of cr2 in element a and -02 in element 6, according to wave theory. If g, is transmitted into both elements, then the average stress in elements a and b is uI , and this can be placed in the first element of the forward Gottschall array. The stress u, - u2 in element b is absorbed by the dash pot, which has an impedance matching that of the rod, leaving the reverse wave uz to transmit back down the rod. Note that element b is not part of the physical model, but serves to separate the forward and reverse waves. The averaging of stresses in elements a and b smooths the numerical oscillation of the finite element model effectively, but means that the stress wave gains about half a time increment on entering the Gottschall arrays. It can be seen from Fig. 4 that with time step matched to element length, there is a delay in the response to a load. For this reason, if, for instance a ramp is fed from the reverse Gottschall array into the row of rod elements, the ramp is transmitted with the correct slope, but at a time delay of about half of a time increment. As an arbitrary load can be regarded as the summation of a series of step loads, one occuring at each time interval, a similar delayed response will occur. These time effects are not very significant in a large model.
elemen!length
c= speed of sound a= tvne mterual
ckmtnt number x to a step load applied to a row of linear rod elements using varying time steps.
Fig. 4. Response
values A are multiplied by 11, Rl values B are mulhpiied by 12.R2
Fig. 5. Gottschall arrays representing a change in section.
Fig. 6. Gottschall array model of the yoke. forward
585
array
MODELS
OF THE SHOCK TUNNEL
The compression tube-shock tube junction
The junction of the compression tube and the shock tube was modelled with two axisymmetric meshes, a fine mesh shown in Fig. 8 and a coarser, QI reverse array regular mesh forming part of the model of the whole Fig. 7. A pair of Gottschall arrays attached to a row of tunnel, Fig. 9. Both used bilinear isoparametric firod elements. nite elements with consistent mass matrices. The pressure loading was separated into three compolustrating the first approach is shown in Fig. 6. The nents (diaphragm, end pressure, and radial pressure loading), as shown in Fig. 8. All three were found three reflection coefficients and six transmission to cause significant wave propagation out of the coefficients, along with time delays for transmission, can be determined by applying ramp loads to joint region. The time step for the mesh of Fig. 8 a finite element model of the yoke. The holding ar- could not be matched to element size. Hence, discompression
tube
IIll farce
I I I I II on
burster
II
I
I
_-
I _-
disk
-_-_-_ Fig. 8. Axisymmetric
finite element model of the joint between the compression tube.
tube and the shock
W. J. T. DANIEL
586
halt of yoke
Gottschall arrays are attached at A.6 and C. 1.mdlcates
nodes with the same displacement. Fig. 9. Model of the complete shock tunnel.
Table I. Maximum stresses and in the shock tunnel Location
Complete assembly
Launch tube Compression
tube
Shock tube (before yoke)
Shock tube (after yoke)
19 MPa at 102 ms - 19 MPa at 72 ms 22 MPa at I I I ms -22 MPa at 100 ms I64 MPa at 80 ms - 172 MPa at I I .5 ms II9 MPa at 96 ms -113MPaat28ms
persion of the wave and some spurous reflections occurred. However, a ramp pressure load, applied to this model and to the coarser mesh, gave very similar coefficients for transmission of axial stress into the compression tube and the shock tube. This indicated that the approximations of the geometry and of areas of contact, which were necessary to obtain a regular mesh, were justified. In the longer section of the compression tube represented in the coarse model, axisymmetric bending waves (or “bulges”) were evident, but were of an order of magnitude less than the axial plane waves. The yoke region
While wave propagation in the yoke should strictly be modelled in three dimensions, a threedimensional model would have slowed a solution of a model of the whole tunnel unacceptably. Hence models of the yoke, approximating the geometry, and assuming either plane strain or plane stress were used. Bilinear isoparametric elements were used, the mass matrices for the plane stress case including transverse inertia terms. There proved to be little difference between the two models, although there was a substantial difference from a one-dimensional solution, which had the effect of forcing equal transmission of stress into the damper rods and into the shock tube downstream. The yoke model in Fig. 9 has a higher mesh refinement than is necessary for the yoke, in order to maintain mesh regularity. The complete
No launch tube
38 -37 360 - 352 257 -241
IMPa at MPa at MPa at MPa at MPa at MPa at
ms ms ms ms ms ms
35 MPa - 30 MPa 25 MPa -28 MPa I44 MPa - 143 MPa 239 MPa -231 MPa
at 29 ms at 86.6 ms at 93 ms at 88 ms at 77 ms at 81 ms at 61 ms at I19 ms
elsewhere. The arrays representing the launch tube and part of the compression tube incorporate a change of cross section, as in Fig. 5. This model was varied to study the effects of removing the dampers or isolating the launch tube. In each case, a period of 120 ms or 12 500 iterations was analysed. The highest stresses occurring in the four parts of the tunnel, together with the times when they occurred, are listed in Table 1. As the short length of shock tube between the burster disk and the yoke tends to trap stress waves, the highest values are found there, except when the damper rods are removed. Their presence halves the maximum stress in the shock tube downstream. Isolation of the launch tube is seen to be unnecessary from the point of view of protecting the launch tube and detrimental to stress levels elsewhere.
CONCLUSION It has been shown that a two-dimensional finite element model can be coupled to a one-dimensional wave propagation model in analysing problems of transient dynamics. An alternative approach of using finite element models to determine transmission and reflection coefficients has also been discussed, as a means of achieving greater computational economy. The technique has enabled validation of the design of a large shock tunnel facility. REFERENCJL!S
shock tunnel
The model of the complete assembly shown in Fig. 9 uses rod elements to link the joint region and the yoke and Gottschall arrays to extend the model
I5 36 I6 30 41 29
No damper rods
1. T. Belytschko,
tegration.
N. Holmes and R. Mullen, Explicit inStability, solution properties, cost. ASME
AMD 14, I-21 (1983).
A coupled finite elementione-dimensional wave model of stress wave propagation 2. 2, P. 3azant. Spurous reflection of elastic waves in nonu~o~ finite element grids. Compur. ,~er~o~s A&. Mechnnics Engng 16, 91-100 (19781. 3. K. F. Graff, Wove Motion in Elastic Solids. Clarendon Press. Oxford f19751. 4. L. Guex and W. Janach, Elastic waves in a short cylinder, impacted by a thin bar. J. Sound Vib. 71, 3894% (1980f.
587
5. K. J. Bathe and E. L. Wilson, ‘VK~ericu~ Metkuds in Finite Efement Analysis. Prentice Hall, Englewood Clifis, N.J. (19761. 6. W. Gottschali, .J. study of high frequency harmonic response of belt amplifiers. Ph.D thesis, University of Queensland (19831.
0045-7949/86 53.00* .oO t I986Pergamon PressL:d.
A COUPLED FINITE ELEMENT/ONE-DIMENSIONAL WAVE MODEL OF STRESS WAVE PROPAGATION IN A SHOCK TUBE W. J. T. Mechanical Engineering Department,
DANIEL
University of Queensland, Australia
(Received
1 I December
St. Lucia 4067, Queensland,
1984
finite element model coupled to a numerical representation of one-dimensional wave propagation is used to model stress wave propagation in the structure of a large shock tunnel. Modelling
Abstract-A
difftculties and the method of coupling the two types of model are described. Results are presented which indicate that stress wave damper rods, attached to the tunnel, are effective as a means of controlling stress levels.
INTRODIJCTION
At the University of Queensland, a large shock tunnel has been designed to test theoretical principles aimed at developing a high speed ramjet engine to power space craft. The shock tunnel is illustrated in Fig. I and has a launch tube, 22 m in length, down which a piston is fired to compress helium in the compression tube. A diagraphm at the end of this tube (the “burster disk”) ruptures, releasing a shock wave which travels down the shock tube, enabling very high pressure levels to be attained momentarily at the test section. A pressure pulse peaking at 200 MPa occurs as the piston nears the end of the compression tube, as shown in Fig. 2. This leads to stress wave propagation in the structure of the shock tunnel. As severe stress levels are possible in the shock tube, two long rods known as stress wave dampers are attached to the shock tube by a yoke as shown in Fig. 1. These are intended to reduce stress levels in the shock tube downstream from the yoke. The structure can be modelled assuming one-dimensional wave propagation in the long rods and tubes. However, the joint between the compression tube and shock tube and the yoke region,. both shown in Fig. 3, are not likely to be modelled successfully in one dimension. The arrows on Fig. 3 indicate alternate paths by which a stress wave could travel between the compression tube and shock tube, demonstrating that the joint is not a simple change of section. A model of stress wave propagation combining two-dimensional finite element analysis for these two regions and a one-dimensional wave propagation model elsewhere, has been successfully used to predict transient stress levels in the shock tunnel. A program “WAVE” was written to implement the analysis.
583
ANALYSISOF STRESSWAVE PROPAGATIONWITH THE FINITE ELEMENT METHOD
The finite element method has to be applied to direct integration of transient dynamics problems with caution. The mesh tends to filter high frequencies from the solution[l] and changes in mesh refinement can cause spurious reflection[2]. Lumped mass and consistent mass matrices have been used, the choice depending on the method of time integration used. To avoid smoothing effects at the wave front, elements based on linear interpolation of displacement are preferred to quadratic ones. An example of variation in the wave front resulting in response to a step load is shown in Fig. 4. A row of rod elements is integrated with time steps of from half to twice the time the wave takes to transit one element. Smoothing of the wave front occurs if the time step is too large and small oscillations ahead of the wave front can occur if the time step is too small. The numerical response continues to oscillate about the solution obtained from simple wave theory. Oscillations of similar amptitude are expected from an exact elastic analysis[3], but the numerical oscillations are not of the correct periodicity. The best representation is obtained when the time step is matched to the time to transit one element[4]. Numerical problems are alleviated when examining the response to ramp loading, as high frequencies contribute less to the response. This is the situation for the present analysis, as the pulse of Fig. 2 can be approximated as a series of ramps. Reference [5] describes and compares numerical integration algorithms. The Newmark algorithm with integration parameters set to represent constant average acceleration has been used in the present work, as in Ref. 141. It is unconditionally
W. J. T.
DANIEL
burster launch
tube,
compression
tube,
disk
I shock
test section
tube,
I I 22m
L.67m
lm’
8m
Fig. I. Plan view of the shock tunnel
Fig. 2. Pressure
pulse exciting
the stress
facility.
wave
propagation.
effect of internal damping in steel and hence not unrealistic. High frequencies were progressively filtered from the pulse as it travelled.
ube
Fig. 3. Detail of the joint between the compression tube and the shock tube and of the yoke. stable, and leads to solutions showing little numerical decay. A test of propagation of a square pulse along a row of one-dimensional rod elements, with consistent mass matrices, showed a decay of e-0.00’r, where x is position along the rod. This decay of average stress values is comparable to the
COTTSCHALL ARRAYS Gottschall[6] has introduced a simple means of analysing one-dimensional wave propagation. Two arrays are used, one for the forward wave and one for the reverse wave. At each time step, stress values in the forward arrays are moved forward one element and those in the reverse arrays moved back one element. At any instant, the resultant stress is the sum of the values in the forward and reverse arrays. Reflection due to a fixed or free end, or reflection and transmission due to a change in cross section, is easily represented. In Fig. 5, a rod with a change of section is modelled, illustrating the technique. The method can also model material damping and friction[6]. This technique is more computationally efficient than direct integration of a finite element model of a rod or tube, as it involves very few calculations. It is ideal for most of the shock tunnel. However, to model the whole system, it is necessary either to (a) determine transmission and reflection coefficients from two- or three-dimensional models of the sections that cannot be modelled in one dimension or (b) couple the two models. The latter approach is preferable for accuracy but extravagant computationally. A Gottschall array model of the yoke region il-
.A coupled finite element’one-dimensional
Iz
wave model of stress wave propagation
rays simulate the time delays that occur while stress waves traverse the yoke. The second option can be implemented as illustrated in Fig. 7. Consider a forward wave with stress ul, propagating through the linear rod elements to node C and a reverse wave of stress u2 arriving from the Gottschall arrays. If the rod is of cross-sectional area A, a force Za2A applied to node C as shown, will cause a stress of cr2 in element a and -02 in element 6, according to wave theory. If g, is transmitted into both elements, then the average stress in elements a and b is uI , and this can be placed in the first element of the forward Gottschall array. The stress u, - u2 in element b is absorbed by the dash pot, which has an impedance matching that of the rod, leaving the reverse wave uz to transmit back down the rod. Note that element b is not part of the physical model, but serves to separate the forward and reverse waves. The averaging of stresses in elements a and b smooths the numerical oscillation of the finite element model effectively, but means that the stress wave gains about half a time increment on entering the Gottschall arrays. It can be seen from Fig. 4 that with time step matched to element length, there is a delay in the response to a load. For this reason, if, for instance a ramp is fed from the reverse Gottschall array into the row of rod elements, the ramp is transmitted with the correct slope, but at a time delay of about half of a time increment. As an arbitrary load can be regarded as the summation of a series of step loads, one occuring at each time interval, a similar delayed response will occur. These time effects are not very significant in a large model.
elemen!length
c= speed of sound a= tvne mterual
ckmtnt number x to a step load applied to a row of linear rod elements using varying time steps.
Fig. 4. Response
values A are multiplied by 11, Rl values B are mulhpiied by 12.R2
Fig. 5. Gottschall arrays representing a change in section.
Fig. 6. Gottschall array model of the yoke. forward
585
array
MODELS
OF THE SHOCK TUNNEL
The compression tube-shock tube junction
The junction of the compression tube and the shock tube was modelled with two axisymmetric meshes, a fine mesh shown in Fig. 8 and a coarser, QI reverse array regular mesh forming part of the model of the whole Fig. 7. A pair of Gottschall arrays attached to a row of tunnel, Fig. 9. Both used bilinear isoparametric firod elements. nite elements with consistent mass matrices. The pressure loading was separated into three compolustrating the first approach is shown in Fig. 6. The nents (diaphragm, end pressure, and radial pressure loading), as shown in Fig. 8. All three were found three reflection coefficients and six transmission to cause significant wave propagation out of the coefficients, along with time delays for transmission, can be determined by applying ramp loads to joint region. The time step for the mesh of Fig. 8 a finite element model of the yoke. The holding ar- could not be matched to element size. Hence, discompression
tube
IIll farce
I I I I II on
burster
II
I
I
_-
I _-
disk
-_-_-_ Fig. 8. Axisymmetric
finite element model of the joint between the compression tube.
tube and the shock
W. J. T. DANIEL
586
halt of yoke
Gottschall arrays are attached at A.6 and C. 1.mdlcates
nodes with the same displacement. Fig. 9. Model of the complete shock tunnel.
Table I. Maximum stresses and in the shock tunnel Location
Complete assembly
Launch tube Compression
tube
Shock tube (before yoke)
Shock tube (after yoke)
19 MPa at 102 ms - 19 MPa at 72 ms 22 MPa at I I I ms -22 MPa at 100 ms I64 MPa at 80 ms - 172 MPa at I I .5 ms II9 MPa at 96 ms -113MPaat28ms
persion of the wave and some spurous reflections occurred. However, a ramp pressure load, applied to this model and to the coarser mesh, gave very similar coefficients for transmission of axial stress into the compression tube and the shock tube. This indicated that the approximations of the geometry and of areas of contact, which were necessary to obtain a regular mesh, were justified. In the longer section of the compression tube represented in the coarse model, axisymmetric bending waves (or “bulges”) were evident, but were of an order of magnitude less than the axial plane waves. The yoke region
While wave propagation in the yoke should strictly be modelled in three dimensions, a threedimensional model would have slowed a solution of a model of the whole tunnel unacceptably. Hence models of the yoke, approximating the geometry, and assuming either plane strain or plane stress were used. Bilinear isoparametric elements were used, the mass matrices for the plane stress case including transverse inertia terms. There proved to be little difference between the two models, although there was a substantial difference from a one-dimensional solution, which had the effect of forcing equal transmission of stress into the damper rods and into the shock tube downstream. The yoke model in Fig. 9 has a higher mesh refinement than is necessary for the yoke, in order to maintain mesh regularity. The complete
No launch tube
38 -37 360 - 352 257 -241
IMPa at MPa at MPa at MPa at MPa at MPa at
ms ms ms ms ms ms
35 MPa - 30 MPa 25 MPa -28 MPa I44 MPa - 143 MPa 239 MPa -231 MPa
at 29 ms at 86.6 ms at 93 ms at 88 ms at 77 ms at 81 ms at 61 ms at I19 ms
elsewhere. The arrays representing the launch tube and part of the compression tube incorporate a change of cross section, as in Fig. 5. This model was varied to study the effects of removing the dampers or isolating the launch tube. In each case, a period of 120 ms or 12 500 iterations was analysed. The highest stresses occurring in the four parts of the tunnel, together with the times when they occurred, are listed in Table 1. As the short length of shock tube between the burster disk and the yoke tends to trap stress waves, the highest values are found there, except when the damper rods are removed. Their presence halves the maximum stress in the shock tube downstream. Isolation of the launch tube is seen to be unnecessary from the point of view of protecting the launch tube and detrimental to stress levels elsewhere.
CONCLUSION It has been shown that a two-dimensional finite element model can be coupled to a one-dimensional wave propagation model in analysing problems of transient dynamics. An alternative approach of using finite element models to determine transmission and reflection coefficients has also been discussed, as a means of achieving greater computational economy. The technique has enabled validation of the design of a large shock tunnel facility. REFERENCJL!S
shock tunnel
The model of the complete assembly shown in Fig. 9 uses rod elements to link the joint region and the yoke and Gottschall arrays to extend the model
I5 36 I6 30 41 29
No damper rods
1. T. Belytschko,
tegration.
N. Holmes and R. Mullen, Explicit inStability, solution properties, cost. ASME
AMD 14, I-21 (1983).
A coupled finite elementione-dimensional wave model of stress wave propagation 2. 2, P. 3azant. Spurous reflection of elastic waves in nonu~o~ finite element grids. Compur. ,~er~o~s A&. Mechnnics Engng 16, 91-100 (19781. 3. K. F. Graff, Wove Motion in Elastic Solids. Clarendon Press. Oxford f19751. 4. L. Guex and W. Janach, Elastic waves in a short cylinder, impacted by a thin bar. J. Sound Vib. 71, 3894% (1980f.
587
5. K. J. Bathe and E. L. Wilson, ‘VK~ericu~ Metkuds in Finite Efement Analysis. Prentice Hall, Englewood Clifis, N.J. (19761. 6. W. Gottschali, .J. study of high frequency harmonic response of belt amplifiers. Ph.D thesis, University of Queensland (19831.