UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Acoustical modeling of rooms with extended-reaction surfaces Wareing, Andrew 2000

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
831-ubc_2000-0607.pdf [ 8.57MB ]
Metadata
JSON: 831-1.0089679.json
JSON-LD: 831-1.0089679-ld.json
RDF/XML (Pretty): 831-1.0089679-rdf.xml
RDF/JSON: 831-1.0089679-rdf.json
Turtle: 831-1.0089679-turtle.txt
N-Triples: 831-1.0089679-rdf-ntriples.txt
Original Record: 831-1.0089679-source.json
Full Text
831-1.0089679-fulltext.txt
Citation
831-1.0089679.ris

Full Text

ACOUSTICAL MODELING OF ROOMS WITH EXTENDED-REACTION SURFACES by ANDREW WAREING B.A.Sc, The University of British Columbia, 1995 A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF APPLIED SCIENCE in THE FACULTY OF GRADUATE STUDIES (Department of Mechanical Engineering) We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA July 2000 © Andrew Wareing, 2000 In presenting t h i s t h e s i s i n p a r t i a l f u l f i l m e n t of the requirements f o r an advanced degree at the U n i v e r s i t y of B r i t i s h Columbia, I agree that the L i b r a r y s h a l l make i t f r e e l y a v a i l a b l e f o r reference and study. I f u r t h e r agree that permission f o r extensive copying of t h i s t h e s i s f o r s c h o l a r l y purposes may be granted by the head of my department or by h i s or her r e p r e s e n t a t i v e s . I t i s understood that copying or p u b l i c a t i o n of t h i s t h e s i s f o r f i n a n c i a l gain s h a l l not be allowed without my w r i t t e n permission. Department of A f £c/wa/lICoJ ' n gg~r1' The U n i v e r s i t y of B r i t i s h Columbia Vancouver, Canada Date Abstract The acoustical modeling of rooms has always been a great challenge, especially when efforts are made to incorporate acoustical phenomena that are complicated to model. Knowledge of the acoustical behaviour of room surfaces is fundamental to predicting the sound field in a room. Surfaces are classified, acoustically, as of either extended or local reaction. All known room-prediction models assume, whether implicitly or explicitly, that surfaces are of local reaction. How can extended-reaction surfaces be incorporated into room-prediction models? What is the significance to predicted steady-state sound pressures in rooms of surfaces modeled as of extended vs. local reaction? This thesis is a detailed account of research dedicated to answering these questions. The main research goal was to develop a computationally efficient room-prediction model that included phase, and was applicable to rooms with extended- and local-reaction surfaces. A literature review concluded that a combined beam-tracing model with phase, and a transfer-matrix approach to model the surfaces, was the best choice. The transfer-matrix model is applicable to extended- or local-reaction surfaces. These surfaces are modeled as multi-layers of fluid, solid or porous materials. Biot theory was used in the transfer-matrix formulation of the porous layer. The new model consisted of the transfer-matrix model integrated into the beam-tracing algorithm. The model is valid for specular reflection only, and calculations were performed in the frequency domain. Both models were validated, and applied to three different room configurations: a 3m x 3m x 3m small office; a 10m x 3m x 3m corridor; and a 10m x 10m x 3m small industrial workroom. The test surfaces consisted of a glass plate, double drywall panels, double steel panels, a carpeted floor and a suspended-acoustical ceiling. Two predictions were made for each test surface in each room configuration - one for extended and one for local reaction. Results showed significantly higher predicted sound-pressure levels in rooms with a suspended-acoustical ceiling modeled as extended reaction, at low frequency. Rooms with walls of double drywall panels modeled as an extended-reaction surface showed significantly higher predicted sound-pressure levels. Sound-pressure levels were shown to be nearly equal in rooms with a single glass panel, carpet, and fibre-glass modeled as an extended-or local-reaction surface. An analysis of the reflection coefficients of these surfaces was performed to explain the results. The results confirm that it can be important to model the extended-reaction nature of room surfaces in some cases. ii Table of Contents A b s t r a c t i i L i s t of T a b l e s • v i L i s t of F i g u r e s v i i A c k n o w l e d g e m e n t s x i C H A P T E R 1 I n t r o d u c t i o n 1 C H A P T E R 2 L i t e r a t u r e R e v i e w 12 2.1 R O O M PREDICTION M O D E L S 12 2.2 S O U N D PROPAGATION THROUGH L A Y E R E D MATERIALS 17 2.3 R E S E A R C H OBJECTIVES 20 C H A P T E R 3 T h eory a n d I m p l e m e n t a t i o n 22 3.1 S O U N D PROPAGATION THROUGH L A Y E R E D MEDIA 2 2 3 . 1 . 1 F l u i d l a y e r s 2 6 3 . 1 . 2 E l a s t i c - s o l i d l a y e r s 2 7 3 . 1 . 3 E l a s t i c - p o r o u s l a y e r s 2 8 3 . 1 . 4 M u l t i p l e l a y e r s o f f l u i d , s o l i d , o r e l a s t i c - p o r o u s 3 1 3 . 1 . 5 S u r f a c e i m p e d a n c e a n d r e f l e c t i o n c o e f f i c i e n t 3 6 3.2 B E A M TRACING 37 3 . 2 . 1 D e f i n i n g s u r f a c e g e o m e t r y 4 0 3 . 2 . 2 M o d e l i n g o f a n o m n i - d i r e c t i o n a l s o u r c e 4 1 3 . 2 . 3 F i n d i n g a r e f l e c t i n g p l a n e 4 5 3 . 2 . 4 F i n d i n g a n i m a g e s o u r c e 4 6 iii 3 . 2 . 5 A n g l e o f i n c i d e n c e a n d d i r e c t i o n v e c t o r s f o r t h e r e f l e c t e d b e a m 4 7 3 . 2 . 6 R e c e i v e r d e t e c t i o n 4 9 3 . 2 . 7 C a l c u l a t i n g t h e s t e a d y - s t a t e s o u n d - p r e s s u r e l e v e l 5 0 3.3 IMPLEMENTATION 51 3 . 3 . 1 T h e t r a n s f e r - m a t r i x p r o g r a m 5 1 3 . 3 . 2 T h e b e a m - t r a c i n g p r o g r a m 5 4 C H A P T E R 4 M o d e l V a l i d a t i o n 57 4.1 VALIDATION OF THE TRANSFER-MATRIX MODEL 58 4 . 1 . 1 S i n g l e p l a t e 5 9 4 . 1 . 2 D o u b l e p l a t e 6 4 4 . 1 . 3 P o r o u s M a t e r i a l 6 8 4.2 VALIDATION OF THE BEAM-TRACING MODEL 77 4 . 2 . 1 I m a g e - p h a s e m o d e l 7 8 4 . 2 . 2 N u m b e r o f b e a m s a n d r e f l e c t i o n o r d e r 7 9 4 . 2 . 3 C o m p a r i s o n o f t h e b e a m - t r a c i n g a n d i m a g e - p h a s e m o d e l s f o r s e l e c t e d r o o m s 8 2 C H A P T E R 5 M o d e l A p p l i c a t i o n s 87 5.1 TEST CONFIGURATIONS 88 5.2 RESULTS 92 5 . 2 . 1 S i n g l e g l a s s p l a t e i n R o o m s 1 , 2 , a n d 3 9 2 5 . 2 . 2 D o u b l e d r y w a l l p a n e l e d w a l l s i n R o o m s 1 a n d 2 9 6 5 . 2 . 3 D o u b l e s t e e l p a n e l e d w a l l s i n R o o m 3 9 9 5 . 2 . 4 C a r p e t e d f l o o r i n R o o m s 1 a n d 2 9 9 5 . 2 . 5 S u s p e n d e d a c o u s t i c a l t i l e c e i l i n g i n R o o m s 1 a n d 2 1 0 3 iv 5 . 2 . 6 D o u b l e s t e e l p a n e l e d c e i l i n g i n R o o m 3 1 0 6 5 . 2 . 7 F i b r e - g l a s s l a y e r w i t h c o n c r e t e b a c k i n g a s t h e c e i l i n g i n R o o m 3 1 0 6 5.3 INTERPRETATION OF RESULTS 109 5 . 3 . 1 S i n g l e g l a s s p l a t e 1 1 2 5 . 3 . 2 D o u b l e - d r y w a l l p a n e l s 1 1 2 5 . 3 . 3 C a r p e t e d floor 1 1 6 5 . 3 . 4 S u s p e n d e d - a c o u s t i c a l - t i l e c e i l i n g 1 1 8 5 . 3 . 5 F i b r e - g l a s s w i t h c o n c r e t e b a c k i n g 1 2 0 5.4 SUMMARY 122 C H A P T E R 6 Conclusions 123 6.1 FUTURE WORK AND IMPROVEMENTS 128 Nomenclature 131 Bibliography 134 APPENDIX A Interface and Transfer Matrices 137 A . l TRANSFER-MATRICES 139 A . 2 INTERFACE MATRICES 142 APPENDIX B Plate Theory Transfer Matrices 146 B . l SINGLE-PLATE BOUNDED BY AIR ON BOTH SIDES 146 B . 2 DOUBLE-PLATE WITH AIR-GAP, BOUNDED BY AIR ON BOTH SIDES 149 APPENDIX C Matlab Code 153 C . 1 M A T L A B M-FILES FOR THE TRANSFER-MATRIX MODEL 153 C.2 M A T L A B M-FILES FOR THE BEAM-TRACING PROGRAM 167 v List of Tables T A B L E 3.1: M A T E R I A L PARAMETERS APPLICABLE TO FLUID, SOLID, A N D POROUS L A Y E R TYPES. . .52 T A B L E 4 .1 : M A T E R I A L PROPERTIES FOR THE GLASS A N D D R Y W A L L PLATES 60 T A B L E 4.2: M A T E R I A L PROPERTIES OF THE POROUS L A Y E R S USED IN VALIDATION CASES C) A N D D) 60 T A B L E 4.3: R O O M CONFIGURATIONS FOR THE B E A M - T R A C I N G VALIDATION TESTS 78 T A B L E 4.4: A V E R A G E PERCENTAGE DIFFERENCE IN COMPARING PREDICTIONS B Y THE IMAGE-PHASE A N D B E A M - T R A C I N G MODELS 83 T A B L E 5.1: R O O M CONFIGURATIONS 9 0 T A B L E 5.2: M A T E R I A L PROPERTIES OF THE STEEL PLATE USED FOR THE W A L L S A N D CEILING IN R O O M 3 90 T A B L E 5.3: INCIDENT A N G L E S F R O M THE FIRST ORDER REFLECTED-SOUND PATHS FOR THE THREE ROOMS I l l T A B L E A . 1: B O U N D A R Y CONDITIONS FOR A N INTERFACE BETWEEN ADJACENT LAYERS IN COMBINATIONS OF FLUID, SOLID, A N D POROUS L A Y E R TYPES 145 vi List of Figures FIGURE 1.1: A N IMPULSE RESPONSE RECORDED IN THE AUDIENCE SECTION OF THE C H A N CENTRE CONCERT H A L L , AT U B C 6 FIGURE 1.2: T H E FREQUENCY RESPONSE (MAGNITUDE A N D PHASE A N G L E ) OF THE IMPULSE RESPONSE SHOWN IN FIGURE 1.1 7 FIGURE 1.3: T H E SOUND-PRESSURE RESPONSE IN THE TIME DOMAIN, A T A RECEIVER, AFTER THE SOURCE IS TURNED ON, A N D THE D E C A Y AFTER THE SOURCE IS TURNED OFF 7 FIGURE 1.4: A N INCIDENT SOUND W A V E IS REFLECTED, ABSORBED A N D TRANSMITTED AT A SURFACE 8 FIGURE 3.1: W A V E PROPAGATION THROUGH A SINGLE L A Y E R OF M A T E R I A L 23 FIGURE 3.2: W A V E PROPAGATION THROUGH MULTIPLE L A Y E R S OF M A T E R I A L 31 FIGURE 3.3: TRAJECTORY-LOOP A L G O R I T H M 3 9 FIGURE 3.4: G E O M E T R Y OF A T R I A N G U L A R - F A C E D B E A M 40 FIGURE 3.5: S U R F A C E DEFINITION 41 FIGURE 3.6: B R E A K D O W N OF A N ICOSAHEDRON F A C E ; THE ICOSAHEDRON FREQUENCY IS 3 43 FIGURE 3.7: SOURCE WITH 3 2 0 B E A M S (ICOSAHEDRON OF FREQUENCY 4) 44 FIGURE 3.8: D I S T A N C E F R O M S TO INTERSECTION POINT ON A P L A N E 47 FIGURE 3.9: DIRECTION VECTOR OF A REFLECTED R A Y 48 FIGURE 3.10: A RECEIVER POINT THAT IS DETECTED B Y A B E A M 50 FIGURE 3.11: T H E TRANSFER-MATRIX A L G O R I T H M 53 FIGURE 3.12: T H E B E A M - T R A C I N G A L G O R I T H M 55 FIGURE 4 .1 : R E F L E C T I O N COEFFICIENT OF A SINGLE GLASS PLATE, 6 M M THICK, PREDICTED WITH PLATE THEORY 61 FIGURE 4.2: PREDICTION OF THE REFLECTION COEFFICIENT OF A 6 M M SINGLE GLASS PLATE, AT N O R M A L INCIDENCE 62 FIGURE 4 .3 : PREDICTION OF T H E REFLECTION COEFFICIENT OF A 6 M M SINGLE GLASS PLATE, AT 4 5 ° 63 FIGURE 4.4: PREDICTION OF THE REFLECTION COEFFICIENT OF A 6 M M SINGLE GLASS PLATE, A T 7 5 ° 63 vii FIGURE 4.5: SURFACE IMPEDANCE OF A SINGLE GLASS PLATE, 6 M M THICK, PREDICTED WITH THE TRANSFER-MATRIX MODEL 64 FIGURE 4.6: REFLECTION COEFFICIENT OF A DOUBLE DRYWALL PANEL, PREDICTED WITH PLATE THEORY 65 FIGURE 4.7: PREDICTION OF THE REFLECTION COEFFICIENT OF A DOUBLE DRYWALL PLATE WITH AIR GAP, AT NORMAL INCIDENCE 66 FIGURE 4.8: PREDICTION OF THE REFLECTION COEFFICIENT OF A DOUBLE DRYWALL PLATE WITH AIR GAP, AT 45° ANGLE OF INCIDENCE 67 FIGURE 4.9: PREDICTION OF THE REFLECTION COEFFICIENT OF A DOUBLE DRYWALL PLATE WITH AIR-GAP, AT 75° ANGLE OF INCIDENCE 67 FIGURE 4.10: SURFACE IMPEDANCE OF A DOUBLE DRYWALL PANEL WITH AIR-GAP, PREDICTED WITH THE TRANSFER MATRIX MODEL 68 FIGURE 4.11: SURFACE CONFIGURATION OF A POROUS LAYER WITH RIGID BACKING 69 FIGURE 4.12: SURFACE IMPEDANCE OF A SINGLE LAYER OF GLASS-FIBRE ATNORMAL INCIDENCE. 71 FIGURE 4.13: SURFACE IMPEDANCE OF A SINGLE LAYER OF GLASS-FIBRE AT VARIOUS ANGLES OF INCIDENCE, USING THE TRANSFER-MATRIX MODEL 72 FIGURE 4.14: REFLECTION COEFFICIENT OF A SINGLE LAYER OF GLASS-FIBRE AT VARIOUS ANGLES OF INCIDENCE, USING THE TRANSFER-MATRIX MODEL 73 FIGURE 4.15: PREDICTION VS. MEASUREMENT OF THE SURFACE IMPEDANCE AT NORMAL INCIDENCE, OF A CARPETED FLOOR 74 FIGURE 4.16: PREDICTION OF THE SURFACE IMPEDANCE AT OBLIQUE INCIDENCE, FOR CARPETED FLOOR 75 FIGURE 4.17: PREDICTION OF THE REFLECTION COEFFICIENT AT OBLIQUE INCIDENCE, FOR CARPETED FLOOR 76 FIGURE 4.18: SOUND-PRESSURE LEVEL PREDICTED BY THE BEAM-TRACING MODEL WITH REFLECTION NUMBERS 1 - 2 5 80 FIGURE 4.19: THE AVERAGE PERCENTAGE DIFFERENCE BETWEEN STEADY-STATE SOUND PRESSURE-LEVEL PREDICTIONS MADE BY THE BEAM-TRACING AND IMAGE-PHASE MODELS : 82 viii FIGURE 4.20: C O M P A R I S O N OF STEADY-STATE SOUND-PRESSURE-LEVEL PREDICTIONS B Y THE IMAGE-PHASE A N D B E A M - T R A C I N G M O D E L , FOR R O O M 1 84 FIGURE 4 .21 : COMPARISON OF STEADY-STATE SOUND-PRESSURE L E V E L PREDICTIONS B Y THE IMAGE-PHASE A N D B E A M - T R A C I N G M O D E L , FOR R O O M 2 85 FIGURE 4.22: C O M P A R I S O N OF STEADY-STATE SOUND-PRESSURE L E V E L PREDICTIONS B Y THE IMAGE-PHASE A N D B E A M - T R A C I N G M O D E L , FOR R O O M 3 86 FIGURE 5.1: T H R E E R O O M CONFIGURATIONS USED IN THE APPLICATION OF T H E R O O M PREDICTION M O D E L TO STUDY THE EFFECT OF E X T E N D E D - VS. L O C A L - R E A C T I O N SURFACES ON THE PREDICTED STEADY-STATE SOUND-PRESSURE L E V E L 91 FIGURE 5.2: F R E Q U E N C Y RESPONSE OF THE STEADY-STATE SOUND-PRESSURE L E V E L PREDICTED IN R O O M 1 WITH A SINGLE GLASS PLATE IN THE Y = 0 P L A N E ; A L L OTHER SURFACES H A V E A N A V E R A G E ABSORPTION OF 0.1 93 FIGURE 5.3: F R E Q U E N C Y RESPONSE OF THE STEADY-STATE SOUND-PRESSURE L E V E L PREDICTED IN R O O M 2 WITH A SINGLE GLASS PLATE PM THE Y = 0 P L A N E ; A L L OTHER SURFACES H A V E A N A V E R A G E ABSORPTION OF 0.1 94 FIGURE 5.4: F R E Q U E N C Y RESPONSE OF THE STEADY-STATE SOUND-PRESSURE L E V E L PREDICTED IN R O O M 3 WITH A SINGLE GLASS PLATE IN THE Y = 0 P L A N E ; A L L OTHER SURFACES H A V E A N A V E R A G E ABSORPTION OF 0.1 95 FIGURE 5.5: F R E Q U E N C Y RESPONSE OF THE STEADY-STATE SOUND-PRESSURE L E V E L PREDICTED IN R O O M 1 WITH DOUBLE D R Y W A L L PANELS (WITH 100 M M AIR GAP) FOR THE W A L L S ; A L L OTHER SURFACES H A V E A N A V E R A G E ABSORPTION OF 0.1 97 FIGURE 5.6: F R E Q U E N C Y RESPONSE OF THE STEADY-STATE SOUND-PRESSURE L E V E L PREDICTED IN R O O M 2 WITH DOUBLE D R Y W A L L PANELS (WITH 100 M M AIR GAP) FOR THE W A L L S ; A L L OTHER SURFACES H A V E A N A V E R A G E ABSORPTION OF 0.1 98 FIGURE 5.7: F R E Q U E N C Y RESPONSE OF THE STEADY-STATE SOUND-PRESSURE L E V E L PREDICTED PN R O O M 3 WITH D O U B L E STEEL PANELS (WITH 100 M M AIR GAP) FOR THE W A L L S ; A L L OTHER SURFACES H A V E A N A V E R A G E ABSORPTION OF 0.1 100 FIGURE 5.8: F R E Q U E N C Y RESPONSE OF THE STEADY-STATE SOUND-PRESSURE L E V E L PREDICTED IN R O O M 1 WITH CARPETED FLOOR; A L L OTHER SURFACES H A V E A N A V E R A G E ABSORPTION OF 0.1 101 ix FIGURE 5.9: F R E Q U E N C Y RESPONSE OF THE STEADY-STATE SOUND-PRESSURE L E V E L PREDICTED IN R O O M 2 WITH CARPETED FLOOR; A L L OTHER SURFACES H A V E A N A V E R A G E ABSORPTION OF 0.1 102 FIGURE 5.10: F R E Q U E N C Y RESPONSE OF THE STEADY-STATE SOUND-PRESSURE L E V E L PREDICTED IN R O O M 1 WITH ACOUSTICAL TILE SUSPENDED 4 5 7 M M F R O M CONCRETE CEILING; A L L OTHER SURFACES H A V E A N A V E R A G E ABSORPTION OF 0.1 104 FIGURE 5.11: F R E Q U E N C Y RESPONSE OF THE STEADY-STATE SOUND-PRESSURE L E V E L PREDICTED IN R O O M 2 WITH ACOUSTICAL TILE SUSPENDED 4 5 7 M M F R O M CONCRETE CEILING; A L L OTHER SURFACES H A V E A N A V E R A G E ABSORPTION OF 0.1 105 FIGURE 5.12: F R E Q U E N C Y RESPONSE OF THE STEADY-STATE SOUND-PRESSURE L E V E L PREDICTED IN R O O M 3 WITH DOUBLE STEEL PANELS O N THE CEILING; A L L OTHER SURFACES H A V E A N A V E R A G E ABSORPTION OF 0.1 107 FIGURE 5.13: F R E Q U E N C Y RESPONSE OF THE STEADY-STATE SOUND-PRESSURE L E V E L PREDICTED IN R O O M 3 WITH 100 M M THICK FIBRE-GLASS L A Y E R LINING A CONCRETE CEILING; A L L OTHER SURFACES H A V E A N A V E R A G E ABSORPTION OF 0.1 108 FIGURE 5.14: D I A G R A M OF R O O M 2, IN THE X - Z P L A N E (TOP) A N D THE X - Y P L A N E (BOTTOM), TO C A L C U L A T E THE INCIDENT A N G L E S FROM THE FIRST-ORDER REFLECTED-SOUND PATH I l l FIGURE 5.15: R E F L E C T I O N COEFFICIENT OF THE SINGLE GLASS PLATE, AS A FUNCTION OF INCIDENT A N G L E 113 FIGURE 5.16: R E F L E C T I O N COEFFICIENT OF T H E D O U B L E - D R Y W A L L P A N E L , AS A FUNCTION OF INCIDENT A N G L E 114 FIGURE 5.17: R E F L E C T I O N COEFFICIENT OF THE CARPET, AS A FUNCTION OF INCIDENT A N G L E . . . 1 1 7 FIGURE 5.18: R E F L E C T I O N COEFFICIENT OF THE SUSPENDED-ACOUSTICAL TILE, AS A FUNCTION OF INCIDENT A N G L E 119 FIGURE 5.19: R E F L E C T I O N COEFFICIENT OF THE FIBRE-GLASS L A Y E R WITH A CONCRETE B A C K I N G , AS A FUNCTION OF INCIDENT A N G L E 121 FIGURE A . 1: L A Y E R D I A G R A M 138 x Acknowledgements First, I would like to express my deepest gratitude to my supervisor, Professor Murray Hodgson, Director of the Acoustics and Noise-Control Laboratory in the Department of Mechanical Engineering. I am very grateful to him for his support, and the time and effort he has committed as my advisor for the duration of this project. I must also acknowledge a very special person in my life, my girlfriend Cynthia Dourasoff. I cannot thank her enough for the motivation and support that she has given me over the time that we've been together. I also appreciate her time reviewing my work and the use of her computer to run my simulations. Finally, I would like to thank all my friends and family for their support throughout the course of my graduate studies. Andrew Wareing June, 2000. xi C H A P T E R 1 Introduction One motivation for engineering research is to develop new techniques or methods that predict the outcome of a design problem given a set of input variables. The outcome may be one or more calculated value(s) of a particular quantity that yield(s) performance criteria to aid in decision-making during the design stage. The inputs to the problem may be design parameters that, if adjusted according to the calculated results, can be optimized to produce the best possible solution to the problem. With the rapid development of computer technology over the last 15 years, the engineer has an increasing ability to solve very complex design problems that at one time were not feasible. Computer-aided design (CAD) involves simulations in the computer environment to produce a numerical solution that leads to answers to problems without building expensive prototypes or performing costly experiments. For the purpose of engineering research, computer simulations have become necessary for rigorous modeling of problems when a solution is not possible or practical using analytical methods. Today, the hardware and software of personal computers (PC) are inexpensive, but powerful enough to run simulations using complex techniques such as the finite-element method (FEM), boundary-element method (BEM), statistical analysis, and ray or beam tracing to solve complicated engineering problems both in practice and for research. Mechanical engineers have used these numerical methods in solid and fluid mechanics for quite some time. Computer modeling has also played an important role in both practice and research in the field of acoustics. Acoustical modeling with computers has been significant in the design of loudspeakers, mufflers, sound barriers, enclosures, noise-control materials and, more recently, rooms. The acoustical design of rooms has always been a great challenge, especially when attempts are made to 1 incorporate acoustical phenomena that are very complicated to model. Room-prediction models vary greatly in their complexity, accuracy and computation times. They take room-acoustical input parameters - such as the room surface properties - into account under varying assumptions. These varying assumptions potentially have a profound impact on the usefulness and applicability of the room-prediction model. The research presented here is a documented solution to the challenge of room-acoustical modeling, with particular attention to the room boundaries or surfaces. In August of 1995, Professor Murray Hodgson, Director of the Acoustics and Noise-Control Laboratory in the Department of Mechanical Engineering, proposed a Masters of Applied Science research project. The general objective of this project was to develop a computationally efficient acoustical room-prediction model that incorporates the sound/structure interaction of the surfaces. Detailed modeling of sound propagation in rooms, and the sound/structure interaction, has been dealt with extensively. However, there is little knowledge of the impact that accurate surface modeling has on predicting the sound field in rooms. Fundamental to acoustical analysis is the study of sound generation and its propagation through media. In a fluid, sound energy is emitted from a source and propagates as fluctuating pressure in the form of a compressional wave. The sound energy is proportional to the square of the pressure magnitude and is detected at a receiver. The propagation of sound in a medium is governed by the wave equation, expressed in a linearized form as: c1 dt2 The regions of compression and rarefaction in a sound wave correspond to pressure disturbances or acoustical pressures that differ from the ambient pressure of the fluid. In the case of room acoustics, the fluid is stationary air at standard temperature and pressure, 20 °C and 101.3 kPa, respectively. Since sound propagates as a wave, it is meaningless to express the acoustical pressure without noting the frequency. Sound waves are a combination of many pure tones, or pressure waves at single frequencies, that may be analyzed in the frequency domain by performing signal processing operations such as the Fast Fourier Transform (FFT). Since human hearing is sensitive to sound in the 20 - 20000 Hz range, it is common to work in frequency 2 intervals such as octave- or third-octave bands, rather than at single frequencies. Octave- and third-octave bands are identified by their center frequency - for example, 125 Hz, 250 Hz, 500 Hz, and 1000 Hz octave bands. The desired frequency response of the acoustical pressure, with The emitted pressure wave exists in different forms depending on the physical nature of the source. A planar source such as a vibrating plate or piston emits a plane wave, with planar wave fronts and pressure magnitude decaying only with absorption in the medium of propagation. A point source such as a human speaker or loudspeaker emits a spherical wave, with spherical wave fronts and pressure magnitude decaying inversely as the square of the distance from the source. In most cases of room-acoustical analysis, and in this research, point sources emitting spherical waves are considered. It should also be pointed out that a typical source is only omni-directional - with sound energy radiating equally in all directions - at low frequency, and is more directional at high frequency. The directivity of a real source is very difficult to model and, for consistency, this research was conducted assuming an omni-directional source at all frequency. A brief explanation of the quantification of sound energy is necessary to understand the analysis of any acoustical problem. Sound is perceived as the pressure fluctuations detected at a receiver, due to a pressure wave emitted from a source. The sound energy is proportional to the square of the pressure p, and is commonly expressed as the sound-pressure level, Lp, in decibels a harmonic time dependence Helmholtz equation: can be obtained by re-writing the wave equation as the A2p + k2p = 0 1.2 where k = — and co = 2nf . c (dB): ( 2 1 Lp = 10 • log , with pref= 20 //Pa. \Pref ) 1.3 3 The initial sound energy emitted by a source depends on the source strength or sound power. Given the sound power W in Watts, the sound-power level, Lw, is expressed in decibels as: 4, =io- log , with Wrer = 1012 W. 1.4 The sound wave that travels from source to receiver may encounter obstacles and travel by more than one path in different environments. The simplest situation is a free-field with no obstacles, where there is only one 'direct' path of sound propagation from source to receiver. A free-field environment can be approximated in an anechoic chamber, a room with extremely sound-absorptive walls. A slightly more complex situation is an outdoor environment with no obstacles, and a partially absorptive ground surface. In this case, two paths of sound propagation from the source contribute to the total energy detected at the receiver. One is the direct path from source to receiver as in the free-field case, and the other is the reflected path from source to ground to receiver. The inclusion of obstacles influences the source-receiver path either by scattering/diffraction or by absorption/reflection/transmission of the sound wave incident from the source. The much more complicated room environment, with no internal obstacles, is the focus of this research. A room is an enclosure containing air, as well as the source and receiver. The air in the room is bounded by partially absorptive surfaces (walls, floor and ceiling). Some or all of the sound energy incident on each surface is reflected and remains in the room. The sound energy detected at a receiver point is the total contribution of those of the many source-receiver paths of sound propagation. The source-receiver paths in a room consist of the direct path, paths due to reflection from one surface, paths due to reflections from two surfaces, and so on. The sound field at a point in a room is completely defined with knowledge of the source signal and the room impulse response at that point. The objective of room acoustics is to predict the room impulse response for given source and receiver positions. The impulse response is the sound pressure received at the receiver, in the time domain, due to an impulsive signal radiated by the source. An impulse response that was measured in a concert hall is shown in Figure 1.1. Information on the arrival of sound at this 4 receiver point is obtained from this response, indicated by the acoustical pressure vs. time. The time delay between t = 0 and the first sharp peak in the impulse response is the propagation time for the sound wave to travel directly from the source to the receiver. The subsequent sharp peaks in the impulse response signal indicate the arrival of sound following surface reflections. As time increases, the time interval between peaks of pressure contributions from surface reflections decreases. Once all of the sound energy has dissipated, the sound pressure approaches the ambient pressure of the room. To evaluate a room acoustically, it is necessary to know the contributions of pure tones to the total sound signal - the frequency response. The frequency response is obtained by an FFT of the impulse response. The result is a complex signal that is the magnitude and phase of the pure tones that contribute to the total sound signal. From the impulse response shown in Figure 1.1, the frequency response was obtained; its magnitude and phase are shown in Figure 1.2. The peaks with the largest magnitudes indicate the pure tones that contribute the most significantly to the total sound signal. With knowledge of the frequencies, or octave band(s), that show the most prominent sound contributions, appropriate decisions can be made on the acoustical treatment of the room. The response to a continuous sound source in a room is obtained from the source signal and the room impulse response. The sound-pressure response to a continuous source is obtained by convolution of the source signal and the impulse response in the time domain, or by multiplication of the source signal and the frequency response in the frequency domain. When the source stops emitting sound energy, the sound pressure at the receiver decays to the ambient pressure in the room. The pressure response in the time domain, due to a continuous source that is turned on then later turned off is represented in Figure 1.3. The sound pressure at a receiver builds up before reaching equilibrium (steady state). The steady-state sound pressure is maintained until the source is shut off. After that point, the sound pressure at the receiver decays to the ambient pressure in the room. The rate of sound-energy decay and the steady-state sound pressure in a room, are two important parameters to evaluate the room acoustically. The rate of sound-energy decay is quantified by the reverberation time, or T 6 0 . The reverberation time is the time over which the sound-pressure level decreases by 60 dB. For a given room and source-receiver position, the reverberation time is an important parameter in characterizing the quality of sound. The steady-5 0.04 r -1.03 -0.02 -Time (ms) Figure 1.1: An impulse response recorded in the audience section of the Chan Centre concert hall, at UBC. The delay between t = 0 and the first peak is the time for the sound wave to travel directly from the source to the receiver. Subsequent peaks indicate the sound contributions of reflected source-receiver paths. Note that the first-order reflections are more prominent and distinctive. Later in time, the sound-pressure contributions from the source-receiver paths of higher reflection order decay to the ambient pressure in the room. state pressure-level determines the loudness. The desired reverberation time and steady-state pressure level vary according to the purpose of the room - i.e. concert hall, classroom, theatre, office, or industrial workspace. Correct modeling of the source-receiver paths involving reflection requires an understanding of the energy and wave effects associated with reflection, absorption and transmission of sound at a surface. Sound energy that meets a boundary is divided by reflection at the surface, absorption by the surface, and/or transmission through the surface. The sound-energy distribution at a surface is illustrated in Figure 1.4, with rays illustrating sound-propagation paths. This distribution of energy is quantified as follows: yO + a + r = l , 0 < / 7 < l , 0 < a < l a n d O < r < l . 1.5 where p = energy reflection coefficient a = energy absorption coefficient T = energy transmission coefficient 6 Frequency (Hz) Figure 1 . 2 : The frequency response (magnitude and phase angle) of the impulse response shown in Figure 1 . 1 . The peaks with the largest magnitude indicate the pure tones that contribute most significantly to the total sound pressure at the receiver point. Pressure Source turned on Figure 1 .3 : The sound-pressure response in the time domain, at a receiver, after the source is turned on, and the decay after the source is turned off. 7 Wave effects resulting from a sound wave incident on a surface involve pressure magnitude and phase changes of the reflected and/or transmitted sound waves. Complex pressure reflection and transmission coefficients are defined, to account for magnitude and phase changes at a boundary, as follows: P P R = -^,andT = 1.6 P P. i i where 0 < R < 1 and 0 < T< 1, real or complex (to include phase). With P, = Pj + Pr a necessary boundary condition, the pressure reflection and transmission coefficients are related by: P.+P 7/ = -* r- = l + R 1.7 P, In wave-based models, it is more common to specify the normal surface impedance as the boundary condition for a room surface. Acoustically, impedance is defined as the ratio of sound pressure to particle velocity of the sound wave. For plane waves propagating through a medium, the characteristic impedance, Zc, is constant with Zc = pQ • c, for density po and wave speed c . For air at standard temperature and pressure, the characteristic impedance is Zc = 415 Pa-s/m or rayl. The normal surface impedance is defined, at a point just on the incident side of 8 the boundary (points in Figure 1.3), as the ratio of pressure to the normal component of particle velocity. The pressure is the sum of the incident and reflected pressures; the normal component of velocity is the component of velocity in the direction of the surface normal. The normal surface impedance is a complex quantity that represents the magnitude of the pressure/velocity ratio and the phase difference between pressure and velocity at a boundary. The surface parameters discussed thus far are valid for specular reflection, and not diffuse reflections. Specular reflection is shown in Figure 1.4 - the angle of incidence is equal to the angle of reflection. The angle of incidence is the angle between the surface normal and the vector of the propagation direction of the incident sound wave. Sometimes the reflected sound wave is not directed along a distinct path but is, rather, scattered randomly in all directions. Such reflections are 'diffuse' and usually occur later in the sound build-up stage. The quantification of diffuse reflections is complicated and beyond the scope of this research. A surface is assumed to be of either local or extended reaction, which has implications on the normal surface impedance of the boundary. The common and often-assumed locally-reacting surface is such that the normal surface impedance is independent of the angle of incidence, and depends only on frequency. This is exactly true when the transmitted wave propagates in a direction along the surface normal, for all angles of incidences. Local-reaction is a reasonable assumption for surfaces when the wave speed of sound in the surface material is much less than the speed of sound in air. For many common building materials this is not always true and the acoustical behaviour of these materials is described more realistically as extended reaction. Extended reaction implies that the normal surface impedance varies with the angle of incidence at the boundary, and with frequency. Locally-reacting surfaces are often assumed for simplicity, but the validity of this assumption is not generally known. This is especially true for room surfaces and regarding the impact it has on determining the impulse response of a room. The impulse response of a room can be obtained by experiment, but many prediction methods have been developed for design purposes. Prediction methods are generally classified as wave-based (including phase) or energy-based models. More recently, however, there has been development toward models that are a combination of both classes. Traditionally, wave-based methods involve solving the wave equation analytically, or numerically with FEM, BEM, or finite differences. Wave methods are applicable to low frequencies where modal separation is 9 noticeable, and generally assume specular reflection. Numerical methods are more attractive for rooms with complex geometries and boundary conditions. However, finite- and boundary-element methods require at least five or six nodes per wavelength and the computation run-times become excessive at higher frequencies. In the mid- to high-frequency range, modal densities are too high to consider wave effects practically and, therefore, energy approaches are preferred. Common energy-based methods such as ray or beam tracing and the method of images, ignore phase and estimate the sound energy at a receiver point in the room. Ray tracing models sound propagation as energy rays that start from the source and propagate throughout the room. The total sound energy at a receiver is estimated by recording the energy contributions of rays that penetrate a small volume surrounding the receiver point. Specular and diffuse reflections at the room surfaces can be modeled with ray tracing. However, one disadvantage is that the divergence of rays at high reflection orders results in rays not penetrating the receiver volume when they should. Beam tracing is similar to ray tracing except that beams are traced rather than rays, allowing a point receiver and requiring fewer beams than rays to estimate the sound energy detected at a receiver. The use of beams solves the ray-divergence problem, but other problems occur when a beam intersects more than one surface, and when the beam front becomes larger than a surface. The method of images is another energy-based model, whereby image sources are calculated for all surfaces to estimate all of the source-receiver paths. The reflected path from source to receiver is the same length as the path from image source to receiver. This is an efficient and accurate model, but is limited to specular reflection and rectangular room shapes. Statistical Energy Analysis (SAE) and radiosity are other energy-based methods not discussed here. There is one common assumption in all of the room-prediction models mentioned - locally reacting surfaces. Why has this assumption been made in all prediction models? What would be the difference in the predicted steady-state sound pressures and reverberation time in a room with extended-reaction surfaces? An attempt at answering these questions was made in this research. A computationally efficient wave-based model was required to predict the response of a room bounded by extended-reaction surfaces. The general goal of this research was to develop this model, and use it to study the importance on modeling room surfaces accurately as of extended reaction. The specific research objectives were determined by a review of literature pertaining to the acoustical properties of extended-reaction surfaces and current room-prediction models. The literature review and specific research objectives are discussed in the next chapter. 10 Theory and implementation of the new room-prediction model are described in the third chapter. The fourth chapter involves model validation and establishes test cases of real rooms and surfaces for the application of the new model. The fifth chapter contains the analysis performed to study the impact of extended vs. locally-reacting surfaces on predicting the room response. The research objectives are achieved with the work completed in chapter five. Finally, chapter six is a summary of the work performed and conclusions drawn from the application of the new room-prediction model. 11 C H A P T E R 2 Literature Review In the previous chapter, a brief discussion of room-acoustical concepts led to the statement of the general research objective. The objective was to develop a computationally efficient acoustical room-prediction model that takes into account extended-reaction surfaces to calculate the steady-state sound-pressure level at a single receiver point. It is now necessary to review the literature that has led to the development of this model. This literature included topics on recent developments in room-prediction models and predicting the acoustical properties of surfaces. The first topic is room-prediction models, and includes analytical and numerical wave, geometric, and hybrid methods. The next topic discussed is sound propagation through layered materials, including Biot theory that describes wave propagation through elastic-porous materials. The review of the literature pertaining to these topics allowed the detailed research objectives to be defined. 2.1 Room-Prediction Models Discussion of the prediction of sound-pressure levels in rooms (or enclosures) begins with 'wave-based' approaches (predicting pressure as a function of time and position, with phase). The basic approach involves solving the wave equation directly with the appropriate boundary conditions. Solutions to the wave equation are well known and have been presented in numerous publications. For example, Kuttruff1 describes the wave equation and its solution for an enclosed space (room). While this is the most fundamental analytical solution and yields the exact answer, the computation is only practical for rectangular rooms at low frequency. Finite-12 and boundary-element methods (FEM/BEM) can be used to solve the wave equation numerically, and work well for rooms with arbitrary geometry2'3'4. However, the use of FEM/BEM is limited to low frequency and small room volumes, since five or six elements are required per wavelength2. With this minimum node requirement, solutions to rooms with large volumes, in the high-frequency range, result in excessive computation times even with modern computers. A sub-region coupling approach has been developed, where FEM/BEM has been applied to large room spaces5'6. Finite- and boundary-element methods are practical, and are used for small cavities such as the driver/passenger compartment in a car, and for the design of mufflers. In contrast to wave-based models, geometric models apply to higher frequencies and are 'energy-based' (ignoring phase), since the sound field in a room is more diffuse and individual room modes are almost indistinguishable. Several geometric approaches are the method of images, ray tracing, beam tracing, and radiosity. The method of images, presented by Allen and Berkley , involves calculating all source-receiver paths by replacing reflecting surfaces with image sources and computing their sound-field contributions at the receiver. The image method is relatively efficient and approaches the exact solution of the wave equation for rooms with rigid, specularly-reflecting surfaces and rectangular shape. For simulations involving high reflection orders and rooms with arbitrary geometry (more than six surfaces that form a rectangular room), the number of image sources increases exponentially and the computation becomes inefficient. Many of the image sources are not necessary, since they represent invalid source-receiver paths. This method has been applied to rooms with arbitrary geometry to model real concert halls, but calculation times are still prohibitive8. The ray-tracing technique, described by Kulowski9, estimates the source-receiver paths by tracing a large number (10000 or more) of one-dimensional rays from the source to a small receiver volume. The receiver cannot be a point, as would be desirable, since one-dimensional rays would never encounter it. The use of a non-point receiver introduces errors due to spatial and temporal smearing. The advantages of ray tracing are its application to rooms with arbitrary geometry, but the disadvantages are its high-frequency limitation and inaccuracies due to ray divergence. Ray tracing traditionally ignores phase and is, therefore, restricted to higher frequency at which modal density is high enough to neglect interference effects. Phase was introduced in an approximate manner by De Geest10 so that the ray-tracing method could be 13 applied to lower frequencies. The ray-tracing model with phase was shown to agree well with a FEM/BEM prediction of the steady-state sound-pressure level in a room. The ray-divergence problem concerns the angle of separation between adjacent rays, leading to inaccuracies at high reflection orders even when using a large number of rays. For a given set of adjacent rays expected to contribute sound energy at the receiver, the ray separation may be such that none of the rays will penetrate the receiver volume. The total sound energy, therefore, will be under-represented at the receiver. Increasing the number of rays improves accuracy, but the computation times are excessive. A better solution to the ray divergence problem is to replace rays with beams. A beam-tracing algorithm11 was developed for computer-graphics applications but adopted by acousticians as an improvement to the ray-tracing technique. Beam tracing is similar to ray tracing computationally, but with much improved efficiency for similar accuracy. For a given error due to ray divergence, fewer beams than rays are required to estimate the source-receiver paths in a room. In addition, a receiver point is used as it can be detected by beams that are three-dimensional. With respect to room acoustics, there have been two main approaches to implementing beam tracing from the ray-tracing model. A beam-tracing model that used conical beams traced by their central-axis ray was introduced by van Maercke and Martin . A spherical source, with an even energy distribution over its surface, was approximated by overlapping beams with a Gaussian energy distribution across each beam face. A relationship between the numbers of beams, rays and the reflection order was established by van Maercke and Martin, such that the ray-divergence error was minimized. This relationship showed that fewer beams were needed and, hence, the computation was more efficient. However, the use of overlapping beams over-estimated the sound energy detected at the receiver point with the possibility of multiple beams illuminating the receiver point. This was compensated for statistically by applying a weighting function to the energy distribution for each beam, so that for all beams detected at a receiver point, the correct sound energy was recorded. The statistical process required for overlapping beams was eliminated and the computational efficiency improved, by using triangular-faced instead of circular-faced beams. Beam tracing with triangular-faced beams has been implemented by Farina13'14 and Lewers15. Geometrically, these beams are pyramids with the vertex at the source and bounded by three planes. The advantages of using pyramidal, triangular-faced beams are two-fold; the 14 energy distribution across each face is constant, and adjacent beams always share a common boundary plane. The latter advantage implies that inaccuracies due to ray divergence and the multiple-beam detection problems with conical beams are no longer an issue. For each beam that contains the receiver point in its traced volume, there will be a unique source-receiver path. Farina initialized beams in his model by approximating a spherical source with an octahedron, and subdivided each face into more triangles to increase the number of beams. A better approximation to a spherical source was shown by Lewers15, who used a subdivided icosahedron (20-sided polyhedron). In both models, each beam was traced along its centre ray. A reflecting plane for a beam was determined by the closest surface that intersected the centre ray of that beam. The tracing of a beam centre ray to find reflecting surfaces introduced a new problem - it was possible for a beam to intersect more than one surface (i.e. a corner in a room), but the entire beam was reflected by the one surface determined by the beam's centre-ray. This resulted in the incorrect amount of sound energy detected at the receiver point; either too little energy was detected at the receiver from potentially valid beam components that were missed, or too much energy was detected from invalid beam components. This problem was discussed by Lewers, who argued that the net result of this error would be small from the cancellation of excess and deficient energy contributions from beams that intersected multiple surfaces. The adaptive beam-tracing algorithm by Drumm and Lam1 6 addressed the problem of a beam intersecting multiple surfaces. An incident (parent) beam that encountered multiple surfaces was subdivided into several reflected (child) beams. The cross-sectional shape of a child beam from one of the reflecting surfaces, adapted to the area of intersection between the parent beam and the surface. For increasing reflection order, the total number of beams to trace increases exponentially; however Drumm and Lam showed that only 20 beams were initially required to achieve good results. The adaptive beam-tracing model was compared to Lewers' beam-tracing and traditional ray-tracing models with similar accuracy and increased computational efficiency. The adaptive beam-tracing model was, however, more difficult to implement than Lewers' beam-tracing model. The radiant-exchange method, or radiosity1'15'17, is used to model the sound field in a room with diffuse surface reflections. The surfaces are subdivided into patches and, along with the source and receiver points, modeled as nodes in a network of energy exchange. The amount of sound energy exchanged between the nodes is calculated to determine the sound field in the 15 room. Radiosity is effective in predicting the reverberant tail of the room impulse response, due to dominant diffuse reflections. However, radiosity is computationally expensive in acoustical applications because of the large number of patches required for sufficient accuracy over a wide range of frequency. All of the individual methods described thus far have their advantages and disadvantages; however, the combination of these methods has led to more successful and complete room-prediction models for acoustical applications. Various numerical geometric room-prediction models have been developed as hybrids of the method of images, ray tracing, beam tracing, and radiosity. One such method is the combined ray-tracing/image-source algorithm of Vorlander18. This method works by calculating possible image sources based on the path history of rays traced from source to receiver. This has been effective in reducing computation times, while maintaining the accuracy of the image method, particularly for rooms with arbitrary geometry. Another model, presented by Dalenback19, incorporates both specularly and diffusely reflecting surfaces into a cone-tracing method. Diffuse reflection is accounted for by a split-up of cones incident on diffusing surfaces. A similar model that includes surfaces that exhibit both specular and diffuse reflection has been developed by Lewers15. This model used a triangular beam-tracing technique (non-adaptive) and the radiosity method for diffuse reflections. The commonality among all of the room-prediction models discussed is the assumption of locally-reacting surfaces. This assumption is explicit in a wave-based model, as the boundary conditions are defined by a normal surface impedance that is independent of incident angle. The local-reaction assumption is implied in energy-based models, which define boundaries with an energy-absorption coefficient. This diffuse energy-absorption coefficient is obtained from the energy-absorption coefficient that is a function of incident angle. The angle-dependent absorption coefficients have been directly implemented into a ray-tracing model by Hodgson . The angle-dependent absorption coefficients used by Hodgson and others, however, were determined from the normal surface impedance that was independent of incident angle. Presently, there is no known room-prediction model that uses extended-reaction surfaces. 16 2.2 Sound propagation through layered materials The theory of sound propagation through layers of homogeneous fluid has been published in many sources. Kinsler et. al. described sound propagation through single and multiple fluid layers of finite thickness and infinite extent. Expressions are given for acoustical properties such as pressure- and energy-reflection coefficient, absorption coefficient, and the normal surface impedance on the inner boundary (incident side) of the layer(s). Sound propagation through layers of elastic solids is more complicated because there are both compressional and shear waves to consider. Wave propagation through a single solid layer, with finite thickness and infinite extent, has been described using a transfer matrix22. The single solid-layer transfer matrix relates normal and tangential velocities, and normal and shear stresses, on each boundary of the layer. The normal surface impedance is calculated with knowledge of the normal stress and particle velocity at the surface of the layer, and the reflection coefficient is obtained from the surface impedance. Calculating the acoustical properties of surfaces that consist of multiple layers of solids is possible by multiplying the transfer matrices for each layer to yield one transfer matrix for the entire composite-layered surface. Layers of elastic-porous materials, such as glass fibre and foams, are considered as a two-phase material with fluid and elastic-solid components. Modeling wave propagation through elastic-porous layers is complicated, with several models developed that vary in complexity and validity. Porous materials are commonly used in noise-control applications to absorb sound. The accurate prediction of the acoustical properties of porous sound-absorbing materials is of great importance to those who work in the field. Various models involve empirical relations, the porous layer as an equivalent fluid layer with a rigid frame (solid component), and a porous 23 material with an elastic frame. The well-known empirical relations of Delany and Bazley calculate the characteristic impedance and wave number for a limited range of fibrous materials. The fundamental material property used was the flow resistivity - the ratio of pressure drop to flow rate of air propagating through the layer of porous material. A more rigorous model was developed by Attenborough24 for rigid-frame porous materials, for outdoor sound-propagation applications. Attenborough developed theory to predict the bulk modulus and surface impedance of porous materials modeled as an equivalent fluid; it did not include the elasticity effects of the frame. However, this model was applicable to more types of porous materials than those 17 described by the Delany and Bazley equations. Further development of Attenborough's model 25 was done by Johnson , who introduced the concept of tortuosity. Tortuosity is the ratio of the average path-length of a pore to the thickness of a porous material. The rigid-frame model was further improved by the rigid-frame model of Allard26, for indoor sound-absorbing materials. The Attenborough-Allard model accounted for irregular pore geometry by the introduction of characteristic dimensions. The characteristic dimensions are used to correct the effective bulk modulus and density of air in the porous material, to match the experimentally measured values for porous materials with irregular pore geometries. Prediction of the normal surface impedance of various porous materials using this model was shown to agree well with experiment23. The Attenborough-Allard model, however, does not include the elasticity effects of the frame, which include wave-coupling effects that can significantly affect the surface impedance of the material. The theory of wave propagation through elastic-porous media was developed by Biot27 in the 1950's, primarily for the purpose of seismology. Biot theory was later adopted by acousticians to predict the acoustical properties of porous materials. The significance of Biot's work applied to acoustics is in the application of the equations of motion derived for waves propagating through elastic-porous media and the coupling of such waves. Biot theory has been implemented into a FEM/BEM model developed by Panneton and Atalla28, and into a transfer-matrix method developed by Allard . Panneton and Atalla modeled surfaces as locally reacting, while Allard's transfer-matrix model included extended-reaction surfaces. The transfer-matrix model is simpler to implement, but only applicable to layers that are finite in thickness and infinite in extent. Nevertheless, transfer matrices have been proven an effective tool to evaluate the acoustical properties of materials, and used in this research. The application of Biot theory to elastic-porous sound-absorbing materials requires parameters to describe the elasticity of the frame and airflow through the pores. The frame is described by the Young's and shear moduli (for porous materials this usually is a complex quantity that includes a loss factor), and Poisson's ratio. These constants are not easily 30 31 determined for porous materials, though measurement techniques have been developed ' . Work has been done by Stinson and Champoux32 and Leclaire et. al.3 3 to determine the characteristic dimensions and shape factors for porous materials with specific pore geometries. Most of these parameters are not well known for many materials used as room surfaces and, therefore, Biot theory is not often applied in practice and remains in the research stage. With 18 improved measuring techniques and further interest, material parameters required for Biot theory should become available for common building materials. When the material data is available, it is expected that Biot theory will increasingly be used to predict acoustical properties of surfaces involving porous materials. The results presented by Allard34 indicate that there is an advantage to predicting the normal surface impedance and, hence, the reflection coefficient, with Biot theory. Room surfaces rarely consist of a single layer type; instead, they are comprised of multiple layers of different types of materials. An example of such a surface is a wall in a house, with a sheet of drywall, followed by an air-gap or insulating (porous) material, followed by, perhaps, a sheet of plywood that may have wood or aluminum siding on the outside. The prediction of acoustical properties, such as normal surface impedance and reflection coefficient, for individual layers of fluid, elastic-solid, and elastic-porous materials have been combined in the form of a transfer-matrix model35. This transfer-matrix model has been applied to a system of layers that are all finite in thickness but infinite in extent. Furthermore, it is applicable only to specular reflection. The transfer-matrix model was used in this research to model extended-reaction surfaces and implement Biot theory for porous materials. This model also has the potential to be directly implemented into a beam- or ray-tracing model that works with specular reflection only. The local- vs. extended-reaction assumption for a surface is a key point that is studied in this research. The review of literature on predicting acoustical properties of surfaces revealed very little information on this subject. The issue was discussed by Attenborough36 who pointed out that although many have argued the local-reaction assumption to be valid for porous materials, experimental measurement has shown otherwise. Experiments performed by Cops showed that the surface impedance measured for various porous materials varied significantly with the angle of incidence. The differences in the predicted steady-state sound fields resulting from local- vs. extended-reaction surfaces in a room, however, remain unknown. 19 2.3 Research objectives The specific research objectives were established upon review of literature on room-prediction models and sound propagation through layered media. These objectives, as reported in this thesis, are as follows: • Develop a transfer-matrix model similar to Brouard et. al.3 5, to predict the acoustical properties of extended-reaction multi-layered room surfaces that consist of layers of fluid, elastic-solid, and elastic-porous materials. Use Biot theory for the elastic-porous material. The last layer of a surface is backed by air or a rigid solid. • Develop a beam-tracing model based on Lewers'15 algorithm, and include phase, for predicting the steady-state sound-pressure level in rooms, with a single source and receiver point. • Implement a computer program that integrates the transfer-matrix and beam-tracing models. The transfer-matrix component is used to calculate a pressure-reflection coefficient of an extended-reaction room surface. This reflection coefficient determines the change in reflected pressure magnitude and phase for a reflecting beam. • Validate both models, separately, by comparing predictions with other theory and experimental measurements for selected test cases. The transfer-matrix and beam-tracing models are validated individually because there is no known combined model of this kind to compare with. • Analyze and draw conclusions on the prediction of the steady-state sound-pressure level in rooms with local vs. extended reaction surfaces. This is to be accomplished by applying the integrated transfer-matrix/beam-tracing model to a number of realistic room configurations. This thesis is organized to follow the above statement of research objectives. The development and implementation of the transfer-matrix and beam-tracing models are discussed in the next chapter. Chapter 4 is a detailed account of the validation process for both models. 20 The final and most important objective is covered in Chapter 5. A summary of the work performed in this research is given in Chapter 6, along with suggestions for future work. 21 C H A P T E R 3 Theory and Implementation In the previous chapter, literature dedicated to the development of room-prediction models and sound propagation through layered media was reviewed. The literature review established the specific research objectives, of which the first three are addressed in this chapter. The first section in this chapter presents the theory of sound propagation through layers of fluids, elastic-solids and elastic-porous materials. This theory provides the formulation of a transfer-matrix model, which is applied to predict the acoustical properties of local- and extended-reaction surfaces. The second section presents the details of an algorithm that was developed to implement a triangular-beam tracing model. The final section presents details of the implementation of the new room-prediction model - the transfer-matrix model integrated into the beam-tracing model - into a computer program. 3.1 Sound propagation through layered media The first requirement was to develop a transfer-matrix model to predict the complex surface impedance and pressure reflection coefficient of extended-reaction surfaces. This model was developed from the theory of sound propagation through a medium of finite thickness (along the surface normal), but infinite extent in the other directions. This medium is referred to as a layer and, for the purpose of this research, the layer may be a stationary fluid, an isotropic elastic-solid, or an isotropic porous-elastic material. The acoustical property to be predicted was 22 the complex pressure-reflection coefficient, to be used in the room-prediction model as the boundary condition for the surfaces. This reflection coefficient describes the pressure amplitude and phase of a reflected plane wave that is specularly reflected from a surface. Calculation of the reflection coefficient requires the normal surface impedance, determined from the normal component of particle velocity and the acoustical pressure just outside the boundary of the surface, on the incident-wave side. The configuration is illustrated in Figure 3.1; a plane wave with angle of incidence dis specularly reflected from the surface, on the outside of boundary A. If the surface has a reflection coefficient of less than unity, the sound wave is transmitted through the layer at an angle of incidence less than 9 if the layer is denser than the outer medium (as is the case for a room occupied by air). Another specular reflection of the transmitted wave occurs on the inside of boundary B, and a reflected wave propagates toward boundary A. Therefore, two waves propagate in a sound-absorbing layer, and two independent acoustical quantities are required to completely define the sound field for that layer. For each pair of transmitted and reflected waves in a layer, the normal component of particle velocity (the component in the direction of the surface normal) and the normal stress were chosen as the independent quantities at each boundary. The normal component of velocity V3 and pressure p, are required at the boundary A in order to calculate the normal surface impedance, Z A : \ 1 • A Nl N2 B e 9 y 1 k r r Figure 3.1: Wave propagation through a single layer of material. Incident wave is denoted by "i" , reflected wave denoted with "r", and transmitted wave indicated by "t". 23 v3(A) The normal velocity and stress components are determined from the properties inside the layer at boundaries Nj and N2 in Figure 3.1. The properties at A are related to the properties at Nj through boundary conditions which are such that continuity of displacements and a force balance exist. The situation is similar at the N2-B boundary; however the impedance at B is the characteristic impedance of air, with normal surface impedance Z B = ZQI COS(#). Therefore, the properties at N2 are known from those at B, and the properties at A are determined from those at Nj. The properties at Nj and N2, expressed as a vector V containing the normal velocity and stress components, are related via a transfer matrix35, T: V { N X ) = T - V ( N 2 ) 3.2 This transfer-matrix function completely describes the sound field within a layer given an incident and reflected sound wave(s) propagating through the layer. The number of independent acoustical quantities chosen at N\ and N2 is equal to the number of waves that propagate within the layer. A minimum of two quantities is chosen at a boundary, the normal component of velocity and stress, which are required to calculate the normal surface impedance. Calculation of the velocities and stresses at a boundary requires a transfer matrix that has been developed for the appropriate types of media - fluid, elastic-solid, and porous-elastic materials - with elasticity theory applied to deformable media. The theory presented here on elasticity and wave propagation through elastic-solid media has been reproduced from Allard38. The deformations are quantified by a displacement vector u, in terms of scalar and vector potentials, #>and \\> = (y/\, y/i), respectively: u = grad^ + curhj/ u - S7q> + V x \j/ With reference to Figure 3.1, the displacements are in the x\Xs plane; the direction of propagation is along the positive xi and X3 axes for the transmitted wave, and along the positive 24 xx and negative x3 axes for the reflected wave in the layer. Furthermore, the only non-zero component of the vector component i|/ is y/2, for displacements parallel to the xix3 plane. The potentials (p and y were obtained by Allard39, and used in this research: The exponential terms in Equation 3.4 contain the wave number components k\ and £ 3 . The component of the incident wave number vector in a direction tangential to the surface of a layer is k\, such that kx = |k| • sin(#) where k is the wave number vector of the incident wave (on the multi-layered surface) in free air and 0 is the incident angle. This component is equal to the k\ component of all sound waves that propagate in each layer. The component in the direction of the surface normal is denoted by £ 3 and is different in each layer. The acoustical influence on room surfaces or structures is such that the particle displacements in deformable media that make up the layers of the surface are due to the stresses in the material induced by the acoustical pressure in the room. The displacement u is the result of two actions, compression and shear. Compression displacements are described by the gradient (grade?) operation; displacements due to shear (rotations) are described by the vector curl (curly) operation of u in Equation 3.3. The stress-strain relations for isotropic elastic media are: <p = AeJ{a,~kl'l~k'Xl) + Bej{m'' y/2 =Cej{f"-k*-kiX>) + Dej{f 3.4 cr.. = A0So + 2fjev 3.5 where 0 = en + e22 + £33 = volumetric strain X and ju are Lame Coefficients, such that bulk modulus, K = A + 2/u/3, and X = vE (l + v X l - 2 v ) >• for an elastic solid, Poisson ratio v, and Young's modulus E. E 2(1 + v) 25 IJ [Oifz*;. The components of the strain tensor, ey, are defined in terms of the displacements w,: 1 ' • = 2 ' +• 1 v dXj dxi 3.6 The wave equations in a deformable media are determined from the stress equation of motion, expressed in terms of the displacement vector u: d2u i \d6 -> p—j- = (/t + ju) + juV2u, + Bt , where i = 1, 2, 3 3.7 dt dx. The wave equations in fluid, elastic-solid, and porous-elastic media, and the components of velocity and stress for such layers, are discussed further in the following sections. 3.1.1 Fluid layers Fluid layers consist of a stationary fluid. This implies that u — V<p, as Vxi|/ = 0 for no rotations. The substitution of u into Equation 3.7 yields the wave equation for fluid layers in terms of the scalar potential (p: V > = ^ ? 3.8 K dt2 The solution to Equation 3.8 is a scalar potential for a single compressional wave traveling through a fluid layer. For the situation indicated in Figure 3.1, two waves exist in the layer - one incident and the other reflected from the boundary B. Thus, the two quantities chosen to represent the sound field in the fluid layer are the normal components of velocity and stress at each boundary. As in Equation 3.2, a fluid-layer transfer matrix40 relates the velocity and stress at each boundary: 26 f U(/v 2 ) j 3.9 where T f = cos(£ 3 • L) k, , y'-sin(*3 -L) \_co-p — j -sin(*3 -Z) *3 COS(^3 • Z) 3.10 Elastic-solid layers Both normal and shear stresses occur in elastic solids, resulting in displacements according to Equation 3.3. Equation 3.3 is substituted into Equation 3.7, and the potentials cp and \\i satisfy the equation of motion for elastic solids if: P dV l + 2ju dt2 3.11 V 2\ |/ = p d2v ju dt2 3.12 Equation 3.11 is the wave equation for compressional waves; Equation 3.12 is the wave equation for a shear wave that propagates through elastic-solids. Therefore, according to Figure 3.1, four waves exist in the elastic-solid layer. There are two incident waves - compressional and shear - with corresponding reflected waves that propagate inside the elastic-solid layer. The four quantities at the layer boundaries are the normal and tangential velocity components, V3 and vi, respectively, and the normal and shear stress components, 033 and on, respectively. The system of equations for the four unknowns are represented by a solid transfer-matrix T s: 3.13 'v,(W2)~ v3(/V,) • - T • -v3(#2) , 0-33 W ) 0-33 (^2) c-13(iV,). ^ 1 3(/v 2). 27 The solid-layer transfer-matrix, T S ] was developed by Folds and Loggins22, but for different time dependence (e"jcot) and a different propagation direction (in the -x 3 direction). The solid transfer matrix has been re-worked, consistent with Allard's development for the porous transfer matrix (time dependence e^ 1 and incident wave traveling in the positive X3 direction), the details of which are given in Appendix A. 3.1.3 Elastic-porous layers Sound propagation through layers of porous-elastic material is described by Biot theory, which is essentially a modification of the stress-strain and energy relations for deformation of an elastic-solid. Porous materials consist of a solid component (frame) and a fluid component. In room-acoustical applications, air-saturated porous materials only are considered; hence, the fluid component is air at standard temperature and pressure. The fundamentals of Biot theory for wave propagation through air-saturated porous media have been reproduced here from Allard34. The stress-strain relations according to Biot are: cr*y= {(P - 2JV>'> Q9f }StJ.+ 2Ne'y cr{ =-<f>p = QOs +R0f In the context of Biot theory, the particle displacements are written for both the fluid and frame: us = Vcps } > displacements for compressional waves 3.15 uf =V(pf\ u ! =Vxi | / ' , [•displacements for shear waves. 3.16 = Vx \|/^ In the above equations, the superscript "s" denotes the solid or frame, and "f identifies the fluid (air). As in the case of elastic solids, the particle displacements and stresses are substituted into the equation of motion of Equation 3.7. The process, while similar to that of elastic-solids, contains many extra elasticity and inertial-coupling coefficients introduced by Biot 28 to account for fluid-frame coupling. The result is three equations for the unknowns <ps, (p\ and v|/s. There are two solutions (for the compressional waves) to the following: -ft»2M"'pc9 = V V 3.17 where M = <P = P Q Q R . A i PM Pn P22. where pn, pn, and p22 are inertial coupling terms and one solution (shear wave) to the following: V 11/ + — /V ^2 A PuPn ~Pn P22 \|/ = 0 3.18 J Equation 3.17 consists of two equations represented as an eigenvalue-vector equation for unknowns 69s and qf, that describes the propagation of compressional waves through the fluid and frame. The two eigenvalues are the squared complex wave numbers, 5i2 and 02 3 4 . The eigenvectors, 91 and (p2, are the vectors of solid and fluid displacement potentials for each wave. Equation 3.18 describes the single shear wave in terms of one vector displacement only for the solid phase - i.e. i|/ s = v|/. As in the case of fluid layers, the fluid in a porous material is assumed stationary and, therefore, shear waves do not propagate in the fluid phase. Depending on the degree of coupling, these waves may travel through both media simultaneously or independently. Strong coupling between the frame and fluid result in a fast and a slow compressional wave. These waves travel through both frame and fluid simultaneously; the frame and fluid displacements are identical for the fast wave and nearly opposite for the slow wave. Furthermore, the damping due to viscosity is much stronger for the slow wave. For ordinary air-saturated porous materials, the coupling is not very strong between 29 the air and frame because the frame is much denser. The two compressional waves exist as frame-borne and air-borne waves. For almost no coupling, the frame-borne wave propagates in the frame only, and the air-borne wave propagates in the air-filled pores only. Partial de-coupling is common, as the frame-borne wave induces vibrations in the air and, therefore, the frame-borne wave may propagate in both the frame and air. The air-borne wave has little effect on the frame, and propagates mostly in the air-filled pores. The shear wave is strictly a frame-borne wave. The acoustical properties at the boundaries of a porous layer are determined by six waves that propagate through the layer. These are the two compressional waves and the one shear wave that are transmitted through the layer, and the same three waves that are reflected from boundary B inside the layer, as indicated in Figure 3.1. The six independent quantities chosen to represent the sound field at the boundaries are the fluid and elastic-solid velocities and stresses - i.e. vis, V 3 S , v3f, 0 3 3 s , <Ti3S and cr33f. The transfer-matrix relation for an elastic-porous layer is: 'v{(N2)~ vl(N2) < v 3 /(/V 1) ' = T D v{(N2) > P 3^3(^ 2) ' o-;3(N2) 3^3 W ) . <(N2\ 3.19 The elements of the transfer-matrix T p are given in Appendix A, from the work of Allard29, who developed the transfer matrix for a porous layer for time dependence e*™1 and incident wave propagation in the positive X3 direction. It is important to note that the transfer matrix for the elastic-solid layer is a subset of the porous-layer transfer matrix, for which those rows and columns that contribute to the calculation of the fluid velocity and stress components are eliminated. Furthermore, some of the coefficients are changed to include only the elastic constants for solids (Young's Modulus and Poisson's ratio) rather than the Biot elasticity and inertial-coupling terms. 30 3.1.4 Multiple layers of fluid, solid, or elastic-porous For an individual layer, the transfer-matrix approach is sufficient to calculate the reflection coefficient, assuming an extended-reaction surface. However, the surfaces of real rooms rarely consist of a single layer or materialtype. Therefore, it was necessary to model the surfaces as a combination of layers that represent a typical construction, such as drywall, insulation (sound absorbing or porous material), plywood and air-gaps. The transfer-matrix method for single layers has been extended to include multiple layers of fluid, elastic-solid, or porous-elastic materials . This model involved the transfer matrices for the individual layer types described in Sections 3.1.1, 3.1.2, and 3.1.3, as well interface matrices that satisfy the boundary conditions at the interface of two adjacent layers, or layer and half-space (the space outside the layers). The situation is illustrated in Figure 3.2, There are 1 ... / number of layers, bounded by air on the incident side, and either air or a solid beyond layer /. If there is a solid after the last layer, the composite layer structure has a rigid backing. For any two adjacent boundaries, at points / and z+1 a set of homogeneous equations describes all boundary conditions for that interface. At the boundary between two adjacent layers, there must be continuity of displacement and force balance. Therefore, the set of homogeneous equations is expressed in terms of the velocity and stress components at the boundaries, and is related by two interface matrices, I and J. Nl N2 N3 N4 N5 N21-2 N21-1 N21 kil ^ ^ k i / >^krl ki2V kr/ kr2 Layer 1 Layer 2 B x3 Layer/ t xl Figure 3.2: Wave propagation through multiple layers of material. Incident wave is denoted by "i", reflected wave denoted with "r", and transmitted wave indicated by "t". 31 With reference to Figure 3.2, the interface matrices define the set of homogeneous equations for the multi-layered structure: • V{Nt)+ 3mn • V(NM ) = 0 within the composite layer structure I f n • V(NA) + 3 f n • V{NX) = 0 boundary between air and first layer ' ^ 0^ 2/) + Jm,n " V(NB) = 0 boundary between last layer and fluid or solid, 3.20 where m, n e {f,s,p} and i = 1,2,3,...21. With V(NM) = Tn • V(Ni+2), Equation 3.20 can be expressed as follows: V{Ni+2)=0 i = A,2,A,6,...2l,B 3.21 The elements of the interface matrices I and J for all combinations of adjacent layer types are given in Appendix A. These have been reproduced from the works of Brouard et. al.35, with the exception of the solid-porous interface matrix that was re-worked because of a mistake in the literature. The development of this interface matrix is shown below; the others were derived in a similar manner. The development of the solid-porous interface matrix begins with a statement of the boundary conditions for a solid-porous interface. With reference to Figure 3.2, layer 1 is a solid and layer 2 is a porous material. The boundary conditions at points N2 and N3 can then be summarized as: V;{N2)=V;(N3)=V({N3) 3.22a v;(/v2)=v;(iv 3) 3.22b (iV 2)=Cr; 3(7V 3) 3.22c 32 as33 (N2) = a{3 (/V3) + a;3 (N3) 3.22d The boundary conditions stated in Equations 3.22a - 3.22d are written as a set of homogeneous equations in matrix form: KP-V(N2) + JS,P-V{N3) = 0 V(N2) = {v;(N2) v°(N2) CT;3(N2) CT;3(/Y2)}rand 3.23 V(N3) = {vt(N3) v;(N3) v{(N3) CT;3(N3) a^3(N3) cr{3(N3)f The interface matrices I S j P and J S > P are, therefore: "1 0 0 0" "1 0 0 0 0 0" 0 1 0 0 0 1 0 0 0 0 0 1 0 0. and J V P = 0 0 1 0 0 0 3.24 0 0 1 0 0 0 0 1 0 1 0 0 0 1 0 0 0 0 1 0 For the case of layer 1 as a porous material and layer 2 as a solid, the interface matrices are interchanged; i.e. I P ; S = J S j P and J P ] S = IS ] P. The acoustical properties of a composite, layered structure with / number of layers are defined throughout the structure from points Nj to N21 by assembling all instances of Equation 3.21 at each boundary. The assembly of all boundary conditions yields a set of homogeneous equations, expressed in matrix form: 3.25 0 0 0 Dm,n=lm,n S n d Em,n = Jm,„ " T/ • D V = 0 12 D 0 0 "12 0 D 23 0 0 '23 D m-\,n-\ 0 0 0 ~"m-\,n-\ 33 For a fluid (air) backing: VD={p(A) v{(A) V{N2) V(NA) ••• V(N2l) p(B) v{ (B)f and for a rigid (solid) backing: VD={p(A) v{(A) V(N2) - V(N2l) v;(B) vf(5) <r'33'B) <r°33(B)}T where i = 1,2,..., / layer number, for / number of layers m = i n = i + l. At this point, the sound field is only defined inside the composite layer, from points Nj to N21. The sound field must also be known just outside the first boundary - at the point A in Figure 3.2 - requiring manipulation of Equation 3.25. To calculate the sound field at A, the boundary conditions at the A-Nj and N21-B interfaces were incorporated into the matrix D in Equation 3.25. The first step was to state the boundary conditions for the last boundary, the interface between points N21 and B. Two situations exist for this boundary, one for which the outside medium is a fluid (air); the other is the case is of the composite layered structure with a rigid backing (rigid solid). For a fluid-backed structure, the normal surface impedance at B is the perpendicular component of the characteristic impedance of the fluid; the last interface matrix in Equation 3.25 is Jmj (where m is either an s for solid, or p for porous material). The case of a rigid backing implies no displacements (and hence the velocities are zero) at the point B, and the last interface matrix is JmiS (where m may be either/for fluid, or p for a porous material). The two cases for the boundary conditions at the N2i-B interface are stated as follows: P(B) Z, c for air backing 3.26a vi (B) COS0 or 34 V3 (B) = v* (5) = 0 for a rigid backing. 3.26b Row(s) are added to the bottom of the D matrix in Equation 3.25 to include the boundary conditions from Equation 3.26a or 3.26b. The new matrix, D', becomes: or D' = D' = D 0 1 -ZB D for air backing 0 ••• 0 1 0 0 0 0 ••• 0 0 1 0 0 for rigid backing. 3.27a 3.27b The boundary conditions for the first interface, between points A and TV/, are similar to Equation 3.26a for an air backing, except that ZA is the unknown that is finally calculated. The boundary conditions at the A-Nj interface are: v{(A) 3.28 The boundary conditions in Equation 3.28 are implemented by adding a row to the top of D' (air or rigid backed), yielding the matrix D" that now completely describes the sound field from points A to B: D" = 1 0 D' 0 3.29 Next, the set of homogeneous equations that completely define the sound field throughout the composite-layer structure from point A to B is represented as D" • VD = 0, with VD defined in Equation 3.25 for / number of layers. From this matrix relation, the normal surface impedance was calculated and, hence, the reflection coefficient was found and used as the input 35 to the beam-tracing algorithm to model the reflection of sound waves from extended-reaction surfaces. 3.1.5 Surface impedance and reflection coefficient The solution for the normal surface impedance at A , ZA, was obtained from D" • VD = 0. As a set of homogeneous equations, a unique solution exists only if the determinant of D" is zero. Therefore, ZA was determined from D" in Equation 3.29 as follows: |r>"| = |z>;|+^ -|£> 2 | = o therefore, 3.30 where \Dj "\ is the determinant of D' with the first column removed and |£>2 '| is the determinant of D' with the second column removed. With knowledge of the normal surface impedance ZA, the reflection coefficient at point A was calculated, for angle of incidence 0, as follows: zA-zj«m Z,+Z,/cos(0) This is the complex pressure reflection coefficient calculated for a set of given material properties for all layers in the surface; it yields the pressure amplitude and phase change for the reflected plane wave when multiplied by the complex pressure for the plane wave incident on this surface. This reflection coefficient was applicable to the beam-tracing algorithm to determine the reflected sound wave from a surface, including both magnitude and phase changes. Furthermore, the reflection coefficient calculated in this way correctly modeled the extended-reaction nature of the surfaces incorporated into this room prediction-model. 36 The transfer-matrix method is applicable to extended-reaction surfaces as the calculated surface impedance is inherently angle-dependent, due to the wave-number components in each layer. The wave-number components are in the normal and tangential directions to the surface, and ki, respectively. The magnitude of these components depends on the angle of incidence; the fact that the kj component is non-zero at angles other than normal incidence implies that the surface is of extended reaction. This is in contrast to the assumption of local reaction, for which the normal surface impedance is independent of the angle of incidence, and the ki component of the wave-number vector is always zero. Only the ks component is non-zero and constant for a given frequency. Theoretically, local reaction is only correct when a wave propagates through a layer in a direction along the surface normal. The local-reaction surface impedance is equivalent to the surface impedance at normal incidence36. Therefore, the transfer-matrix model can be used to model locally-reacting materials by predicting the surface impedance at normal incidence. The local-reaction surface impedance, as with an extended-reaction surface impedance, is frequency dependent. The reflection coefficient is still angle-dependent and is obtained from Equation 3.31, but the surface is of local reaction. 3.2 Beam tracing The beam-tracing algorithm was chosen to model sound propagation in an enclosed space occupied by air at standard temperature and pressure, referred to here as a typical room. Beam tracing was identified as a computationally-efficient method to model sound propagation in a room, although it is typically used as an energy-based method valid only at higher frequencies. Extended to include phase, and combined with the transfer-matrix method outlined in Section 3.1, beam tracing was the basis of the new room-prediction model applicable to rooms with extended-reaction surfaces, over a broad frequency range. As stated from the research objectives, the triangular-beam tracing method used was the one described by Lewers15. Only the concepts and issues associated with Lewers' beam-tracing method were discussed in his paper and, therefore, the algorithm and geometric calculations involved needed to be developed. The algorithm was developed by first identifying the individual processes involved with beam tracing in rooms. The main components to the beam-tracing algorithm consist of: 37 a) Surface geometry-definition b) Source generation c) Trajectory loop for a beam d) Receiver-point detection, and e) Calculation of the steady-state sound-pressure level. Surface definition, whereby the room geometry is defined, is discussed in Section 3.2.1. The source, discussed in Section 3.2.2, is an omni-directional point source radiating spherical waves, with the same sound power at all frequencies. The fundamental component of the beam-tracing algorithm is the trajectory loop, illustrated by the block diagram in Figure 3.3. The trajectory loop, executed for each beam, combines the process of determining the nearest plane, calculating the new image source, and calculating the new direction vectors and angle of incidence of the reflected beam. Each beam originates at the defined source, and the trajectory loop is executed for each reflection number, until the maximum reflection order is reached. Upon completion of the trajectory loop, the beam number is incremented and the next beam is traced, with direction vectors reset to the original source. The components to the trajectory loop are discussed in Sections 3.2.3 - 3.2.5. After each execution of the trajectory loop, a test must be made to determine if the receiver point is contained within the path of the beam from the previous source to the new image source. The receiver-point test is described in Section 3.2.6. If the receiver point is contained within the beam's path, the pressure contribution of the current beam is calculated at the receiver point, described in Section 3.2.7. Upon completion of tracing all the beams, the complex pressure contributions of each beam are added and the results converted to decibels, to yield the steady-state sound-pressure level at the receiver point. All calculations are performed in the frequency domain. This beam-tracing algorithm used beams with a triangular cross-section, as illustrated in Figure 3.4. Beams with triangular cross-sections were used because a point source radiating spherical waves could be modeled with no overlapping beams and an even energy distribution across the beam face. Furthermore, the triangular cross-section of a beam represented a plane wavefront, which is assumed for an incident sound wave on a surface in the application of the transfer-matrix method to calculate the complex reflection coefficient of that surface. A single 38 Beam orig inates at source Si, with direction vectors bi for vertices, a for center Find the reflecting plane i Determine new source, Si+i, that is an image of the previous source, mirrored about the reflecting plane. Determine the angle of incidence and new direction vectors bi+i and ci+i for the reflected beam Reset source: Si = S i+l Reset direction vectors: bi = bi+i and ci = ci+i ^ Yes Increment the beam number and reset source and direction vectors to their starting values Figure 3.3: Trajectory-loop algorithm, beam is identified by three bounding rays, (vectors bi, b2 and b3), center vector c, three bounding planes (Bi, B2 and B 3 ) , and cross-section plane X. The vertex or origin of the beam is the single point, S. A beam's vertex is the real source in the room if it corresponds to the direct source-39 receiver path. Otherwise, for a beam corresponding to a source-receiver path due to reflections, the vertex of that beam is the image source that was last calculated. 3.2.1 Defining surface geometry The basic room is constructed with rectangular surfaces, each defined by a set of four points in Cartesian coordinates. For surfaces that are non-rectangular in shape, more coordinates are needed to define the boundaries as necessary; for example, six points are needed to define a hexagon. In order to perform the ray/beam-tracing calculations, the normal vectors and plane coefficients of each surface must be calculated. Several rules for defining the surface geometry are required to ensure consistency in the ray/beam calculations: • The room surfaces, indicated by Th with surface number are defined in the positive quadrant of the Cartesian coordinate axes x\,X2, and x 3. • For TL, the unit normal vector n,- points outward (away from the inside of the room, or from the source). • The surfaces Tii are defined as a plane equation of the form AjX\ + B&i + Qx 3 = D„ where Ah Bt, C,, and £>,• are the coefficients of plane i. The coefficients A,, Bh and C, are also the components of the normal vector N,- of plane i. The coefficient £>,• is the perpendicular distance from plane i to a parallel plane that contains the origin. Figure 3.4: Geometry of a triangular-faced beam. 40 Regardless of the irregularity of the surface, three points are required to calculate the surface geometry necessary to define a room boundary in the beam-tracing algorithm. The required geometric calculations for surface i are the surface normal N„ and plane coefficients Ah Bh d, and £),•• The plane configuration is shown in Figure 3.5, relative to the coordinate axes of the x/, X2, and X3 system and the source S. The normal vector to the plane and the plane coefficients are calculated as follows: N , = ( P 2 - P 1 ) x ( P 3 - P 1 ) n N ' n; = 1—r N Ai=Nn,Bi = iV , 2 ,C , . = iV / 3 ,A=N , - P 1 Note that the plane is defined in the positive xj, x2, and X3 axes, and the normal vector points away from the source (interior of the room). These rules, as described above, facilitate the computations to produce the correct distances and angles, and to test for surfaces, as described in the following sections. 3.2.2 Modeling of an omni-directional source For the purpose of this room-prediction model, a single source was used that was omni-directional at all frequencies. Furthermore, an ideal spherical source was modeled that approximates most sources found in room-acoustical applications, such as a loudspeaker or a human speaker (though they are not omni-directional at high frequency). The base icosahedron Figure 3.5: Surface definition. The unit normal N points outward, away from the source S, such that N»S<0. 41 (20-sided geometric object) was used to model a spherical source with center S, with each face representing the triangular plane X of a beam. The geometry of the icosahedron was developed by Campbell4 1 who used the golden ratio, r, to calculate the vertices: 3.33 With the 12 vertices located at coordinates given by: (0,±r,±l),(±l,0,±r),(±r,±l,0) 3.34 The icosahedron geometry was modified to model a point source with a spherical wave-front that consisted of more than the 20 faces (beam fronts) of the base icosahedron. To do this, the base face of the icosahedron was defined in the coordinate system %\, £2 and £3, and each edge was divided into equal lengths, the number of which is referred to as the frequency,/(not to be confused with the frequency of sound waves); this results in a polygon with 20xf2 faces or beam fronts. In this coordinate system, a single face is represented by vertices at (0, 0,f), if, 0, 0) and (0 , / 0). The method developed by Campbell was further refined to calculate the vertices of just one face - the base face - and perform coordinate transformations to get the other faces. An example of the icosahedron with frequency / = 3 is shown in Figure 3.6, in the coordinate system of £1, £2, and £3. The icosahedron face shown in Figure 3.6 was transformed to its real representation, described in the x\, x%, and X3 coordinate system, with the following relations: Equations 3.35 defined the base face of the icosahedron in the x\, xi, and X3 coordinate system. The base-face coordinates were further manipulated so that all vertices were points on the surface of a unit sphere centered at the origin, and one edge of the base face was lying along the x\ axis. The remaining icosahedron faces were calculated by a coordinate transformation x, = £ sin(72°) x2 = £ 2 +£,cos(72°) x 3 = / / 2 + & / r 3.35 42 (3, 0, 0) (2, 1, 0) (1,2,0) (0, 3, 0) Figure 3.6: Breakdown of an icosahedron face; the icosahedron frequency is 3. using the Euler rotation angles y/, 6 and <f>. The base face n was rotated by the appropriate angles to find the new face m in the x\, X2, and X3 coordinate system: = E- 3.36 where cosi/cos^ - sin I/cos#sin^ sin^cos^ + cos^cos^sin^ sin#sin^ E = I-cos^sin^ - sin 1//cos cos ^  -sin^sin^ + cosi/ycos#cos^ sin#cos^ sin ^  sin # -cosy sin # cos (9 A further refinement was made by eliminating all duplicate vertices and forming a matrix that identified each face (along a row) with a set of three coordinate numbers (columns) corresponding to the elements of the vertex matrix that contained the vertices for that face. The coordinates were also scaled so that all vertices lie on a unit sphere. Together, the vertex matrix Vs and face matrix Fs for the source identified the starting coordinates for the vertices bi, b2 and b3 for all beams. A source constructed with 320 faces (Vs [162x3] with 162 unique vertices, and Fs [320x3]) is shown in Figure 3.7. 43 Figure 3.7: Source with 320 beams (Icosahedron of frequency 4). The center ray, c, for each beam originating at the source was calculated as the unit vector that was the unit normal to the plane that contained the vertices bi , b 2 and b 3. This plane was parallel to a plane tangent to the unit sphere, at the point that was intersected by the center ray for that beam. The starting coordinates of the center ray for a given beam were calculated: N c = ( b 2 - b , ) x ( b 3 - b , ) 3.37 C = 1— 7^ K I The following check is made to ensure that the centre ray points outward: fc i f c - b , > 0 [-c if c b , < 0 The geometry of the source was completely described by the vertices bi , b 2 and b3j and the centre ray c for all beams. With this geometric information, beams were traceable, but further acoustical information was required. The beams initially modeled a point source radiating spherical waves and, therefore, the initial pressure amplitude was required. A source is characterized by its radiated sound power; this is independent of the environment of its location. 44 The initial sound pressure was calculated from the initial sound power of the source, in a free field. A spherical vsave,p(r, f) is radiated from the source and expressed as: p{r,t) = --e^ 3.38a r The amplitude constant A, is determined from the sound intensity of the source. At distances sufficiently larger than the wavelength, such that kr » 1, the sound intensity I (in units ofW/m2), is: , p2 A2 I = = — 3.38b Po •c r • Po •c The sound power, W, is found by multiplying the area of a sphere by the sound intensity: -ITT A ? T ^ * 7t ' A W = 4-7r-r2-I = 3.38c Po • c Therefore, A = 3 V 2-n The initial phase angle was set to zero, completing the information required to begin the beam tracing. The beam tracing essentially models a spherical wave propagating from a point source. The sound pressure of the spherical wave decays as \IL, at a distance L from the source, with phase angle -kL, where k is the wave number. 3.2.3 Finding a reflecting plane One of the fundamental components of the beam-tracing algorithm, and of the implementation of extended-reaction surfaces, was to determine the surface of reflection for a given beam with direction vectors and source point. Determining the surface of reflection involved two steps: first, the total number of room surfaces was reduced to include only the 45 subset of permissible surfaces (i.e. those surfaces that are in the path of the centre ray); then, the distance from S to the point of intersection with each of these surfaces, in the direction of c (the centre vector), was calculated. The surface that corresponded to the minimum distance s, was the surface that was intersected by the beam. For a beam b, with center direction vector c* and start point S, the condition for permissible surfaces was cb -N,. > 0, for surface i with normal vector N,. For each permissible surface, the distance from S to each point of intersection, in the direction of c^ , was calculated, with the aid of Figure 3.8 as follows: d s = • cos#2 3.39a d = D + \S\-cos6x N,.-S n,-S costf, = - , ', , , = - j ^ ^ 3.39b N ' I S I l n <H S l :.d = D + \S\--. • =D-n. S K l - N N, • c, n, • cA COS02 =—'•—b— = ^ — = 11,. -ch N,. • cb n,. • c„ 3.39c Therefore, , = D - - - 8 3.39d 3.2.4 Finding an image source For a beam originating at point S, and intersecting surface IT, the image source Sj+i is found as follows: S i + 1 =S, +2-sP -n n 3.40 where sp is the perpendicular distance from the source S; to the image-source mirror plane. 46 Figure 3.8: Distance from S to intersection point on a plane. The geometry of a line (center ray of a beam) intersecting a plane (the room surface) in the X1X2 plane is shown. This is shown in two dimensions for clarity; however the calculations are applied to three dimensions. In three dimensions, the vectors S, c*, and n/ are all co-planar. The distance D is the perpendicular distance from the plane to a parallel plane that contains the origin. D is the plane coefficient from Axi+Bx2+Cx3=D. In Equation 3.40, the perpendicular distance sp is obtained by substituting the unit-normal vector n, for the center unit vector c^ , in Equation 3.39d. The image source becomes the new start point for the reflected beam; the direct path from the image source to the beam front is equivalent to the reflected path from the previous source to the beam front. 3.2.5 Angle of incidence and direction vectors for the reflected beam The angle of incidence of a beam that intersects a surface was required for use as input into the transfer-matrix model to calculate the complex reflection coefficient. The angle of incidence was calculated as the angle between the center ray and the surface normal, as shown in Figure 3.9. From Figure 3.9: u n cost/ = = u n u • n Therefore, 0 = cos _ I ( i i - i i ) 3.41 47 The angle of incidence and the frequency, together with the material properties of each layer in the surface, provide the input parameters required to calculate the reflection coefficient by the transfer-matrix method. The complex reflection coefficient was multiplied by the pressure of the incident wave at the point of intersection, and the pressure magnitude and phase of the reflected wave were determined. In addition to the angle of incidence, the direction vectors were calculated for the reflected beam. The reflected beam was the beam originating from the image source calculated in Equation 3.40, with new direction vectors b i , i , b i ,2 , b i ,3 , and Ci that represent the first-order reflection vectors of the incident direction vectors b n , i , bo,2, bo,3, and cn. Therefore, for reflection order r, the reflected direction vectors are obtained from the incident direction vectors as the set { b r - / , i — » b , - , i , b r . / , 2 — » b r , 2 , b r./,3—>b r,3, cr.i—>cr}, and are represented in general as u—»v where u is the incident ray and v is the reflected ray. The new (reflected) unit direction vectors represented individually as the vector v, are calculated with the aid of Figure 3.9, as follows: v = u + |u| • cos 0 • (-n) + |v| • cos 6 • (-n) Since u and v are unit direction vectors, lui = Ivl = 1, and cos# = " ° = u • v and, |u| • |v n Figure 3.9: Direction vector of a reflected ray. 48 therefore, v = u-2 (u n)n 3.42 The new reflected beam with the new direction vectors and image source are used to define a new incident beam, and the reflection order is increased by one. The processes described in Sections 3.2.3 - 3.2.5 combine to form the trajectory loop for one beam; they are repeated a number of times equal to the maximum reflection order that is pre-defined. The trajectory loop may also be repeated less than the maximum reflection order if another termination criterion, such as beam frontal area or energy level, reaches a critical value, also pre-defined. Within the trajectory loop, there is a test to detect the receiver point, as described in the next section. At the receiver point, for a positive detection for a given beam, the complex acoustical pressure is calculated and recorded for the current reflection order, as described in Section 3.2.7. 3.2.6 Receiver detection The motivation in developing any room-prediction model is to determine the sound field at a specified receiver position, R, in the room. The possibility that the receiver point lies somewhere inside the volume of a beam that propagates from the start point S to the intersecting surface n must be tested for; that is, a test is done on the propagating beam to determine if it is received. A test algorithm for the receiver point detection was developed by Campbell41. Three tests were performed - one for each boundary plane of a beam, to determine on which side the point R was located. Each test evaluated to true if the point R was located on the side of the plane such that it was inside the volume that was traced by the beam. For all three tests evaluated as true, the receiver point was contained within the beam. For the successful detection, the complex pressure amplitude was recorded for the path length L of the beam - from the point S to the point of intersection by the center vector c and the plane X that contained the receiver point R. The geometry of the beam/receiver detection configuration is illustrated in Figure 3.10. The beam/boundary plane receiver-point test is performed by evaluating a set of three equations within each instance of the trajectory loop. The receiver point R is contained within a beam if all three conditions are true: 49 B, Figure 3.10: A receiver point that is detected by a beam. The receiver point R lies in the plane X. (R-S) .(R-S) NB3-(R-S)J N N Bl B2 >0 3.43 N B 1 = N B 2 = N B 1 f b x x b 2 i f b , x b 2 > 0 [ -b j x b 2 otherwise [ b 2 x b 3 i f b 2 x b 3 > 0 l - b 2 x b 3 otherwise f b 2 x b 3 i f b 2 x b 3 > 0 [ - b 2 x b 3 otherwise. After a positive receiver-detection result, the complex pressure amplitude for that beam was calculated at R and the trajectory loop continued, until a termination criterion was satisfied. Therefore, it was possible for the same beam to hit the receiver point more than once (but at different reflection orders). 3.2.7 Calculating the steady-state sound-pressure level Following a positive receiver-point test, the complex pressure amplitude for that beam was calculated. The pressure is calculated for the beam starting at S, to the point of intersection 50 between the center ray c and the beam-front plane X that contained the receiver point R. The beam's contribution to the acoustical pressure at the receiver point R is calculated as follows: ('A\ pr(L,co) = — \ L J •Reffe~JkL 3.44 The result is complex pressure (magnitude and phase) at the point R for a single frequency and for a single beam between reflections r and r + 1. For the same beam, all pressures recorded at the receiver point are added to yield the contribution to the steady-state pressure of that beam. The total steady-state acoustical pressure is calculated by summing the total pressures of all beams detected at R. The steady-state sound-pressure level is calculated to express this result in decibels (dB). 3.3 Implementation Both the transfer-matrix and beam-tracing models were implemented via a computer program in the Matlab environment, on a Windows-based P C . The programs were written using Matlab language and functions as a collection of m-files. The transfer-matrix and beam-tracing programs were written independently so that they could be used separately. The transfer-matrix program is a function that outputs acoustical surface properties for a multi-layered surface, given the material properties of the layers. The beam-tracing program is a script file that performs the trajectory loop for each beam, and computes the sum of the pressures recorded at the receiver point. The transfer-matrix function is called within the beam-tracing program to calculate the reflection coefficient for a surface, at the calculated angle of incidence, and over the specified frequency range. The code for all of the m-files is given in Appendix C . 3.3.1 The transfer-matrix program The transfer-matrix function implements the theory presented in Section 3.1, to model a surface as a series of unbounded layers with finite thickness. The outputs are the acoustical surface properties - surface impedance, reflection coefficient, absorption coefficient, transmission coefficient (if applicable) and transmission loss (if applicable) - for the modeled 51 surface. Each predicted surface property is a matrix, in which an element with indices [i, j] corresponds to the surface property at the / t h frequency and the / h angle. The inputs to this function are the frequency, angles of incidence, material properties of each layer, and a rigid-backing indicator. The material properties for each layer are input as a structure-variable type, with dimensions equal to the number of layers, and with separate fields that contain the solid and fluid properties of each layer. The solid and fluid properties of each layer type - fluid, elastic-solid, or porous-elastic - are given in Table 3.1. The rigid-backing indicator is set to a value of '1' to indicate there is a rigid backing, or '0' to indicate no rigid backing. A surface that consists of one or more layers backed by a thick concrete slab can be considered to have a rigid backing. The algorithm for the transfer-matrix function is shown in Figure 3.11. The heart of this function is a double-loop that computes the surface impedance at all angles, at each frequency. The surface impedance is calculated from the assembled set of homogeneous equations that define the boundary conditions at each interface, using Equation 3.30. The reflection coefficient is calculated from the surface impedance, for each frequency and angle, using Equation 3.31. The Table 3.1: Material parameters applicable to fluid, solid, and porous layer types. A layer with properties marked with a ">^ " indicate a required property for that layer. Those properties marked with a "*" are not required for that layer. Material Property Fluid Layer Solid Layer Porous Layer Thickness (mm) • • Density (kg/m3) • • • (fluid and frame) Sound speed (m/s) X X Std air pressure (kPa) X X •(fluid) Dynamic viscosity (Ns/m ) X X •(fluid) Specific heat ratio X X •(fluid) Prandtl number X X •(fluid) Young's Modulus (Gpa) X • X Shear Modulus (Mpa) X X • (frame) Poisson ratio X • • (frame) Flow resistivity (Ns/m4) X X • (frame) Porosity X X • (frame) Tortuosity X X • (frame) Viscous dimension (m) X X • (frame) Thermal dimension (m) X X • (frame) 52 • Import layer material properties • Initialize internal variables For each frequency, f For each angle of incidence, 6 For each interface in the surface, w i= l:(number of layers + 1) * -Determine the left and right layer types -Find the appropriate interface matrices -Construct the transfer matrix for the right layer type 1 -Construct the D matrix with the sub-matrices IV(n) and JTV(n+2), for the current layer and interface. Return matrix for each acoustical property, for all frequencies and angles Figure 3 . 1 1 : The transfer-matrix algorithm. assembled set of boundary conditions contained the interface and transfer matrices given in Appendix A, and the fluid-layer transfer matrix in Equation 3.10. It is important to note that the computation of the solid and porous transfer-matrices, involving matrix inversion, required a pseudo-inverse function, as these matrices were near singular. The transfer-matrix program is used to calculate the reflection coefficients for the extended-reaction surfaces, but also those required for local-reaction surfaces. This was necessary to compare the effects of extended- vs. local-reaction surfaces in the prediction of the steady-state sound-pressure level in rooms. Predicting the surface impedance at normal incidence, and calculating the reflection coefficient using Equation 3.31 for all angles, easily accomplishes this. 3.3.2 The beam-tracing program Programming of the beam-tracing algorithm involved implementing three fundamental components - source-generation, trajectory-loop and receiver-detection. The source was constructed as an icosahedron with a specified number of faces or beams, as shown in Section 3.2.2. The trajectory loop was implemented using the theory developed in Sections 3.2.3 - 3.2.5; a function for the receiver detection was programmed from the theory presented in Section 3.2.6. The beam-tracing algorithm, illustrated in Figure 3.12, begins with source generation, loading data describing room geometry and surface properties, and initialization of sound pressures. The sound pressures that represent the total pressure contribution of each beam at the receiver point are initialized to zero. Beams are traced via a beam loop that consists of a receiver-detection test and trajectory loop, and are repeated for the specified reflection order. A detected receiver point results in the calculation of the pressure contribution of the current beam for the current reflection number, according to Equation 3.44, and its addition to the total pressure contribution by that beam. Upon completion of the beam loop, the pressure contributions by each beam are added to yield the frequency response of the steady-state sound-pressure level. The steady-state pressure is determined from the effective reflection coefficient - the product of all reflection coefficients calculated at each reflection. The calculation of the reflection coefficient for a given surface depends on how that surface was defined. For multi-layered surfaces, whether of extended- or local-reaction, the reflection coefficient is calculated via the transfer-matrix function. Surfaces may also be defined by an impedance that varies with 54 /- Import input data from file - Initialize pressure contributions of each beam to zero - Initialize other internal variables 5 fe For beam number b, w b = 1 to Number of beams * For reflection order r, r = 1 to Number of reflections Receiver \. Yes detected? Calculate pressure contribution of current beam No Execute trajectory loop < Yes * Calculate steady-state sound-pressure level from contributions of all beams Figure 3.12: T h e beam-tracing a lgor i thm. frequency only (local reaction), or with an angularly-constant reflection coefficient. Specifying the reflection coefficient explicitly is the most direct way to define the acoustical property of a surface, in the context of the beam-tracing model. This is useful in outdoor or free-field applications, in which surfaces representing open air are modeled with reflection coefficients of zero. The combined beam-tracing and transfer-matrix program was applicable to comparing steady-state sound-pressure level predictions in rooms with extended- vs. local-reaction surfaces - the research objective. Extended-reaction surfaces are always assumed in the transfer-matrix model; however, locally-reacting surfaces can be modeled as well. Locally-reacting surfaces were modeled with the transfer-matrix model by predicting the surface impedance at normal incidence. The surface impedance calculated at normal incidence was used to calculate the reflection coefficient at all angles. The reflection coefficients were input into the beam-tracing algorithm for each beam reflection, to calculate the effective reflection coefficient. The effective reflection coefficient described the cumulative effect of magnitude and phase changes in the pressure contribution of a beam, for all reflections, from source to receiver point. The steady-state sound-pressure level at the receiver point was computed by summing the contributions from all beams. The next step in applying the combined beam-tracing and transfer-matrix model to meet the research objectives was to validate the models using published or known theory. In the next chapter, several validation cases are predicted for various surfaces using the transfer-matrix model, and for several room configurations using the beam-tracing model. After validation, the combined beam-tracing and transfer-matrix modelis applied to realistic rooms, to meet the research objective, as described in Chapter 5. 56 C H A P T E R 4 M o d e l V a l i d a t i o n In the last chapter, theory was developed to implement a transfer-matrix model and a beam-tracing algorithm. The transfer-matrix model predicts the acoustical properties of extended-reaction surfaces. Locally-reacting surfaces can also be modeled by the transfer-matrix method by predicting the surface impedance at normal incidence. The surfaces are modeled as a series of unbounded layers of finite thickness and made of materials classified as one of fluid, elastic-solid or elastic-porous. The beam-tracing algorithm simulates the propagation of a spherical wave, with pressure magnitude and phase, emitted by a point source in a room. The predicted result is the steady-state sound-pressure level at a receiver point. The beam-tracing model is applicable to specular reflections only, and for a reflection order suitable to predict the steady-state sound-pressure level. The transfer-matrix model was incorporated into the beam-tracing model to calculate the reflection coefficient of multi-layered surfaces that were either extended- or locally-reacting. A combined model was, therefore, in place to satisfy the main research objective - to quantify the differences between extended- vs. locally-reacting surfaces on steady-state sound-pressure level in rooms. Before the new model could be applied, however, validation was necessary. Validation of both the transfer-matrix and beam-tracing models was necessary to apply the new room-prediction model to real rooms with confidence. Since there are no other known room-prediction models of this kind, the transfer-matrix and beam-tracing models were validated separately. The validation of the beam-tracing model was also necessary to determine the number of beams to use for a given room size. With confidence from separate validation tests 57 that both models were accurate, the combined model was assumed valid. A discussion of the validation process followed to satisfy the fourth research objective is given in this chapter. 4.1 Validation of the transfer-matrix model The validation of the transfer-matrix model involved predicting the reflection coefficient and surface impedance of surfaces that are commonly found in rooms. Surfaces were chosen to best represent real surfaces in rooms that are of acoustical interest, and for which measurement data or comparative models exist. For the chosen surfaces, the surface impedance was investigated for its angle dependence, to confirm the surface as extended reaction. As these surfaces were shown to be clearly of extended reaction, they were used in rooms studied in the application of the room-prediction model, in Chapter 5. The following surface test cases were chosen for validation, which was done as described: a) Single glass plate - transfer-matrix predictions were compared to plate-theory predictions. b) Double-drywall plate structure with air-gap - transfer-matrix predictions were compared to plate-theory predictions. c) Glass-fibre porous material with rigid-backing - transfer-matrix predictions were compared to measurements performed by Allard42, and to predictions by the Allard-Attenborough26 rigid-frame porous model. d) Carpeted floor, modeled as a double porous-layer with rigid-backing - transfer-matrix predictions were compared to measurements performed by Brouard et. al.43. The material properties of the plates and porous materials are given in Tables 4.1 and 4.2, respectively. The data for cases involving plates were easily selected from data available for a broad selection of elastic-solid materials. The choice of porous materials, however, was limited because of the limited use of Biot theory for sound-absorbing materials at the present time. However, Biot parameters were obtained for the fibre-glass material that was used in case c)42, and for the carpet used in case d)43. The desired acoustical property to compare predictions in all 58 cases was the reflection coefficient. This was possible for the two cases involving the single and double plates using plate theory, but not possible for the two porous-material cases. The transfer-matrix application to porous materials was validated by comparing the surface impedance predicted at normal incidence to experimental results. The reflection coefficient, however, is easily obtained from the surface impedance using Equation 3.31. For all cases, the surface impedance was predicted at several angles of incidence using the transfer-matrix model, in order to illustrate the extended-reaction nature of these surfaces. The same frequency range and angles of incidence were used in all validation tests. The frequency range used was 50 - 4000 Hz, the frequencies of interest in room acoustics, with 10 Hz increments. For each frequency in this range, the reflection coefficients and surface impedances for all cases were predicted at angles of incidence ranging from 0° - 90° using the transfer-matrix model. For the first two cases involving plates, the reflection coefficients were predicted with plate theory for these angles as well. Only the reflection coefficient and surface impedance predicted at angles of 0°, 45°, and 75°, at all frequencies, are shown in subsequent figures. Upon inspection of the results, it was determined that these three angles best displayed, as clearly as possible, the angle and frequency dependence of the predicted acoustical properties. 4.1.1 Single plate The application of the transfer-matrix model to a single plate was validated by a comparison of transfer-matrix and plate-theory predictions of the reflection coefficient, at various angles of incidence. Plate theory is commonly used to calculate the reflection coefficient and transmission loss of plates and, therefore, is appropriate for the validation of the transfer-matrix model in this case. Plate-theory prediction of plates as layers with finite thickness and infinite extent was conveniently expressed in a matrix form similar to the formulation of Brouard et. al.35 for a porous material bonded to a plate. The matrix formulation for plate-theory prediction was modified for the single-plate case and is described in Appendix B. The plate-theory predictions of the reflection coefficient are shown in Figure 4.1, to illustrate the acoustical behaviour of this surface. The transfer-matrix predictions of reflection coefficient are compared to the plate-theory predictions for angles of incidence of 0°, 45°, and 75°, in Figures 4.2 - 4.4, respectively. 59 Table 4.1: Material properties of the glass and drywall plates. Material Property Glass Drywall (Gypsum board) Thickness (mm) 6 12 Density (kg/m3) 2500 1200 Young's modulus (GPa) 60 7 Poisson ratio 0.33 0.25 Table 4.2: Material properties of the porous layers used in validation cases c) and d). Both layers have a rigid backing; the carpet consists of two porous layers. Properties for the air-saturated pores of all porous materials are as follows: Air density = 1.2 kg/m3, sound speed = 344 m/s, standard air pressure = 101.3 kPa, dynamic viscosity = 18.2xl0'6 Ns/m2, specific heat ratio = 1.4, and Prandtl number = 0.710. Material Property Fibre-glass Carpet Layer 1 Layer 1 Layer 2 Thickness (mm) 56 3.5 3.5 Frame density (kg/m3) 130 60 . 60 Frame shear modulus (MPa) (2200 +j)x 10"3 (10+>5)x 10"3 (10+y5)x 10"3 Poisson ratio 0 0 0 Flow resistivity (Ns/m4) 40000 20000 5000 Porosity 0.94 0.99 0.99 Tortuosity 1.06 1 1 Viscous dimension (m) 0.56 x 10"4 1.5 x 10"4 2.3 x lO - 4 Thermal dimension (m) 1.12 x 10"4 2.2 x 10"4 2.8 x 10"4 It is expected that the reflection coefficient of a single plate is mostly real and approaches a value of one. This is confirmed upon inspection of Figure 4.1, for most angles, except at 75° at which a sharp spike appears at approximately 2200 Hz for both the real and imaginary components. The real part decreases rapidly towards zero; the imaginary part increases rapidly from 0 - 0.5, with a sharp decrease down to -0.5 followed by an increase to zero, over a narrow frequency band. This indicates almost total sound transmission through the plate, with a 180° phase change, for the incident sound wave with an angle of incidence of 75° at 2200 Hz. The phenomenon is known as 'coincidence', a condition that occurs in plates for an incident wave at oblique incidence, at a single frequency. The frequency at which coincidence occurs, known as 60 the coincidence frequency, depends on the angle of incidence and the material properties of the plate. The coincidence frequency is the frequency at which the projected wavelength of the wave, obliquely incident on a plate, matches the bending waves in the plate due to the acoustical pressure loading. The reflection coefficient of the glass-plate was predicted using the transfer-matrix model and was compared to the plate-theory predictions. The results, shown in Figures 4.2 - 4.4, indicate very good agreement between the two models. The only noticeable difference was at 75°, where the transfer-matrix model predicted a coincidence frequency slightly higher than that predicted by plate theory. A possible reason for the discrepancy is in the assumptions made in plate theory. Plate theory assumes a neutral stress plane at the center plane of the plate, and that plate cross-sections remain plane and undistorted. On the other hand, in the transfer-matrix 500 1000 1500 2000 2500 Frequency (Hz) 3000 3500 4000 0:5 i o o c o -0:5 ; • •  i: i i . 1 j ^ J i i i 6=75° 6=0°, 6=45° 500 1000 1500 2000 2500 Frequency (Hz) 3000 3500 4000 Figure 4.1: Reflection coefficient of a single glass plate, 6 mm thick, predicted with plate theory. Real parts are shown in A, imaginary parts in B. For angles of 0°and 45°, the curves are almost superimposed. 61 formulation of an elastic-solid layer, plate cross-sections can be distorted and there is a difference in the tangential displacements as predicted by the two theoretical models. The transfer-matrix results would be expected to be more accurate as it models the surface displacements more accurately. Using the transfer-matrix model the surface impedance was predicted as a function of angle and frequency, to illustrate that the surface is of extended reaction. The surface-impedance curves for the three angles used in the validation cases are shown in Figure 4.5. The impedance has been normalized with respect to the characteristic impedance of air, 415 rayls. If the glass plate were locally reacting, the impedance curves in Figure 4.5 would be superimposed, indicating no variation with angle of incidence. This is, however, clearly not the case for the glass plate. A strong angle-dependence of the surface impedance is seen in Figure 4.5. The real part is almost constant with frequency, and increases with increasing angle. The imaginary part is dominant, especially at higher frequency, with angle-dependence at frequencies above 1000 Hz. At higher frequencies, the imaginary part decreases rapidly with increasing angle. These results confirm the extended-reaction nature of the single glass plate. 0.2 Figure 4.2: Prediction of the reflection coefficient of a 6 mm single glass plate, at normal incidence. Prediction with plate theory: —. Prediction with transfer-matrix model: — . Real part is shown in A, imaginary part in B. Note that the curves are superimposed. 62 1.02" Frequency (Hz) Figure 4.3: Prediction of the reflection coefficient of a 6 mm single glass plate, at 45°. Prediction with plate theory: —. Prediction with transfer-matrix model: — . Real part is shown in A, imaginary part in B. Note that the curves are almost superimposed. 1-2: 0> 5=o:e; 8 "0.6 C (L) . * 0,2 0.5: f ) f I. . . " -0.5 500. 1000 1500 2000 2500 Frequency (Hz) 3000 3500 4000, 500 •1000. 1500 2000 2500 Frequency (Hz) :3000. 3500 4000 Figure 4.4: Prediction of the reflection coefficient of a 6 mm single glass plate, at 75°. Prediction with plate-theory: —. Prediction with transfer-matrix model: — . Real part is shown in A, imaginary part in B. Note that the curves are almost superimposed, except at the coincidence frequency. 63 500' 1000 1500 2000 2500 Frequency (Hz) 3000 3500 4000: 1,000 500' 0; -500 •b N m -1000 -2000 k i i i — - 1 :l 1 B 6=45° -i i i 6=75° _ i i i 500 1000 1500 2000 2500 Frequency (Hz) 3000 3500 4000: Figure 4.5: Surface impedance of a single glass plate, 6 mm thick, predicted with the transfer-matrix model. Real parts are shown in A, imaginary parts in B. Impedance is normalized with respect to the characteristic impedance of air, 415 rayls. 4.1.2 Double plate A common room surface for an interior wall consists of two drywall panels fixed to either side of 2" x 4" (50 x 100 mm) studs. This surface was modeled using transfer matrices as a solid layer, a fluid layer and another solid layer. Each solid layer represents identical drywall panels, with properties as in Table 4.1. The fluid layer is 100 mm thick, but the studs cannot be taken into account in the modeling. Nevertheless, the transfer-matrix method is effective at predicting the sound propagation through the layers that make up this surface. As with the single glass-plate example, the transfer-matrix prediction of the reflection coefficient was compared to plate-theory predictions. The reflection coefficients predicted using plate-theory are shown in Figure 4.6 for the three angles, to show the acoustical behaviour of the double plate. The transfer-matrix and plate-theory predictions are shown together, for each of the three angles, in Figures 4.7-4.9. 64 As for the single plate, the real part of the reflection coefficient approaches one and the imaginary part is near zero. This is expected for a hard surface, as most of the sound energy will be reflected. The double-plate structure, however, exhibits a mass-air-mass resonance at a particular frequency that is angle dependent. In addition, coincidence effects also occur with double plates. The mass-air-mass resonance and coincidence phenomena are both evident in Figure 4.6. The mass-air-mass resonance is a low-frequency phenomenon, indicated by the first series of sharp peaks in both real and imaginary parts. The real part decreases from one to almost zero, and the imaginary part fluctuates rapidly from 0 to -0.5 to 0.5 and back to 0. This indicates almost total sound transmission, with a 180° phase change, through a double panel at the mass-air-mass resonance frequency. It is important to note, from Figure 4.6, that the mass-air-mass resonance frequency increases with increasing angle. Coincidence effects are noticeable in the double drywall panel, as with the single glass plate, at an angle of incidence of 75°. The coincidence frequency for smaller angles occurs at higher frequencies, not visible in Figure 4.6. o:5: $ o -0.5 '6= 75° t- 0°, 6- 45° 6 = 75° 500: 1000: .1500 .2000 2500 Frequency (Hz) :3000 3500 4000 6 £ 0 1 S -0-5 i i i i 6 - 75° j C J 1- 1 1. <; = 75= v. 0°. (, = 45° I I 1 1 r .• 500 : 1000i .1500 2000 2500 Frequency (Hz) :3000 3500 4000 Figure 4.6: Reflection coefficient of a double drywall panel, predicted with plate theory. The drywall plates are 12 mm thick, and the air gap is 100 mm. Real parts are shown in A, imaginary parts in B. For angles of 0°and 45°, the curves are almost superimposed. 65 The reflection coefficients shown in Figure 4.6 were compared with transfer-matrix predictions, with very good agreement. Comparison of plate-theory and transfer-matrix predictions for angles of 0°, 45°, and 75° are shown in Figures 4.7 - 4.9, respectively. As in the case of the single glass plate, a discrepancy occurred between the predictions of the two models at the coincidence frequency for an angle of 75°. As both models were developed theoretically, the differences in prediction are likely due to the assumptions made in the respective models, as pointed out in the single glass-plate case. The surface impedance for the double drywall plate, predicted with the transfer-matrix model, is shown in Figure 4.10. It is clear that the surface impedance varies with the angle of incidence and, hence, that this surface is of extended reaction. At low frequencies, the extreme points in the real and imaginary surface-impedance curves occur at the mass-air-mass resonance frequencies. These peaks corresponded to the peaks in the reflection coefficient curves in Figure 4.6, at the mass-air-mass resonance frequency; where the real part of the reflection coefficient approached zero, the real part of the surface impedance increased rapidly. Similar trends occur for the imaginary part of the surface impedance and reflection- coefficient, at the mass-air-mass resonance frequency. At the coincidence frequency for an angle of 75°, there is a very slight * 1 hr~— : — : — : — : — : — : — I 0.5-o i o-1 0 * A .0^ 1 1 1 1 _ _ l i 1 L 0 500 1000. 1500 20C0 2500 3000 3500 4000 Frequency (Hz) B • -~"'70 500 1000; 1500 2000: 2500 3000 3500 4000 Frequency (Hr) Figure 4.7: Prediction of the reflection coefficient of a double drywall plate with air gap, at normal incidence. Each drywall plate is 12 mm thick, and the air gap is 100 mm thick. Prediction with plate theory: —. Prediction with transfer-matrix model: — . Real part is shown in A, imaginary part in B. Note that the curves are superimposed. 4 0.2 ? 5 -0:2 66 1000 1500 .2000 2500 Frequency (Hz) 4000 1500 2000 2500 Frequency (Hz) Figure 4.8: Prediction of the reflection coefficient of a double drywall plate with air gap, at 45° angle of incidence. Each drywall plate is 12 mm thick, and the air-gap is 100 mm thick. Prediction with plate theory: —. Prediction with transfer-matrix model: — . Real part is shown in A, imaginary part in B. Note that the curves are almost superimposed. X "1 I T 500 1000, 1500 2000 2500 Frequency (Hz) 4000. i ° 5B "0:5 500 1000, 1500 2000 2500 Frequency (Hz) 4000 Figure 4.9: Prediction of the reflection coefficient of a double drywall plate with air-gap, at 75° angle of incidence. Each drywall plate is 12 mm thick, and the air gap is 100 mm thick. Prediction with plate theory: —. Prediction with transfer-matrix model: — . Real part is shown in A, imaginary part in B. Curves are superimposed except at the coincidence frequency. 67 p — 6 = 0°, i, = 45" Frequency (Hz) Figure 4.10: Surface impedance of a double drywall panel with air-gap, predicted with the transfer matrix model. Each panel is 12 mm thick, and the air-gap is 100 mm thick. Real parts are shown in A, imaginary parts in B. Impedance is normalized with respect to the characteristic impedance of air, 415 rayls. peak in the real part of the surface impedance, negligible compared to the peaks at the mass-air-mass resonance frequency. The imaginary part appears to be angle-independent between the peaks at the mass-air-mass resonance frequency and about 1000 Hz. Above 1000 Hz, however, there is a large variation in the normalized impedance. At normal incidence, the imaginary part of the surface impedance increases at a constant rate to a large, positive value. At near-grazing incidence, the imaginary part of the impedance decreases at a more rapid rate to a large, negative value. The peaks in the imaginary-part curves above 1000 Hz indicate phase changes between the pressure and particle velocity, and occur at higher resonance frequencies of the double plate. 4.1.3 Porous Material The most interesting feature of the transfer-matrix model is its application to modeling layers of porous materials. Porous materials such as glass-fibre and foam are used extensively in 68 rooms and sound enclosures to absorb sound and reduce noise levels to acceptable levels. Common materials found in rooms, such as carpet, acoustical tiles and drapery, are also sound absorbing, and can be considered porous. Therefore, the accurate modeling of such materials is important for a robust room-prediction model to be used to evaluate a room or noise enclosure acoustically. The common assumption regarding porous materials is local reaction; however measurements have been performed 3 7 that concluded that the surface impedance varies with angle of incidence for many porous materials and, hence, that they are extended-reaction materials. Two surface configurations involving single and double layers of porous materials were used for validation cases. For these materials, Biot parameters and measurement data for the surface impedance at normal incidence were obtained from literature. The first surface was a single layer of glass-fibre with a rigid backing, illustrated schematically in Figure 4.11, studied by A l l a r d 2 6 ' 4 2 . The second surface was a carpeted floor, modeled as a double porous-layer with rigid-backing, studied by Brouard et. al.43. The configuration of the double porous-layer is similar to Figure 4.11, except that there are two porous layers instead of one. For both surfaces, the validation was performed by comparing transfer-matrix prediction to the measurements of the surface impedance at normal incidence. Other than the transfer-matrix model, the only known theoretical model applicable to elastic-porous materials is the finite-element method 2 8 . A s this model was not available, the measurement data was chosen to incident ^ plane wave air \ 6 A \ e / reflected / plane wave 3 incident porous Biot waves layer 9 \ 1 3 reflected I Biot waves ////////////////// rigid backing X 3 Figure 4.11: Surface configuration of a porous layer with rigid backing. 69 compare to transfer-matrix predictions. A rigid-frame porous model developed by Allard-Attenborough26 was also used to compare surface impedance predictions at normal incidence. /. Glass-fibre layer with rigid backing The surface impedance at normal incidence is shown for the glass-fibre case in Figure 4.12. The frequency range of 350 - 1400 Hz was used to be consistent with the frequency range of Allard's measurements. The measurements and transfer-matrix prediction show reasonable agreement, particularly in predicting the dip at 800 Hz, and the peak at 850 Hz, in both the real and imaginary parts, respectively. This local minimum/maximum of the surface impedance was identified by Allard as resonance in the frame of the porous material. The frame resonance is possible from the partial coupling that occurs between air-borne and frame-borne waves that propagate through the fluid and solid parts, respectively, of the porous material. This effect was predicted using Biot theory for the porous material, as implemented in the transfer matrix for such a material, because the frame is considered elastic. From Figure 4.12, the surface impedance predicted at normal incidence using the transfer-matrix model and the rigid-frame model, can be seen to agree well at low and high frequency. However, the real component of the transfer-matrix prediction deviates significantly from the rigid-frame prediction in the 400 -1200 Hz range, and the imaginary component of the transfer-matrix prediction departs from the rigid-frame prediction in the 800 - 1200 Hz range. The transfer-matrix prediction of the surface impedance has a real-part minimum at 800 Hz, and an imaginary-part maximum at about 850 Hz. The rigid-frame model ignores this peak in the surface impedance prediction; however, it is predicted by Biot theory, as the transfer-matrix model results confirm. The frame-resonance behaviour in this material is confirmed with the experimental results shown. This example illustrates the effectiveness of using Biot theory for predicting the surface impedance of porous materials. With the transfer-matrix-model prediction of the surface impedance at normal incidence shown to be accurate, the next step was to consider predictions at oblique incidence. Surface impedance predictions by the transfer-matrix model for angles of 0°, 45°, and 85° are shown in Figure 4.13. The 85° angle was chosen to best represent the results for the range of higher angles in this case. The multiple curves in Figure 4.13 are clearly not superimposed. Therefore, this 70 Figure 4.12: Surface impedance of a single layer of glass-fibre at normal incidence. Material properties are listed in Table 4.2, thickness is 56 mm. Real part is shown in A, and imaginary part is shown in B. Transfer-matrix prediction: —. Rigid-frame model prediction: —. Measurements reproduced from Allard1: xxx. Impedance is normalized with respect to the characteristic impedance of air, 415 rayls. material is of extended reaction because the surface impedance varies with the angle of incidence. The frame-resonance frequency of this material increases with increasing angle of incidence, varying from 800 - 850 Hz for the real part and 850 - 900 Hz for the imaginary part. The surface impedance at the frame-resonance frequency is smallest at 0° and increases with increasing angle of incidence. For the real part, the normalized surface impedance at the frame-resonance frequency ranges from about 1.24 to 1.44, and from -0.7 to -0.6 for the imaginary part. The occurrence of a frame resonance for a porous material modeled as an elastic-porous layer is expected because the frame has mass and stiffness. Furthermore, for a coupled or partially-coupled porous material, the frame behaves as a fluid-loaded elastic-solid. The reflection coefficient was also predicted using the transfer-matrix model and is shown in Figure 4.14 for angles of incidences of 0°, 45°, 75°, and 85°. The real part of the 71 Frequency (Hz) Figure 4.13: Surface impedance of a single layer of glass-fibre at various angles of incidence, using the transfer-matrix model. Material properties are listed in Table 4.2, thickness is 56 mm. Real part is shown in A, and imaginary part is shown in B. Impedance is normalized with respect to the characteristic impedance of air, 415 rayls reflection coefficient is greatest at normal incidence and decreases at all frequencies with increasing angle of incidence. There is also a slight minimum in the real part, at all angles of incidence, at the frame-resonance frequency. The imaginary part of the reflection coefficient peaks at the frame-resonance frequency and increases with increasing angle of incidence, but only at frequencies above the resonant frequency. At frequencies below the frame-resonance frequency, the imaginary part of the reflection coefficient decreases for angles of incidence from 0° to 75°, and then increases for greater angles. ii. Carpeted floor — double porous layer with rigid backing The carpeted floor was studied in a similar manner to the single glass-fibre porous layer. The transfer-matrix model was used to predict the surface impedance at normal incidence and 72 -1:1 J L 1 1 I | 200 400 600 800 1000 1200 1400 Frequency (Hz) P •e . •5 -0:3 *g - • • J I. :| 6 = 85° > 6 = 0° _ — — • — ~ 2.- ^*<Z^~J/ - (, = 45° - ^ ^ ^ ^ -(,= 75°' i i i B 51 il L 1 1 1 1 200 400 600 800 1000 1200 1400 Frequency (Hz) Figure 4.14: Reflection coefficient of a single layer of glass-fibre at various angles of incidence, using the transfer-matrix model. Material properties are listed in Table 4.2, thickness is 56 mm. Real part is shown in A, and imaginary part is shown in B. compared to measurements performed by Brouard et. al.43. Predictions of surface impedances at oblique angles were performed with the transfer-matrix model to identify the surface as of extended reaction and suitable for inclusion in the rooms studied in Chapter 5. Furthermore, at oblique incidence, the reflection coefficients were predicted, as they are directly applied to the beam-tracing model. The surface impedance of the carpeted floor at normal incidence is shown in Figure 4.15. The impedance was predicted over the 50 - 4000 Hz range; however, measurement data was only available for 200 - 1300 Hz. The measurements are shown to agree well with prediction, for both real and imaginary parts. At almost all frequencies, the transfer-matrix model predicted slightly lower impedances than were measured. Measurements were also performed at frequencies above 4 kHz, and agreed well with the transfer-matrix predictions; however, they are outside the frequency range of interest here, and are not shown. 73 5 f --35:(f- =" -40 ;- ' <- c L — 1 1 1 1 0 500 1000 1500 2000 2500 3000 3500 4000 Frequency (Hz) Figure 4.15: Prediction vs. measurement of the surface impedance at normal incidence, of a carpeted floor. The carpet is modeled as two porous layers, each 3.5 mm thick, with properties listed in Table 4.2. Real part transfer-matrix prediction: -—. Real part measurements: Imaginary part transfer-matrix prediction: —. Imaginary part measurements: ooo. Impedance is normalized with respect to the characteristic impedance of air, 415 rayls. The surface impedance of the carpeted floor for angles of incidence of 0°, 45°, and 75° is shown in Figure 4.16. The impedance scale in Figure 4.15 has been expanded in Figure 4.16 for clarity. The real part of the normalized impedance is shown to increase by a factor of 2 from an angle of 0° to 75° up to about 1000 Hz. Above 1000 Hz, the rate of increase in impedance drops between 0° and 45°, but remains constant from 45° to 75°. The imaginary part of the impedance is dominant below 500 Hz, with little change between 0° and 75°. Above 500 Hz, there is a small difference in the impedance between angles of 0° and 45°, and a more rapid increase from 45° to 75°. Below 500 Hz, the much more dominant imaginary part varies little with angle of incidence and the local-reaction assumption for this surface is good. However, above 500 Hz, the real and imaginary parts of the impedance are closer in magnitude and vary significantly with 74 0 500 : 1000, 1500 2000 Frequency (Hz) 2500 3000 3500 4000: .0' 6= 75° -.10: 500 1500 2000 2500 Frequency (Hz) 3000 3500 4000 Figure 4.16: Prediction of the surface impedance at oblique incidence, for carpeted floor. Predictions were performed with the transfer-matrix model. The carpet is modeled as two porous layers, each 3.5 mm thick, with properties listed in Table 4.2. Real part is shown in A, and imaginary part is shown in B. Impedance is normalized with respect to the characteristic impedance of air, 415 rayls. angle. Therefore, this surface is indeed of extended reaction at frequencies above 500 Hz and is a good example of a real room surface that will be applied to the rooms studied in Chapter 5. The reflection coefficients of the carpeted floor predicted by the transfer-matrix model at angles of 0°, 45°, and 75°, are shown in Figure 4.17. Upon inspection of the curves in Figure 4.17, the real part exhibits a general decrease - more rapid at greater angles - with increasing frequency. This corresponds to typical porous materials for which the absorption generally increases with increasing frequency. The imaginary part also decreases with increasing frequency but, again, at different rates for different angles. For this material, at frequencies below 1300 Hz, the imaginary part decreases at a more rapid rate with increasing angle. However, at frequencies above 1300 Hz, the rate of decrease of the imaginary part decreases with increasing angle, with a slight increase occurring at an angle of 75°. The strong relative contribution of the imaginary part at high frequency suggests a phase change in the reflected sound wave at this surface. 75 0 500 1000 1500 2000 2500 3000 3500 4000 Frequency (Hz) Figure 4.17: Prediction of the reflection coefficient at oblique incidence, for carpeted floor. Predictions were performed with the transfer-matrix model. The carpet is modeled as two porous layers, each 3.5 mm thick, with properties listed in Table 4.2. Real part is shown in A, and imaginary parts are shown in B. In summary, the transfer-matrix model was tested on four different surfaces commonly found in real rooms. The tests results were confirmed by plate-theory predictions of the reflection coefficient in cases involving a single glass plate and a double drywall panel with air gap. With no other theoretical model available for the two porous-material cases, the transfer-matrix tests were compared to experiment. In these cases, the surface impedance at normal incidence was compared. In all cases, the transfer-matrix results agreed well with the comparative theory or measurements. The surface impedance at oblique incidence was predicted for each surface using the transfer-matrix model, and a variation with incident angle was shown. Therefore, these surfaces were of extended reaction and suitable to use in the rooms studied in Chapter 5. 76 4.2 Validation of the beam-tracing model The beam-tracing model calculated the steady-state sound-pressure level at a receiver point in a room, accounting for phase. The applicability of this model was tested by comparing it to another phase model, valid over the same frequency range as that used in the transfer-matrix model validation. The phase model used to compare predictions was a method-of-images model that included phase. Guo and Hodgson44 developed this model for active noise-control applications in rooms. The image-phase model was chosen for its availability and applicability to the same rooms and frequency range as the beam-tracing model. The validation of the beam-tracing model was a two-step process. The first step involved determining the number of beams, and a reflection order, that could predict the steady-state sound-pressure level accurately. The next step involved testing the beam-tracing model in the cases of three rectangular-shaped rooms, with angularly independent absorption on all surfaces. The three rooms chosen are detailed in Table 4.3. The dimensions and source and receiver positions of the three rooms were chosen to represent a variety of common rooms that were of different rectangular shapes: a cubic, 3 m x 3 m x 3 m high r o o m that represents a small office; a 10mx3mx3m high room that represents a corridor; and alOmx 10ra><3m high room that represents a small industrial workroom. The source and receiver positions were chosen such that they remain at a constant height, and were located at half the room width in all room configurations. The distance between the source and receiver varied with the length of the room. In all tests, the sound power of the source was chosen arbitrarily as 80 dB. All tests were performed over the frequency range of 50 - 4000 Hz, with 10 Hz intervals. The first step of the validation process was performed for the cubic room with no absorption on its surfaces. This configuration corresponded to the worst-case scenario; the sound energy detected at the receiver point was expected to be greater than in the other configurations. Therefore, the errors from these predictions were expected to be greater than errors from predictions involving larger rooms and with absorption. The second step in the validation process was performed for the three rooms, assuming an average diffuse-field absorption coefficient of 0.1 for all surfaces. This was chosen as a typical value corresponding to the typical ambient absorption of surfaces in rooms. The energy-absorption coefficient is related to the pressure-reflection coefficient as follows: 77 a(e) = l-\R(0f 4.1 The pressure reflection coefficient is obtained from the surface impedance, using Equation 3.31. The diffuse-field absorption coefficient is related to the surface impedance by1: a = Ja(6>) • sin(6>) • cos(0) dd 0 K = j(l - \R(0f )• sin(#) • cos(#) d6 4.2 if 1-Z-Zc/cos(8)V^ Z + Zc/cos(d) •sin(0)-cos(0)rf0 Solving the integral symbolically and then solving for the impedance in terms of the absorption, using Matlab, determined the surface impedance Z, in Equation 4.2. For a = 0.1, the surface impedance was calculated to be 29681 rayls. The beam-tracing program was executed for the three rooms, with an impedance of 29681 rayls defined for all surfaces. 4.2.1 Image-phase model An image-phase model was used to validate the beam-tracing model in both steps of the validation process. This model also calculated the steady-state sound-pressure level at a receiver point in a rectangular room, with one source. The pressure was calculated at the receiver point Table 4.3: Room configurations for the beam-tracing validation tests. The first room, with no absorption was used for beam number and reflection order tests. The comparisons involving all three rooms used an average Sabine absorption of 0.1 for all surfaces. Sound power of source is 80 dB, constant at all frequencies. Room Number Length (m) Width (m) Height (m) Source Position (m) Receiver Position (m) 1 3 3 3 {0.5, 1.5,2.0} {2.5, 1.5, 1.8} 2 10 3 3 {0.5, 1.5,2.0} {9.5, 1.5, 1.8} 3 10 10 3 {0.5,5.0, 2.0} {9.5,5.0, 1.8} 78 via an impedance transfer function between the source and receiver, and multiplied by the pressure amplitude constant, A, determined from Equation 3.38, as follows: p(r) = A-Zpp(rl,r0) 4.3 Surfaces in this model were assigned a single, constant surface impedance that is independent of frequency and angle. To model a surface with no absorption, an infinite impedance is specified. A surface defined by an average diffuse-field absorption coefficient is implemented in the same way as in the beam-tracing model, using Equation 4.2. 4.2.2 Number of beams and reflection order The reflection order required to predict the steady-state sound pressure with sufficient accuracy needed to be determined. The required reflection order is defined as the number of reflections of each beam at which the sound pressure build-up reaches a steady-state level, such that the level remained constant with increasing reflections. This was determined for each room using the beam-tracing model. The beam-tracing model was executed with 1280 beams, for 25 iterations of the reflection number ranging from 1-25. The reflection number sought was the one for which the predicted sound-pressure levels were equal for higher reflection numbers. If the levels did not reach steady-state after 25 reflections, the beam-tracing model was executed for higher reflection numbers until a suitable one was found. The beam-tracing predictions with reflection numbers 1-25, for the three rooms, are shown in Figure 4.18. These results indicate that the minimum numbers of reflections that are required for Rooms 1, 2 and 3 are 20, 21, and 21, respectively. In all rooms, for a number of reflections beyond 21, the predicted sound-pressure level was constant. By way of comparison, Guo and Hodgson44 used a reflection order of 21 in the image-phase model for rooms similar in size and surface absorption to the rooms studied here. It was decided to increase the number of reflections to 25, a round number that was a little extra conservative. The reflection order of 25 was used for all validation tests, for both models. This reflection order was also used in all of the room simulations performed in Chapter 5. The number of beams was determined by a series of comparisons of the beam-tracing and image-phase predictions, iterating over the icosahedron frequency. The icosahedron frequency, 79 ©120 5. 1110 g> 100 3.: W jS 90 ?• * 80 3 5 70 i i 10 15 T — r — T — r B 20 25 10 15 Number of reflections 20 25 Figure 4.18: Sound-pressure level predicted by the beam-tracing model with reflection numbers 1-25. A -Room 1, B = Room 2, C = Room 3. f, determines the number of beams, TVb, used to model the source, with JVb = 20x/ 2 . The image-phase model was used once to predict the steady-state sound-pressure level for Room 1 with no absorption on the surfaces. The beam-tracing model was applied to the same room for icosahedron frequencies ranging from 1 to 16, corresponding to a range of 20 to 5120 beams. 80 This range was chosen as it was expected that the number of beams required would be less than about 5000. However, the accuracy was tested further by executing the beam-tracing model with icosahedron frequencies of 23 and 32, corresponding to 10580 and 20480 beams. These numbers of beams were used to determine an upper limit on the relative difference between sound-pressure levels predicted by the two models. For each beam-tracing prediction, the percentage difference at each frequency between the two prediction methods was calculated: ALp =100%-^ p _image ^'p _beam ^ p _ image 4.4 The percentage differences were averaged over all frequencies and used to compare the predictions by the two models, as shown in Figure 4.19. The average percentage differences were calculated arithmetically from the decibel values of sound-pressure level. With the computer run-time directly proportional to the number of beams used in a simulation, it was desired to use the least number of beams, while maintaining reasonable accuracy, when predicting with the image-phase model. The results in Figure 4.19 show that the percentage difference between the predictions by the two models decay exponentially. A very large number of beams is required to get less than a 3% difference in decibel levels between beam-tracing and image-phase model predictions. Again, the disadvantage of using a large number of beams is the computation run-time of the beam-tracing program. The run-time for the beam-tracing model was almost linearly related to the number of beams, ranging from 2 minutes for 20 beams to about 8 Vz hours for 5120 beams, on a Windows-based PC with a 200MHz Pentium processor. The compromise choice was made of an icosahedron frequency of 8, corresponding to 1280 beams. This resulted in an average 5% difference between the predictions by the two models, and a reasonable computer run-time. It is expected that the steady-state sound-pressure levels predicted by the beam-tracing and image-phase models would not agree exactly when using a finite number of beams. This is because of the inherent error in the calculation of the non-exact source-receiver path in the beam-tracing model (see Section 3.2.7). A marginal improvement, to get a 3% difference between predictions by the two models, can be achieved by increasing the number of beams to 5120. Quadrupling that number of beams to 20480 beams achieves only a 1.7% difference. An infinite 81 20 A 8 18 V 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 icosahedron frequency, f B C D E F G 32 H 23 ... Figure 4.19: The average percentage difference between steady-state sound pressure-level predictions made by the beam-tracing and image-phase models. The room was cubic 3mx3mx3m with no absorption (R = 1) for all surfaces. Reflection order was 25 for all tests. Average percent differences were performed arithmetically using decibel values. The number of beams, Nb, is related to the icosahedron frequency, / , as Nb = 20*/2. A = 20 beams, B = S0 beams, C = 180 beams, D = 1280 beams, E = 2000 beams, F = 5120 beams, G = 10580 beams and H= 20480 beams. number of beams would be required for predictions by the beam-tracing model to compare exactly to the image-phase-model predictions. Evaluated for the worst-case scenario, the beam-tracing model was used with 1280 beams to model the source, and a reflection order of 25, in all future tests performed for this research. The next step in validating the beam-tracing model involved comparing predictions by the two models for the three test rooms. 4.2.3 Comparison of the beam-tracing and image-phase models for selected rooms The beam-tracing and image-phase models were applied to the three test rooms described in Table 4.3 to compare their predictions. The average percentage differences are listed in Table 4.4. It was surprising to see an average percent difference greater than the average difference 82 Table 4.4: Average percentage difference in comparing predictions by the image-phase and beam-tracing models. The percent difference is calculated arithmetically from the decibel values using Equation 4.4. The beam-tracing predictions were performed with 1280 and 2000 beams. Room No. Average % difference in prediction (1280 beams used in beam-tracing prediction) Average % difference in prediction (2000 beams used in beam-tracing prediction) 1 13.5 3.7 2 5.6 4.7 3 8.6 5.4 shown in Figure 4.18. The average differences were as high as 13.5% for the cubic room, a value that was considered too high. To improve the beam-tracing-model predictions for these three rooms, the number of beams was increased in two steps, to 2000 beams. The average percentage difference between predictions by the two models is shown in Table 4.4, with 2000 beams used in the beam-tracing model. These values were more satisfactory, and future beam-tracing predictions in this research used 2000 beams. A reflection order of 25 was still used. Run times increased to approximately 20 minutes. A comparison of steady-state sound-pressure-level predictions by the two models is shown in Figures 4.20 - 4.22, using 2000 beams in the beam-tracing model. Comparison of the predicted steady-state sound-pressure level by the two models shows good overall agreement. Generally, all curves are very irregular in shape with many local minima and maxima. Since phase models performed these predictions, it was expected that the peaks and dips in the curves represent excited modes in the room. The modal density of a room, however, is too dense to excite individual modes. These local minima and maxima are actually a close grouping of many excited modes. The maxima correspond to a grouping of modes that are similar in phase, the minima are a grouping of modes that are nearly opposite in phase. There are noticeable differences in the predictions by the two models at some of these points for all rooms. The beam-tracing model appears to under-estimate the sound energy predicted by the image-phase model at some of the peaks and dips, particularly at high frequency. Furthermore, there appear to be some minima predicted by the beam-tracing prediction that are not predicted by the image-phase model. Again, this appears to occur mostly at high frequency (above 2000 Hz) although, for the corridor example (room 2), it is noticeable at low frequency. 83 Figure 4.20: Comparison of steady-state sound-pressure-level predictions by the image-phase and beam-tracing model, for Room 1. The beam-tracing model was used with 2000 beams. The reflection order for both models was 25. Prediction with image-phase model: —. Prediction with beam-tracing model: —. The slight under-statement of sound energy, compared to the image-phase model, predicted by beam tracing can be explained by the differences in calculating the path length. The source-receiver paths calculated by the image-phase model are exact. The source-receiver paths calculated by the beam-tracing model are, however, approximate. As presented in Chapter 3, the source-receiver path calculated by the beam-tracing model is the path from the last image source to a point on a plane that contains the receiver point. Despite the differences in predictions by the two models, the majority of the peaks and dips were close for the two models and occurred at the same frequency. The beam-tracing-model predictions were considered satisfactory, and the model was applied with confidence, to meet the final research objective. 84 90 80 70 60 50 40 - .1 1, 1 1 1 , 1 1 1: I I 01 8! w 8-a 'ti' C : 3 O CO 80 70 60; 50 40 3000 3100 3200 3300 3400 3500 3600 Frequency (Hz) 3700 3800 3900 1000 2000 ..i 2100 2200 2300 2400 2500 2600 2700 2800 2900 3000 4000 Figure 4.21: Comparison of steady-state sound-pressure level predictions by the image-phase and beam-tracing model, for Room 2. The beam-tracing model was used with 2000 beams. The reflection order for both models was 25. Prediction with image-phase model: —. Prediction with beam-tracing model: —. In summary, both the transfer-matrix and beam-tracing models were shown to produce expected results for a series of cases involving real surfaces and rooms. The transfer-matrix model was shown to predict the reflection coefficient and surface impedance of common room surfaces - a sheet of glass, a double drywall-panel with air-gap, a fibre-glass material against a rigid, hard surface, and a carpeted floor - accurately. The number of beams and the reflection order required for the beam-tracing model to predict the steady-state sound-pressure level in a small, cubic room were determined. The beam-tracing model was then applied to three common types of room shapes, and compared to predictions by an image-phase model. The beam-tracing model displayed good overall agreement with the image-phase model, but it was noticed that the 85 90 : r *! i. rr 1 1 i r r SO-TO -60 -50 -4Q I I ' - j , •••••• • •••• ••••••l'. •• ••• . . ••• ••• : i v •_ . . . il.- •••• • I • • •• • .....I 0 100: 200 300 400 500 600 700 800 900 1000 40 I . h ii U L J L J L L 1 2000 2100 2200 2300 2400 2500 2600 2700 2800 2900 3000 50 V U 40 . L 1 i - J L 1 ~L J 1 L 3000 3100 3200 3300 3400 3500 3600 3700 .3800 3900 4000 Frequency (Hz) Figure 4.22: Comparison of steady-state sound-pressure level predictions by the image-phase and beam-tracing model, for Room 3. The beam-tracing model was used with 2000 beams. The reflection order for both models was 25. Prediction with image-phase model: —. Prediction with beam-tracing model: —. beam-tracing model slightly under-predicted the sound energy. With confidence in applying both models to their respective validation cases, the combined transfer-matrix and beam-tracing model was applied to rooms with surfaces of extended and local reaction. These rooms had the shapes used in the beam-tracing validation, and the surfaces used in the transfer-matrix validation. The application of the combined model to these rooms enabled the study of the difference in predicting the steady-state sound-pressure level in rooms modeled with extended vs. local-reaction surfaces. In the next chapter, applications of the combined model are discussed to achieve the final research objective. 86 C H A P T E R 5 M o d e l A p p l i c a t i o n s In the last chapter, the fourth research objective was addressed by validating the transfer-matrix and beam-tracing models. The transfer-matrix model was applied to surfaces that are commonly found in rooms: a single glass plate, double drywall plate, a carpeted floor and a fibre-glass layer. For each surface, the transfer-matrix prediction of the surface impedance and reflection coefficient was compared to prediction by other known theory and experiment, with very good agreement. The validation process for the beam-tracing model consisted of determining the reflection order and the number of beams to use, for three different room configurations. The minimum reflection order was found to be 20 reflections, at which the sound-pressure level approached steady-state. The accuracy of the beam-tracing model is defined by the number of beams used to model the source. When 2000 beams were used, decibel-level predictions made by an image-phase model and the beam-tracing model differed by an average of 5%. This difference was attributed to the non-exact source-receiver paths, and to using a single reflecting plane for an incident beam on intersecting surfaces, in the beam-tracing model. At this point, the new room-prediction model had been developed, implemented and validated, and was ready to be applied to achieve the final research objective. The new room-prediction model was applied to three room configurations to study the effect that an extended- vs. local-reaction room surface has on the predicted steady-state sound-pressure level in the rooms. One or more of the room surfaces were designated as a test surface. A designated test surface was a multi-layered structure that approximated a real room surface. Furthermore, the test surface was modeled in two ways in each room configuration - once as an extended-reaction surface, and the second time as a local-reaction surface. For each room 87 configuration, the sound-pressure level was predicted for each test-surface implementation. The sound-pressure levels were compared in the frequency domain to analyze the significance of extended- vs. local-reaction room surfaces on predicting the steady-state sound field in the room. 5.1 Test configurations Three room configurations approximating a small office, a corridor and a small industrial workshop were studied. Details of the room geometries are given in Table 5.1, with accompanying diagrams in Figure 5.1. Note that Room 1 (the office) has length equal to width equal to height. Room 2 (the corridor) has length greater than height and width. Room 3 (the workshop) has height less than length and width. Four different test surfaces were applied to each of the three room configurations that were used in the beam-tracing-model validation. Thus, twelve room configurations were considered for the study of extended- vs. local-reaction room surfaces. These room shapes were chosen as they are suitable for this research and are commonly found in the real world. A cubic room is of interest because a true diffuse sound field can only be realized in a room of this shape. The most realistic cubic-shaped room was considered to have dimensions of3mx3mx3m;a small office was the most common type of room that was considered to share these dimensions. The dimensions of the other two rooms were modifications of the cubic-shaped room. With constant height for all three rooms, the other two room shapes were chosen to be 'long and narrow', and 'low and square'. The long dimension was chosen to be 10 m and, therefore, the dimensions were 10mx3mx3m and 10 m x 10 m x 3 m for Rooms 2 and 3, respectively. Two realistic rooms were associated with these dimensions - a corridor in the case of Room 2, and a small industrial workroom in the case of Room 3. All surfaces, with the exception of the double-steel panels and the suspended-acoustical ceiling, had been involved in the transfer-matrix-model validation discussed in Chapter 4. The material properties of the glass and drywall are given in Table 4.1; the porous-layer properties of the carpet are listed in Table 4.2. The fibre-glass material properties are given in Table 4.2, except that a thickness of 100 mm (4") was used here. The material properties of the steel plate are given in Table 5.2. The suspended-acoustical ceiling consisted of two layers with a rigid 88 (concrete) backing: the first layer was a 12 mm (1/2") thick layer of the fibre-glass described in Table 4.2, and the second layer was an air gap 457 mm (18") thick. An average diffuse-field absorption coefficient of 0.1 that did not vary with incident angle or frequency, was assigned to all non-designated test surfaces in the rooms. The test surfaces were chosen as surfaces relevant to the three types of rooms. The glass plate was assigned to one wall in each room; this corresponded to a large window in one of the walls for each room. The double drywall panel was used to model all four walls in Rooms 1 and 2. This type of surface is commonly found in offices (especially interior offices in a building) and corridors. Double steel panels were used to model the walls in Room 3, because an industrial workshop is commonly constructed with metal-cladding walls and ceiling (or roof). In this research, it was important to include some examples of porous materials, because of their particular acoustical properties and their importance in rooms. Carpet is a porous material, common in many rooms including offices and corridors, and was an ideal candidate as a floor surface in Rooms 1 and 2. The carpet was omitted from Room 3 because a workshop does not usually have a carpeted floor. It is very common for offices and corridors to have acoustical tiles suspended from a concrete ceiling. Thus, the suspended-acoustical ceiling was used in Rooms 1 and 2. Two types of test surfaces were applied to the ceiling in Room 3. A common type of ceiling in a workshop is metal paneling or concrete that is treated acoustically (or thermally) with sound-absorbing material. The first type was approximated by the double steel panels previously used as the walls in Room 3. The second ceiling type was modeled by a 100 mm (4") fibre-glass layer with a rigid backing. All of these materials were chosen for their relevance to their respective rooms and to demonstrate the functionality of the transfer-matrix model. All tests were performed with the common set of parameters, as follows: • 2000 beams to approximate the source. • Reflection order of 2 5. • Frequency domain: 50 - 4000 Hz with 10 Hz increments. • Source sound-power level was constant at 80 dB for all frequencies. • Source positioned at 0.5 m in front of, and at the centre of, the wall at x = 0, and at a constant height of 2 m. • Receiver positioned at 0.5 m in front of, and at the centre of the wall at x - L, and at a constant height of 1.8 m. 89 Table 5.1: Room configurations Room No. Test surface Location of test surface Dimensions [L, W, H] (m) Source position {SXJ Sy, Sz} Receiver position {Rx, Ry, Rz} 1 Single glass plate Wall at y = 0. 3 x 3 x 3 {0.5, 1.5,2.0} {2.5, 1.5, 1.8} Double drywall panels Walls, x = 0, x = 3,y = 0, and y = 3. Carpeted floor Floor at z = 0. Suspended-acoustical ceiling Ceiling at z = 3 2 Single glass plate Wall at y = 0. 10 x 3 x 3 {0.5, 1.5,2.0} {9.5, 1.5, 1.8} Double drywall panels Walls, x = 0, x = 10, y = 0, and y = 3. Carpeted floor Floor at z = 0. Suspended-acoustical ceiling Ceiling at z = 3. 3 Single glass plate Wall at y = 0. 10 x 10x3 {0.5, 5.0, 1.5} {9.5,5.0, 1.8} Double steel panels Walls, x = 0, x = 10, y = 0, and y = 10. Double steel panels Ceiling at z = 3. Fibre-glass with concrete backing Ceiling at z = 3. Table 5.2: Material properties of the steel plate used for the walls and ceiling in Room 3 Material Steel Thickness 1.6 mm (1/16") Density 2500 kg/m3 Young's modulus 210 GPa Poisson ratio 0.33 90 Room! Room 2 Room 3 y(m) x(m) Figure 5.1: Three room configurations used in the application of the room prediction model to study the effect of extended- vs. local-reaction surfaces on the predicted steady-state sound-pressure level. 91 5.2 Results The predicted variation of the steady-state sound-pressure level with frequency is shown for each room configuration, with surfaces modeled as extended and local reaction, as a set of two curves (narrow band with 10 Hz interval) and as a bar chart of octave-band and total levels. Throughout the discussion and analysis of results, the two sound-pressure levels predicted in each room configuration are referred to as extended- or local-reaction levels. Levels are indicated as extended-reaction when an extended-reaction surface impedance was used to calculate the reflection coefficient of the test surface; a local-reaction level indicates that a locally-reacting surface impedance was used to obtain the reflection coefficient for the same test surface. With reference to comparing decibel levels, it should be noted that a 3 dB increase in sound-pressure level is equivalent to a doubling of the sound energy. With respect to the human perception of sound, an increase of 10 dB in sound-pressure level is equivalent to an increase in the loudness by a factor of 2. A difference of 0.5 dB or less will be regarded as insignificant here. 5.2.1 Single glass plate in Rooms 1, 2, and 3 For all three rooms, a single glass plate was used as a test surface to approximate a large window in the wall at y = 0. The results of the predicted steady-state sound-pressure level in the rooms with extended- and local-reaction surfaces are shown in Figures 5.2-5.4 for Rooms 1 -3, respectively. The narrow-band results also indicate equal levels predicted with the extended-and local-reaction surfaces. With the exception of Room 2, between 2000 Hz to 2500 Hz, these curves appear to be almost superimposed. In the 2000 - 2500 Hz range in Room 2, the extended- and local-reaction levels differ by 3 - 5 dB at the local minima of each curve. In all three rooms, the octave-band results indicate an insignificant difference in levels predicted with extended- and local-reaction surfaces. 92 30 I i I I I J .... 1 .J L L I 3000 3100 3200 3300 3400 3500 36003700 3800 3900 4000: Frequency (Hz) 120 2-110 1 •S 100 3 W w aj 80 c I 70 90 63 125 250 500 1000 Octave Band (Hz) 2000 Total Figure 5.2: Frequency response of the steady-state sound-pressure level predicted in Room 1 with a single glass plate in the y = 0 plane; all other surfaces have an average absorption of 0.1. Narrow-band results shown on top, filtered octave-band results are shown at the bottom. The glass plate was modeled as extended and local reaction. Extended-reaction results: —, grey bar. Local-reaction results: —, white bar. Source position: {0.5,1.5,2.0}. Receiver position: {2.5,1.5,1.8}. Room dimensions: 3mx3mx3m. 93 0 100 '200 300 400 500 600 -700 800 900 1000 jg | L I I L j 1 J I L I 3000 3100 3200 3300 3400 3500 3600 3700 3800 3900 4000 Frequency (Hz) 120 -2-110 1 :£ 100 I 90 a> Q. •6 80 I 70 63 125 250 500 1000 Octave Band (Hz) 2000 Total Figure 5.3: Frequency response of the steady-state sound-pressure level predicted in Room 2 with a single glass plate in the y = 0 plane; all other surfaces have an average absorption of 0.1. Narrow-band results shown on top, filtered octave-band results are shown at the bottom. The glass plate was modeled as extended and local reaction. Extended-reaction results: —, grey bar. Local-reaction results: —, white bar. Source position: {0.5,1.5,2.0}. Receiver position: {9.5,1.5,1.8}. Room dimensions: 10 m x 3 m x 3 m. 94 90 r 1 r r —'—i r 1 ' r r 5 0 -• 40 I L : i' : iL : I L 1 : • I : 1 0 100 200 300 400 500 600 700 800 900 1000 90 r r 1 r r f r 40 I 1. i i L 1 1 J i _ — 1 . 1 3000 3100 3200 3300 3400 3500 3600 3700 3800 3900 4000 Frequency (Hz) 120 : | 110 100 s> 3 W on <n 8i 80 c 3 0 ; <0 70 63 125 250 500 1000 2000 Total Octave Band (Hz) Figure 5.4: Frequency response of the steady-state sound-pressure level predicted in Room 3 with a single glass plate in the y = 0 plane; all other surfaces have an average absorption of 0.1. Narrow-band results shown on top, filtered octave-band results are shown at the bottom. The glass plate was modeled as extended and local reaction. Extended-reaction results: —, grey bar. Local-reaction results: —, white bar. Source position: {0.5, 5.0, 2.0}. Receiver position: {9.5, 5.0,1.8}. Room dimensions: 10 m x 10 m x 3 m. 95 5.2.2 Double drywall paneled walls in Rooms 1 and 2 The small office and corridor had drywall-paneled walls that were approximated by a drywall-air-drywall layered structure, with no rigid backing. In both rooms, this test surface was applied to all four walls. The predicted sound-pressure levels for the test surface modeled as of extended and local reaction are shown in Figures 5.5 and 5.6 for Rooms 1 and 2, respectively. The narrow-band results in Room 1 indicate some differences in the minima and maxima of the two curves. In the 50 - 450 Hz range, the local-reaction levels are lower than the extended-reaction levels by 10 - 15 dB at some of the minima in each curve. Other noticeable differences occur in the 2300 - 2600 Hz range. In this range, the peaks and dips of both curves occur at similar frequencies, but the levels are different. The local-reaction level is mainly higher by 5 - 8 dB at some frequencies, but a minimum less than the extended-reaction level by 10 dB occurs at 2480 Hz. In the case of Room 2, there are fewer differences between the extended- and local-reaction levels in the narrow-band results. Most of the differences are evident in the 50 - 400 Hz range, and at higher frequencies between 2400 - 2900 Hz. In many instances, the differences appear at the minima and maxima in both curves, and occur at similar frequencies. There is a slight shift toward higher frequency, however, in the local-reaction levels predicted in the low-frequency range. In both rooms, the results show significant differences in the extended- and local-reaction levels at low frequency. In Room 1, the extended-reaction level was 9 dB higher than the local-reaction level in the 63 Hz octave band. However, the local-reaction level was higher than the extended-reaction level by 3 dB in the 125 Hz octave band. In the higher-frequency octave bands, the two levels showed insignificant differences. In Room 2, there was less difference between the two levels at low frequency. In the 63 Hz octave band, the local-reaction level is higher than the extended-reaction level by approximately 4 dB. The extended-reaction level is slightly higher than the local-reaction level by 1 dB in the 125 Hz octave band. There is less than a 0.5 dB difference in the levels in the 250 - 2000 Hz octave bands. 96 <D a B; 3 . 0 : 2000 2100 3000 2200 3100 3200 3300; 3400 3500; 3600 Frequency (Hz) 1000: 1900 2000 2300 2400 2500 2600 2700 2800 2900 3000 3700 3800 3900 4000 120 I 1 110 100 I gj 90 8-a •6 80 c I 70 63 125 250 500 1000 Octave Band (Hz) 2000 Total Figure 5.5: Frequency response of the steady-state sound-pressure level predicted in Room 1 with double drywall panels (with 100 mm air gap) for the walls; all other surfaces have an average absorption of 0.1. Narrow-band results shown on top, filtered octave-band results are shown at the bottom. The drywall panels were modeled as extended and local reaction. Extended-reaction results: —, grey bar. Local-reaction results: —, white bar. Source position: {0.5, 1.5, 2.0}. Receiver position: {2.5,1.5,1.8}. Room dimensions: 3mx3mx3m. 97 50 -401 1 1 1 ii : _ J , ii J 1. l;: I 0 100 200 300 4001 500 600; 700 800: 900 1000 2000 2100 2200 2300 2400 2500 2600 2700 2800 2900 3000 63 125 250 500 1000 2000 Total Octave Band (Hz) Figure 5.6: Frequency response of the steady-state sound-pressure level predicted in Room 2 with double drywall panels (with 100 mm air gap) for the walls; all other surfaces have an average absorption of 0.1. Narrow-band results shown on top, filtered octave-band results are shown at the bottom. The drywall panels were modeled as extended and local reaction. Extended-reaction results: —, grey bar. Local-reaction results: —, white bar. Source position: {0.5, 1.5, 2.0}. Receiver position: {9.5,1.5,1.8}. Room dimensions: 10 m x 3 m x 3 m. 98 5.2.3 Double steel paneled walls in Room 3 The two steel panels used as the walls in Room 3 are separated by a 100 mm (4") air gap. The pressure levels predicted for this room with the test surface modeled as of extended- and local-reaction are shown in Figure 5.7. In the narrow-band results, the extended- and local-reaction levels differ at low frequency, in the 50- 350 Hz range. In this range, the peaks and dips in the two curves do not occur at the same frequencies. In fact, the shapes of the extended- and local-reaction curves are almost the opposite of one another in the 50 - 200 Hz range. At frequencies above 400 Hz, the two curves are superimposed. The differences in the extended- and local-reaction levels in this room configuration are less than the differences in levels for the first two rooms with drywall panels. The extended-reaction level is higher than the local-reaction level by 3 dB in the 63 Hz octave band. The local-reaction levels are higher by 3 and 2 dB in the 125 and 250 Hz octave bands, respectively. At higher frequencies, the two levels differ insignificantly. 5.2.4 Carpeted floor in Rooms 1 and 2 The carpeted-floor results are shown in Figures 5.8 and 5.9 for Rooms 1 and 2, respectively. In both rooms, a 1-2 dB difference between the pressure levels predicted with the test surface modeled as of extended and local reaction is evident at high frequency. The narrow-band results show no differences between the extended- and local-reaction levels at frequencies less than 700 Hz, in both rooms. In Room 1, there are noticeable differences in the extended-and local-reaction levels at most frequencies in the 1000 - 4000 Hz range. The local-reaction level generally appears to be higher, with some of the peaks and dips occurring at the same frequency. The largest differences are noticed in the 3400 - 4000 Hz range. At some frequencies in this range, the two levels differ by 5 - 20 dB. Some of these differences occur at minima or maxima in the two curves, at the same frequency. The largest difference occurs at 3600 Hz, where the local-reaction level reaches a local maximum and the extended-reaction curve shows a minimum. The narrow-band results for Room 2 are similar to the narrow-band results for Room 1. The differences in the extended- and local-reaction levels occur mainly at high frequency, 99 50f-40 I L 1 L j '£ 1. L I I 0 100 200 300 400 500 600 700 800 900: 1000 01 40 • L J 1 L 1 I J L L I a 1000 1100 1200 1300 1400 1500 1600 1700 1800 1900 2000 50' 40 1 L 1 L L. 1 L J t L I 3000 3100 3200 3300 3400 3500 3600 3700 3800 3900 4000 Frequency (Hz) 120 I 1 •a 100 e 3 (fl (fl f I 70 80 r 110K 90 ^ 63 125 250 500 1000 C^taye Band (Hz) 2000 Total Figure 5.7: Frequency response of the steady-state sound-pressure level predicted in Room 3 with double steel panels (with 100 mm air gap) for the walls; all other surfaces have an average absorption of 0.1. Narrow-band results shown on top, filtered octave-band results are shown at the bottom. The steel panels were modeled as of extended and local reaction. Extended-reaction results: —, grey bar. Local-reaction results: —, white bar. Source position: {0.5, 5.0, 2.0}. Receiver position: {9.5, 5.0,1.8}. Room dimensions: 10 m x 10 m x 3 m. 100 I 1. •S 8; 3 Ui w S J : i 3 . o; 1000 100 90 80 70 60 50 1100 1200 1300 1400 1500 '1600 1700 1800 3000 3100 3200 3300 3400 3500 3600 Frequency (Hz) 1900 1000 2000 & A I; I h I, ;l 1. I . 1. 1 -I': I I: l> 'I 1 i 1 |.. 3000 3700 3800 3900 4000 120 110 I I •8 100 I 90 t a •A 80 c: I 70 63 125 250 500 1000 Octave Band (Hz) 2000 Total Figure 5.8: Frequency response of the steady-state sound-pressure level predicted in Room 1 with carpeted floor; all other surfaces have an average absorption of 0.1. Narrow-band results shown on top, filtered octave-band results are shown at the bottom. The carpet was modeled as of extended and local reaction. Extended-reaction results: —, grey bar. Local-reaction results: —, white bar. Source position: {0.5, 1.5, 2.0}. Receiver position: {2.5,1.5,1.8}. Room dimensions: 3mx3mx3m. 101 I '8-3; w. in: 8; 1 c 3 O W -3000 3100 3200 3300 3400 3500 3600 Frequency (Hz) 3700 3800 1000 2000 3000 3900 4000 120 I 1 110k * 100 & 90 a to to * I 70 80 63 125 250 500 1000 2000 Total Octave Band (Hz) Figure 5.9: Frequency response of the steady-state sound-pressure level predicted in Room 2 with carpeted floor; all other surfaces have an average absorption of 0.1. Narrow-band results shown on top, filtered octave-band results are shown at the bottom. The carpet was modeled as of extended and local reaction. Extended-reaction results: —, grey bar. Local-reaction results: —, white bar. Source position: {0.5, 1.5, 2.0}. Receiver position: {9.5,1.5,1.8}. Room dimensions: 10 m x 3 m x 3 m. 102 particularly from 1500 - 1600 Hz and above 2000 Hz. In these ranges, the local-reaction level is generally higher than the extended-reaction level. Several minima in the local-reaction level are, however, less than the extended-reaction level - by as much as 30 dB at 3390 Hz, for example. In Room 1, the local-reaction levels were higher than the extended-reaction levels by approximately 1 - 2 dB in the 500 Hz and higher octave bands. In Room 2, the local-reaction levels were 2 - 3 dB higher in the 500 Hz and higher octave bands. 5.2.5 Suspended acoustical tile ceiling in Rooms 1 and 2 The results for the suspended acoustical ceiling are shown in Figures 5.10 and 5.11 for Rooms 1 and 2, respectively. The largest differences in the pressure levels were observed in these room configurations. Sound-pressure-level predictions were much higher in both rooms with the test surface modeled as of extended reaction, at low frequency. In the narrow-band results for Room 1, differences are seen between the extended- and local-reaction levels throughout most of the frequency range. At several frequencies, there is a local maximum in the local-reaction curve and a local minimum in the extended-reaction curve, or vice-versa. In some cases, local minima and maxima occur at similar frequencies, but the levels differ by as much as 10 dB. In the octave-band results, the extended-reaction level was 10 dB higher than the local-reaction level in the 63 Hz octave band. In the 125 Hz octave-band, the extended-reaction level was higher than the local-reaction level by approximately 2 - 3 dB. The two levels were equal in the 250 Hz octave band; the local-reaction level was slightly higher in the higher-frequency octave bands. The narrow-band results for Room 2 yield larger differences between the extended- and local-reaction curves than were seen in the narrow-band results for Room 1. The differences between the two levels occur at many frequencies throughout the frequency range. The local-reaction levels are noticeably lower in the 50 - 150 Hz range. At several frequencies between 1750 and 2790 Hz, the local-reaction curve has local minima that are much lower than the extended-reaction levels - by as much as 25 dB at 2600 Hz, for example. At many frequencies, there are occurrences of a minimum in the local-reaction curve where there is a maximum in the extended-reaction curve, or vice-versa. Furthermore, there are local minima or maxima for both curves that occur at similar frequencies, but with different levels. In the octave-band results, the sound- pressure-level predictions followed the same trends as in Room 1. The extended-reaction 103 1 £ a-in w •a:-Q. 0} 100 90 80 70 60 50 I il I: i i If I I. I • • • i " - -i • .|-. • . . .. ,((. .. . •"::f 1000 2000 3000 3000 3100 3200 3300 3400 3500 3600 Frequency (Hz) 3700 3800 3900 4000 120 110 1 1 'M100 I I 90 •6 80 c I 70 63 125 250 500 1000 2000 Total Octave Band (Hz) Figure 5.10: Frequency response of the steady-state sound-pressure level predicted in Room 1 with acoustical tile suspended 457 mm from concrete ceiling; all other surfaces have an average absorption of 0.1. Narrow-band results shown on top, filtered octave-band results are shown at the bottom. The ceiling was modeled as of extended and local reaction. Extended-reaction results: —, grey bar. Local-reaction results: —, white bar. Source position: {0.5,1.5,2.0}. Receiver position: {2.5,1.5,1.8}. Room dimensions: 3mx3mx3m. 104 120 250 500 1000 Octave Band (Hz) 2000 Total Figure 5.11: Frequency response of the steady-state sound-pressure level predicted in Room 2 with acoustical tile suspended 457 mm from concrete ceiling; all other surfaces have an average absorption of 0.1. Narrow-band results shown on top, filtered octave-band results are shown at the bottom. The ceiling was modeled as of extended and local reaction. Extended-reaction results: —, grey bar. Local-reaction results: —, white bar. Source position: {0.5,1.5, 2.0}. Receiver position: {9.5,1.5,1.8}. Room dimensions: 10 m x 3 m x 3 m. 105 level was 15 dB higher than the local-reaction level in the 63 Hz octave band. In the 125 Hz octave band, the extended-reaction level was 5 dB higher than the local-reaction level. The local-reaction level is about 4 dB higher in the 500 Hz octave band. In the 250, 1000 and 2000 Hz octave bands, there is an insignificant difference between the two levels. 5.2.6 Double steel paneled ceiling in Room 3 The sound-pressure-level predictions for this configuration are shown in Figure 5.12. The narrow-band results show differences in the extended- and local-reaction levels at low frequency, in the 50 - 800 Hz range. In this range, the local-reaction levels are generally higher than the extended-reaction levels. There are also peaks and dips in both curves that occur at different frequencies. The broad-band results reveal higher levels with surfaces modeled as of local reaction - by 2 - 3 dB in all but the 250 Hz octave bands. These results were in contrast to those for the industrial-workroom configuration with the double steel paneled walls. In the previous configuration with steel walls, there were higher extended-reaction levels in the 63 Hz octave band, and equal levels in the octave bands above 250 Hz. Furthermore, the local-reaction levels changed very little between the two room configurations, but the extended-reaction levels fluctuated slightly. This suggests that the positioning of the extended-reaction test surface, the geometry of the room or the source-receiver position, influences the sound-pressure-level prediction. 5.2.7 Fibre-glass layer with concrete backing as the ceiling in Room 3 The results for this last configuration are shown in Figure 5.13. Both the narrow- and octave-band results indicate equal levels with extended- and local-reaction surfaces. It should be noted that this was the same material used in the suspended acoustical tiles; in the present configuration, however, the material was much thicker and there was no air gap beyond the fibre-glass layer. This suggests that the material thickness in porous layers is not a significant factor in predicting the sound-pressure level in rooms with surfaces modeled as of extended- vs. local-reaction. However, it appears that the inclusion of a large air-gap is significant. 106 I 1 s. 3 « c 8 90: 80 70 60 50 1000 90 1100 i 300 400 500 600 700 800 1200 1300 1400 1500 1600 1700 1800 900 1000 I: I I J I I;--I 'I l: I I I I 1900 2000 2000 2100 2200 2300 2400 2500 2600 2700 2800 2900 3000 90 r 3000 3100 3200 3300 3400 3500 3600 3700 3800 3900 4000 Frequency (Hz) 120 I 110 •2: * 100 •E 3 W : W a: I I 70 90 80 63 125 250 500 1000 Octave Band (Hz) 2000 Total Figure 5.12: Frequency response of the steady-state sound-pressure level predicted in Room 3 with double steel panels on the ceiling; all other surfaces have an average absorption of 0.1. Narrow-band results shown on top, filtered octave-band results are shown at the bottom. The ceiling was modeled as of extended and local reaction. Extended-reaction results: —, grey bar. Local-reaction results: —, white bar. Source position: {0.5, 5.0,2.0}. Receiver position: {9.5, 5.0,1.8}. Room dimensions: 10 m x 10 m x 3 m. 107 i to 3 In in « c 3 O to iooo; 1100 1200 1300 1400 1500 1600 1700 1800 1900 1000 2000 2000 2100 2200 2300 2400 2500 2600 2700 2800 2900 3000 3000 3100 3200 3300 3400 3500 3600 3700 3800 3900 4000 Frequency (Hz) 120 110 •8 100 I 90 t •6 80 | 70 63 125 250 500 1000 Octave Band (Hz) 2000 Total Figure 5.13: Frequency response of the steady-state sound-pressure level predicted in Room 3 with 100 mm thick fibre-glass layer lining a concrete ceiling; all other surfaces have an average absorption of 0.1. Narrow-band results shown on top, filtered octave-band results are shown at the bottom. The ceiling was modeled of extended and local reaction. Extended-reaction results: —, grey bar. Local-reaction results: —, white bar. Source position: {0.5, 5.0,2.0}. Receiver position: {9.5,5.0,1.8}. Room dimensions: 10 m x 10 m x 3 m. 108 5.3 Explanation of results The significant differences in the extended- and local-reaction levels need to be explained to achieve the final research objective. Why are there differences in the steady-state sound-pressure-level predictions in rooms with surfaces modeled as of extended vs. local reaction? What can these differences be attributed to? A closer look at the reflection coefficients of the test surfaces yielded the answers to these questions. An analysis of the pressure reflection coefficients of the test surfaces was required to explain the differences between the extended- and local-reaction levels. The beam-tracing algorithm requires the pressure reflection coefficient of a multi-layered test surface, as predicted by the transfer-matrix model. It was, therefore, appropriate to study the differences between the reflection coefficients of each test surface, modeled as of extended and local reaction. The analysis was performed by calculating the two reflection coefficients for each surface - one calculated with an extended-reaction surface impedance, and the other obtained from a local-reaction surface impedance. The extended- and local-reaction surface impedances were calculated using the transfer-matrix model, and the reflection coefficients were obtained from equation 3.31. The reflection coefficients were calculated at the centre frequency of each octave band, for angles of 0° - 90°, with 1° increments. A set of six reflection-coefficient curves is shown for each test surface: single glass plate, double drywall panel with 100 mm air gap, carpeted floor, suspended-acoustical ceiling, and fibre-glass layer with a concrete backing. In addition to the reflection-coefficient information, knowledge of the incident angles associated with the first-order reflection paths was required. The largest contribution to the sound pressure at the receiver is from the direct and first-order reflection paths. With the reflection coefficients predicted for the applicable incident angle, the influence of the test surface on the sound pressure at the receiver can be estimated from a knowledge of the incident angle of the first-order reflected-sound path. Therefore, the difference in the sound-pressure levels in rooms with extended- and local-reaction surfaces may be explained by the differences in the extended- and local-reaction coefficients at the angle of incidence of the first-order reflected-sound paths. 109 The incident angle of the first-order reflected-sound path is easily calculated for all three rooms. In each room there are three unique first-order reflected-sound paths; the source-receiver paths with the reflecting surface as the ceiling, the walls at y = 0 and y = Lw, and the floor. Therefore, there are nine angles in total; these are listed in Table 5.3. The calculation of the angles for Room 2 is shown with the aid of Figure 5.14, as follows: The angle B\ is the incident angle for the source-receiver path with reflection from the ceiling. From the triangle defined by the points R, A, and S' in Figure 5.14, 0\ is obtained as follows: 6>, = tan" (AR) f 9 1 = tan 1 {AS'J U.2j = 76.3° Similarly, the angles &i and 9$ are found from Figure 5.14 as follows: 6 2 = tan t93 = tan -1 yAS"j \SS" j = tan" = tan" -1 (9 -1 0 , = 67.1° = 71.6° The three incident angles are calculated for the other two rooms in the same way. In the following discussion, the term 'coefficient' implies the pressure reflection coefficient. Furthermore, the term 'extended-reaction coefficient' refers to the pressure reflection coefficient of a test surface that was obtained from the extended-reaction surface impedance of that surface. The term 'local-reaction coefficient' refers to the pressure reflection coefficient of a test surface that was obtained from a locally-reacting surface impedance. The term 'first angle of incidence' refers to the angle of incidence of the first reflection on a room surface. 110 Figure 5.14: Diagram of Room 2, in the x-z plane (top) and the x-y plane (bottom), to calculate the incident angles from the first-order reflected-sound path. Table 5.3: Incident angles from the first-order reflected-sound paths for the three rooms. Room No. 0i (ceiling) 0z (floor) 03 (wall) 1 42.3° 27.8° 33.7° 2 76.3° 67.1° 71.6° 3 76.3° 67.1° 42.0° 111 5.3.1 Single glass plate In all three rooms, with a single glass plate along the y = 0 wall, the sound-pressure levels were the same for the test surface modeled as of extended and local reaction. The reflection coefficient for the glass plate, calculated with an extended- and local-reaction surface impedance, is shown in Figure 5.15. These curves reveal why there was no appreciable difference in pressure levels due to the test surface modeled as of extended or local reaction. At each octave-band centre frequency, the extended- and local-reaction coefficients are almost identical in their real parts. There are some differences, however, in their imaginary parts; the local-reaction coefficient is greater than the extended-reaction coefficient at angles close to grazing incidence. At the 63 Hz and 125 Hz centre frequencies, the two curves are identical up to 75°, beyond which the local-reaction coefficient increases more rapidly toward a value of one. At 250 Hz and 500 Hz, the local-reaction coefficient increases more rapidly toward a value of one, for an incident angle greater than 85°. The two curves are almost identical at 1000 Hz. At 2000 Hz, the extended-reaction curve increases more rapidly towards a value of one, for incident angles greater than 60°. It is clear that the extended- and local-reaction coefficients are very similar at incident angles less than 75°, at all of the octave-band centre frequencies. From Table 5.3, the largest first angle of incidence on the wall is 76.3°. Therefore, it is not surprising that the difference in the test surface modeled as of extended and local reaction was not significant in predicting the sound field in any of the rooms. 5.3.2 Double-drywall panels The drywall panels in Rooms 1 and 2 showed a significant difference between the predicted sound-pressure levels with surfaces modeled as of extended and local reaction. The extended- and local-reaction coefficients, as a function of incident angle for the octave-band centre frequencies, are shown in Figure 5.16 for the drywall plates. Significant differences are seen in the reflection coefficient of this test surface modeled as of extended and local reaction. The real part of the extended-reaction coefficient is greater 112 Angle of Incidence (") Angle of Incidence (") Figure 5.15: Reflection coefficient of the single glass plate, as a function of incident angle. Curves are shown for each centre frequency of the six octave bands: 63 Hz, 125 Hz, 250 Hz, 500 Hz, 1000 Hz and 2000 Hz. Reflection coefficient obtained from extended-reaction surface impedance: —. Reflection coefficient obtained from local-reaction surface impedance: —. 113 Angle of Incidence (a) Angle of Incidence (") Figure 5.16: Reflection coefficient of the double-drywall panel, as a function of incident angle. Curves are shown for each centre frequency of the six octave bands: 63 Hz, 125 Hz, 250 Hz, 500 Hz, 1000 Hz and 2000 Hz. Reflection coefficient obtained from extended-reaction surface impedance: —. Reflection coefficient obtained from local-reaction surface impedance: —. 114 than the local-reaction coefficient at 63 Hz. At higher frequencies, the real parts of the extended-and local-reaction coefficients are very similar; there are, however, narrow minima that occur in the extended-reaction curve that do not appear in the local-reaction curve. At 125 Hz, the minimum occurs at an angle of incidence of 55°. With increasing frequency, these minima occur at greater angles of incidence. At 2000 Hz, the minimum in the extended-reaction curve occurs at near grazing incidence. The imaginary parts of the extended- and local-reaction coefficients are significantly different. The most dramatic difference occurs at 63 Hz, where the two curves appear to be the opposite of one another; the extended-reaction curve is mostly positive, and the local-reaction curve is mostly negative. For the octave-band centre frequencies of 125 Hz and higher, there appears to be a phase change predicted by the extended-reaction coefficient only. This is indicated by the sharp local minima and maxima in the extended-reaction curves; they are absent in the local-reaction curve. At a given frequency, the local minima/maxima appear at the same angle of incidence as the minima in the real part of the extended-reaction coefficient. The imaginary parts of the extended- and local-reaction coefficients are similar from normal incidence to the angle at which the phase change is indicated in the extended-reaction curve. Above this angle, however, the imaginary part of the local-reaction coefficient is greater than the imaginary part of the extended-reaction coefficient. It is not entirely clear that the first angle of incidence can account for the difference between the sound-pressure levels with the test surface modeled as of extended and local reaction. In Room 1, the extended-reaction level was higher by 9 dB in the 63 Hz octave band and the first incident angle of the wall was 33.7°. In room 2, the local-reaction level was higher by 4 dB in the 63 Hz octave band and the first incident angle of the wall was 71.6°. The result for Room 2 was unexpected, because the magnitude of the extended-reaction coefficient is greater than the local-reaction coefficient, particularly for increasing incident angle. Therefore, phase effects associated with the local-reaction coefficient must account for the 63 Hz octave-band result in Room 2. The phase effects cause the interference of sound waves that affect the resulting sound pressure; phase is also more significant at low frequency. The same is true in the 125 and 250 Hz octave bands; in Room 1 there were higher sound-pressure levels predicted with surfaces modeled as of local reaction - in Room 2 the extended-reaction level was slightly higher. In both rooms, the first angle of incidence of the wall indicates that the local-reaction 115 level should be equal to or greater than the extended-reaction level. As this is not the case, the difference in the levels is likely due to phase differences between the extended- and local-reaction coefficients. The fact that this occurs at low frequency further suggests that interference effects are the cause of the difference in levels. The extended- and local-reaction levels are almost equal at high frequency, where phase is less of an issue. It should also be noted that the fluid layer has a significant effect on the difference between the extended- and local-reaction reflection coefficients. A fluid layer is by no means of local reaction; except in the case of an incident wave at normal incidence, the wave(s) transmitted through a fluid layer occur at some angle 9, where the impedance is Zc/cos(9). It is quite possible, depending on the angle, that the impedance is very different from the impedance at normal incidence, Zc. Therefore, the fluid layer influences the difference between the extended- and local-reaction reflection coefficients significantly. 5.3.3 Carpetedfloor In Rooms 1 and 2 with carpeted floor, the sound-pressure level predicted with surfaces modeled as of local reaction was slightly higher than levels predicted with extended-reaction surfaces, at high frequency. The reflection coefficient of the carpeted floor is shown in Figure 5.17. The results seem to confirm the equal levels at low frequency, but indicate that the extended-reaction level should be higher at high frequency. Differences in the extended- and local-reaction coefficients appear at high frequency, predominantly in the imaginary part. The real and imaginary parts of the reflection-coefficient curves in Figure 5.17 are almost superimposed at the octave-band centre frequencies of 63 - 250 Hz. At higher frequencies, the imaginary part of the extended-reaction coefficient is noticeably greater than the imaginary part of the local-reaction coefficient. This difference increases with increasing frequency; the difference is greatest for an incident angle of 85° at 500 Hz, and greatest for an incident angle of 75° at 2000 Hz. The magnitudes of the extended- and local-reaction coefficients are almost equal in the 63 - 250 Hz range, which suggests the reason for the equal sound-pressure levels in both rooms with surfaces modeled as of extended and local reaction. At higher frequencies, however, the local-reaction levels were higher; the greatest difference was 3 dB in Room 2 in the 500 Hz octave band. The first angle of incidence of the floor in Room 2 is greater than it is in Room 1; it 116 0 15 30 45 60 75 90 Angle of Incidence (") c :oeffick 0 w c -0.5 •a o -1 0 OI tj a -0.5 (B Imagii -1 0 -0.5 -1 0 -0.5 0 15 30 45 60 75 90 Angle of Incidence (") Figure 5.17: Reflection coefficient of the carpet, as a function of incident angle. Curves are shown for each centre frequency of the six octave bands: 63 Hz, 125 Hz, 250 Hz, 500 Hz, 1000 Hz and 2000 Hz. Reflection coefficient obtained from extended-reaction surface impedance: —. Reflection coefficient obtained from local-reaction surface impedance: —. 117 also is within the range where the imaginary parts of the local- and extended-reaction coefficients differ the most at 500 Hz. In terms of magnitude, the extended-reaction coefficient is greater and it was expected that the extended-reaction levels would be higher. As they were not, it is possible that interference effects were responsible for the 500 Hz octave-band result in Room 2. The slightly higher local-reaction levels in both rooms at higher frequencies also suggests that interference effects are significant at high frequency for this surface. 5.3.4 Suspended-acoustical-tile ceiling The most significant difference between the extended- and local-reaction levels occurred with the suspended-acoustical-tile ceiling. The extended- and local-reaction reflection coefficients for this test surface are shown in Figure 5.18. It is not surprising, given the room-prediction results, that the two reflection coefficients are very different in both their real and imaginary parts, at all frequencies. The greatest difference in the room-prediction results appears at low frequency, particularly in the 63 Hz octave band. The extended-reaction level was higher by 10 and 15 dB in Rooms 1 and 2, respectively. At 63 Hz, the real part of the extended-reaction coefficient is much greater than the real part of the local-reaction coefficient, for all angles except normal incidence at which they are the same. The imaginary part of the extended-reaction coefficient is greater up to an angle of 80°, where it reaches a maximum and decreases toward -1 for increasing angle; the local-reaction coefficient continues to increase toward 0 for increasing angle. The extended-reaction level was also higher in the 125 Hz and 250 Hz octave bands; this is confirmed by a larger magnitude of the extended-reaction coefficient at these frequencies. At higher frequency, the extended-reaction coefficient is somewhat sinusoidal in its variation with incident angle. The local-reaction coefficient varies as a smooth monotonic curve with increasing incident angle; it follows the minima of the real part of the extended-reaction coefficient at 500 - 2000 Hz, and follows the maxima of the imaginary part of the extended-reaction curve at 1000 Hz. The local-reaction levels were equal to, or slightly higher than, the extended-reaction levels in both rooms at high frequency. It is difficult to relate the difference in levels in the higher-frequency octave bands to the reflection coefficients at the first angle of incidence. The first angles of incidence of the floor are 27.3° and 67.1° for Rooms 1 and 2, respectively. At these angles, the local-reaction coefficients appear to be slightly greater in 118 0 15 30 45 60 75 90 0 15 -30? 45 60 75 90 Angle of Incidence (a) Angle of Incidence (") Figure 5.18: Reflection coefficient of the suspended-acoustical tile, as a function of incident angle. Curves are shown for each centre frequency of the six octave bands: 63 Hz, 125 Hz, 250 Hz, 500 Hz, 1000 Hz and 2000 Hz. Reflection coefficient obtained from extended-reaction surface impedance: —. Reflection coefficient obtained from local-reaction surface impedance: —. 119 magnitude than the extended-reaction coefficients at frequencies of 500 - 2000 Hz. This suggests the slightly higher local-reaction levels at higher frequency in both rooms. However, the differences in levels due to interference effects should not be ignored. Furthermore, the extended-reaction coefficient appears very different at a single frequency; the comparison of single-frequency reflection coefficients in this case is difficult to relate to octave-band results for the sound-pressure level. As for the case of the double drywall plates, there is an air gap of significant thickness in this test surface. Since the local-reaction assumption for air is known to be very inaccurate, it is not surprising that the extended- and local-reaction coefficients are very different. Therefore, modeling this surface as of extended reaction is significant in predicting the sound-pressure level in a room, due to the air gap present. 5.3.5 Fibre-glass with concrete backing The industrial-workroom configuration with fibre-glass covering a concrete ceiling is an example of acoustical treatment in an industrial setting. The reflection coefficients for this surface are shown in Figure 5.19. The extended- and local-reflection coefficient curves in Figure 5.19 confirm the almost equal sound-pressure levels predicted with surfaces modeled as of extended and local reaction. The extended- and local-reflection coefficients are very similar at all frequencies and angles for the fibre-glass with concrete backing. A small difference is noted in the real part at 500 Hz where the local-reaction coefficient is slightly greater. This is likely due to the frame resonance that occurs near this frequency for this type of porous material. Small differences occur in the imaginary part at 500 Hz, 1000 Hz and 2000 Hz. At these frequencies, the imaginary part of the extended-reaction coefficient is greater than the imaginary part of the local-reaction coefficient; the largest difference occurs at 2000 Hz. These results confirm the influence of an air gap on the difference between the extended-and local-reaction coefficients. In the acoustical-tile surface, the same material was used as here, but with an air gap that was much thicker. The reflection coefficients of these test surfaces were very much different for surfaces modeled as of extended reaction. It is also interesting to note that the local-reaction coefficients of the acoustical tile (Figure 5.18) and fibre-glass layer (Figure 5.19) appear to be very similar, except for the imaginary part at high frequency. For the 120 0 15 30 45 60 75 90 Angle of Incidence (°) 0) o u o <s (J 0) o •fc CO a ra •8" 0 15 30 45 60 75 90 Angle of Incidence (") Figure 5.19: Reflection coefficient of the fibre-glass layer with a concrete backing, as a function of incident angle. Curves are shown for each centre frequency of the six octave bands: 63 Hz, 125 Hz, 250 Hz, 500 Hz, 1000 Hz and 2000 Hz. Reflection coefficient obtained from extended-reaction surface impedance: —. Reflection coefficient obtained from local-reaction surface impedance: —. 121 fibre-glass, the imaginary part of the local-reaction coefficient increases rapidly toward 0 at incident angles near grazing. For the acoustical tile, this occurs at 63 Hz and 125 Hz; at higher frequencies, the local-reaction coefficient increases slightly. 5.4 Summary The application of the room-prediction model that was developed in this research provided insight into the significance of predicting the sound-pressure level in rooms with surfaces modeled as of extended and local reaction. The following concluding remarks can be made with respect to predicting the sound-pressure level in rooms with surfaces modeled as of extended and local reaction: 1. The difference between modeling a room surface as of extended or local reaction is not significant when the surface is a single plate, or a single layer of material (solid or porous) with a rigid backing. 2. The difference between modeling a room surface as of extended or local reaction is significant when the surface consists of multi-layers of solid or porous material, and includes a layer of fluid with a large thickness relative to the other layers. 3. For surfaces in which the reflection coefficient is different due to using an extended-or local-reaction surface impedance, the extended- and local-reaction assumption may be significant depending on the source and receiver positions. This becomes significant when the positions are such that the first angle of incidence is near grazing. 122 C H A P T E R 6 Conclusions This thesis pertains to the prediction of sound fields in rooms. The focus of this research was to determine the sound field in rooms with an unprecedented modeling of the room surfaces. The concepts of the room response and the acoustical properties of surfaces were introduced. It was revealed that there is a lack of a room-prediction model that incorporates extended-reaction surfaces; i.e. surfaces for which the surface impedance varies with angle of incidence, as well as frequency. Existing room-prediction models assume that the surfaces are locally reacting. Therefore, the significance of extended- vs. local-reaction surfaces on the prediction of the room-response is unknown; this is directly related to the inaccuracy associated with the local-reaction assumption. These considerations led to the general research objective - to develop and apply a new room-prediction model, with phase and the ability to model surfaces as of extended and local reaction. The model was to be applied to rooms with surfaces modeled as of extended and local reaction, to determine the benefit of extended-reaction modeling on the predicted sound field. This was to be a unique room-prediction model, and its development required a review of literature pertaining to existing room-prediction models and to the prediction of the acoustical properties of surfaces. A transfer-matrix approach29'35 was adopted to predict of the acoustical properties of extended- or local-reaction surfaces. This model calculates the surface impedance and pressure reflection coefficient of multi-layered surfaces; these surfaces consist of a series of layers with finite thickness and infinite lateral extent, and materials classified as either fluid, elastic-solid or elastic-porous. Biot theory is used in the transfer-matrix formulation of the porous layer. Of the room-prediction models reviewed, the model chosen to base this research 123 on was a triangular beam-tracing approach proposed by Lewers1 5. The advantages of the beam-tracing model were its computational efficiency, applicability to rooms with arbitrary geometry, potentially high accuracy, ease of incorporating phase, and ability to integrate the transfer-matrix model for the surfaces. From the literature review, the specific research objectives were established: 1. Develop a transfer-matrix model for multi-layered surfaces with an air or rigid backing. Surface layers were assumed to have a finite thickness, with no edge constraints; they consisted of fluid, elastic-solid or elastic-porous materials. 2. Develop a triangular beam-tracing model with phase. This model would be valid for specular reflection only and for a single point source and point receiver; it would predict the steady-state sound-pressure level. 3. Implement the new room-prediction model by integrating the transfer-matrix model into the beam-tracing model. The transfer-matrix model is dedicated to predicting the pressure reflection coefficient of a surface. 4. Validate the new model using test cases chosen to represent real rooms and surfaces. 5. Apply the new model to determine the significance of modelling room surfaces as of extended reaction on the prediction of the steady-state sound-pressure level in a room. The first three objectives were accomplished as discussed in Chapter 3. The first part focused on a detailed development of the transfer-matrix model. The transfer-matrix formulations for fluid, elastic-solid and elastic-porous layers were demonstrated. Biot theory was introduced and used in the transfer-matrix formulation of the porous layer. The complete model was formulated as an assembly of the boundary conditions at each layer interface; this involved interface matrices and the transfer matrix of each layer. It was shown how the surface impedance and pressure reflection coefficient of a multi-layered surface, modeled as of extended or local reaction, was obtained from the transfer-matrix model. The second part of Chapter 3 focused on the development of the triangular beam-tracing model, with phase. The entire algorithm along with the implementation of phase was developed from first principles. The theoretical development involved the source and beam initialization, beam tracing, source-receiver path estimation and the calculation of the sound-pressure level at 124 the receiver. A spherical wave was approximated by a point source surrounded by an icosahedron with subdivided faces. A beam was defined by the point source as the vertex, and the three vertices of an icosahedron face. Each beam was propagated through the room by tracing its centre ray up to a specified reflection order, in an attempt to find a valid source-receiver path. The beam face represented a portion of the spherical sound wavefront as a complex pressure (with magnitude and phase). This was in contrast to the energy-based form of this model, for which only the sound energy is calculated at the beam front. With each surface reflection, the associated complex pressure reflection coefficient was multiplied by the incident beam's complex pressure to find the pressure at the reflected beam front. The sum of the complex pressures at the beam face for each occurrence of a valid source-receiver path is that beam's contribution to the sound pressure at the receiver. The sum of the pressure contributions from all beams yields the steady-state sound-pressure level at the receiver point. The model works in the frequency domain; the time response can be obtained using an Inverse Fast-Fourier Transform (IFFT). The two models were integrated to form the new room-prediction model required for this research. The models were programmed separately in Matlab code. The transfer-matrix model was implemented into a function that output the complex reflection coefficient and surface impedance. The inputs to this function were frequency, incident angle and the material properties of the layer. The function was called in the beam-tracing program at each occurrence of a surface reflection from a multi-layered surface, to calculate the complex reflection coefficient. Both models could be executed separately; the transfer-matrix model calculated the acoustical properties of surfaces, and the beam-tracing model was used for rooms with surfaces defined by a single or frequency-dependent impedance, reflection or absorption coefficient. The new model had to be validated before it could be applied to achieve the final research objective. Validation was the fourth research objective, addressed in Chapter 4. Since the combined model was unique, the beam-tracing and transfer-matrix models had to be validated separately. The transfer-matrix model was tested in the case of surfaces that are commonly found in rooms. The chosen test surfaces were: a single glass plate, double drywall panel, a carpeted floor (as a double porous layer with rigid backing), and a fibre-glass layer with rigid backing. The transfer-matrix predictions in these cases agreed well with the comparative theory (plate examples) and experiment (carpeted floor and fibre-glass examples). 125 The beam-tracing model was applied to three room configurations: a small office (3 m x 3 m x 3 m), a corridor (10 m x 3 m x 3 m), and a small industrial workroom (10 m x 10 m x 3 m). Two issues were involved in the validation of this model; the required reflection order and the number of beams had to be determined in order to predict the steady-state sound-pressure level with sufficient accuracy. The required reflection order was defined as the number of reflections for which the predicted sound-pressure level approached a constant value, for each room. It was found that a reflection order of 25 was suitable for predicting the steady-state sound pressure. The minimum number of beams was determined by comparing predictions made by the beam-tracing model and an image-phase model. A minimum of 2000 beams was needed for the sound-pressure level predicted by the beam-tracing model to be within 5% of the image-phase-model prediction. With both models yielding the expected results for the test cases, the new room-prediction model was applied with confidence to achieve the final research objective. The new room-prediction model was applied to various room configurations with surfaces modeled as of extended and local reaction. This application was an attempt to determine the significance of modeling a surface as of extended or local reaction on predicting the steady-state sound field in a room. The model was applied to twelve room configurations; four different test surfaces were considered in each of the three rooms used in the beam-tracing model validation. The steady-state sound-pressure level was predicted separately with the test surface modeled as of extended reaction and then with the test surface modeled as of local reaction, for each room configuration. The frequency domain was 50 - 2000 Hz, with 10 Hz increments. Comparison was made between the two levels, filtered into octave bands, for each room configuration. The results are summarized as follows: • All three rooms with a single glass plate along one of the walls showed no difference between levels predicted with surfaces modeled as of extended or local reaction. • With double drywall panels on the walls of Rooms 1 and 2, there were differences between the extended- and local-reaction levels at low frequency. In Room 1, the extended-reaction level was 9 dB higher in the 63 Hz octave band; the local-reaction level was 3 dB and ldB higher in the 125 Hz and 250 Hz octave bands, respectively. In the higher-frequency octave bands, the levels were nearly equal. In Room 2, the 126 local-reaction level was higher by 4 dB in the 63 Hz octave band. The levels were nearly equal in the higher-frequency octave bands. • With double steel plates on the walls of Room 3, the extended-reaction level was higher by 2 dB in the 63 Hz octave band. The local-reaction level was higher by 3 dB and 2 dB in the 125 Hz and 250 Hz octave bands, respectively. Both levels were nearly equal in the higher-frequency octave bands. • With carpet covering the floor of Rooms 1 and 2, local-reaction levels were 1 - 3 dB higher in the 500 Hz - 2000 Hz frequency octave bands. • The most significant difference between extended- and local-reaction levels was found in Rooms 1 and 2 with acoustical tile suspended from a concrete ceiling. In both rooms, the extended-reaction levels were higher at low frequency and the local-reaction levels were slightly higher at high frequency. In Room 1, the extended-reaction level was higher by 10 dB and 3 dB in the 63 Hz and 125 Hz octave band, respectively. In Room 2, the extended-reaction level was higher by 15 dB and 5 dB in the 63 Hz and 125 Hz octave band, respectively. The local-reaction level was higher by 4 dB in the 500 Hz octave band. • In Room 3 with a double steel paneled ceiling, the local-reaction level was higher by 2 - 3 dB in the 63 Hz, 125 Hz, and 500Hz - 2000 Hz octave bands. The extended-reaction level was slightly higher in the 250 Hz octave band. • In Room 3 with a fibre-glass-lined concrete ceiling, there was very little difference between the extended- and local-reaction levels in all octave bands. The results were explained by an analysis of the reflection coefficients of the test surfaces. For each test surface, two reflection coefficients were predicted for incident angles of 0° - 90°, with a 1° increment, at each octave-band centre frequency. One of the reflection coefficients was predicted with the extended-reaction surface impedance, the other was obtained with the local-reaction surface impedance. The results were correlated with the differences between the angular variations of the extended- and local-reaction surface reflection coefficients, with reasonable success. The following remarks summarize the findings on the significance of predicting the steady-state sound field in a room with surfaces modeled as of extended or local reaction: 127 1. The difference in modeling a surface as of extended or local reaction is not significant when the surface is a single plate, or a single layer of material (solid or porous) with a rigid backing. 2. The difference in modeling a surface as of extended or local reaction is significant when the surface consists of multi-layers of solid or porous material, and includes a layer of fluid with a large thickness relative to the other layers. 3. For surfaces for which the reflection coefficient is different when obtained with an extended- or local-reaction surface impedance, the difference between making the extended- or local-reaction assumption may be significant, depending on the source and receiver positions. It may be significant when the positions are such that near-grazing incidence occurs in a source-receiver path that includes one reflection. This research was a first attempt to study the behaviour of extended- and local-reaction surfaces on predicting the sound field in rooms. A new room-prediction model was developed in order to accomplish this task. The results obtained by applying this model to various room configurations suggest that the local-reaction assumption of a room surface may be significantly inaccurate in some cases. 6.1 Future work and improvements Improvements to the model are recommended for future work on applying this model to a wider variety of room configurations. Improvements to the model could be made to increase its accuracy, computational efficiency and functionality. The suggested improvements are as follows: 1. Accuracy - Use the exact source-receiver path instead of estimating it as the path from the last image source to the point on a plane that contains the receiver point. This would require a data structure to store all of the reflected surfaces and image sources calculated for the source-receiver path. After the source-receiver path is found, the angles o f incidence and, hence, the reflection coefficients for all surfaces would be determined. 128 2. Accuracy and computational efficiency - Implement an adaptive beam-tracing algorithm16. Currently, the model traces a beam by its centre ray, thus determining a single reflecting surface for the beam. The entire beam is reflected from this surface, even though there is the possibility of a beam intersecting multiple surfaces. The result is that the receiver may not be detected when it should be, or vice-versa. The adaptive beam-tracing algorithm subdivides a beam that intersects multiple surfaces into child beams that adapt to the polygon formed by the intersecting beam and one surface. This would improve both accuracy and computational run-time because fewer beams would be required. The current model required 2000 beams and, combined with the time to compute using the transfer-matrix model, the computation run-time for a single simulation was very slow; in some cases, the run-time was 2 - 3 days on a Pentium 200 MHz computer. 3. Computational efficiency - Pre-calculate the reflection coefficients of multi-layered surfaces at appropriate frequencies and angles, and use an appropriate interpolation scheme in the beam-tracing model. In the current implementation, the transfer-matrix function was called each time the reflection coefficient had to be computed, for the applicable incident angle and frequency-domain vector. This operation increased the computation run-time of the room-prediction model substantially. A suitable interpolation scheme would increase the computational efficiency considerably. Accuracy, however, might be lost. 4. Computational efficiency while maintaining accuracy - As in 3, pre-calculate the properties of a multi-layered surface with the transfer-matrix model, but as a function that takes frequency and incident angle as inputs. This may involve computing the transfer-matrix formulation symbolically prior to beam tracing, and then substituting for frequency and angle within the beam-tracing program. This would improve execution speed without a penalty in accuracy. 5. Functionality - Incorporate diffuse reflections into the model and predict the impulse response of a room. It is not clear how to incorporate diffuse reflection from an extended-reaction surface; this requires further research. By including diffuse reflections, however, the model would be able to predict the impulse response of a 129 room. In addition to predicting the steady-state sound-pressure level, reverberation times and other room-acoustical quantities could be obtained. The above list is not exhaustive, but it does identify the main areas of improvement recommended as future work on this model and to continue this research. The model is unique in its application to rooms with extended-reaction surfaces, but the improvements suggested will increase the practicality of the model. The model could then be applied to study other aspects of sound in rooms. This could involve, for example, comparing predictions made by energy-based models to establish their accuracies and the significance of phase at high frequency. It could also be applied to study sound transmission from one room or structure to another, with accurate modeling of the multi-layered transmitting surfaces; the sound transmission coefficient is predictable by the transfer-matrix model. As a phase model, it could be applied to active noise-control applications in rooms. The application of this model is not to limited rooms, however. It could be applied to the design of noise-control enclosures, to outdoor acoustical environments or to underwater environments. 130 Nomenclature c = wave speed (m/s) e,y = strain tensor elements / = frequency (Hz) h = plate thickness (mm) k = wave number (m"1) / = layer thickness (mm) p = acoustical pressure (Pa) pref= reference pressure, 20 uPa. r - reflection order r0 = source position in Cartesian coordinates (m) ri = receiver position in Cartesian coordinates (m) t = time (s) ui = particle displacement in the xj direction (m) U2 = particle displacement in the x2 direction (m) u3 = particle displacement in the xs direction (m) v/ = particle velocity in the xj direction (m/s) v2 = particle velocity in the x2 direction (m/s) V3 = particle velocity in the X3 direction (m/s) E = Young's modulus (GPa) K = bulk modulus (Pa) L = |R - S| = path length from source S, to receiver R. Lojeam ~ sound-pressure level predicted with the beam-tracing model. L p = sound-pressure level (dB) Lpjmage = sound-pressure level predicted with the image-phase model. iV = shear modulus (Pa) P, Q, and R are Biot elasticity coefficients (Q is a coupling coefficient) Pi = incident pressure at a boundary (Pa) Pr = reflected pressure at a boundary (Pa) P, = transmitted pressure at a boundary (Pa) R = pressure reflection coefficient Reff = Ri><R2x...xRr, the effective reflection coefficient as a multiple of all the reflection coefficients from reflected surfaces, up to the r t h reflection. T = pressure transmission coefficient A = pressure amplitude constant for a spherical wave Bt = component of body force B . Kb = bulk modulus of the frame in a porous material (MPa) K/= bulk modulus of the air in a porous material (MPa) N = shear modulus of the frame (MPa, complex to include viscous losses) Pr = Prandtl number W = sound power (Watts) Wref= reference power, 10"12 Watts. Z = impedance (kg/m -s, or rayls) Zc = characteristic impedance of air = 415 rayls Zpp(ri,ro) = impedance transfer function between source and receiver a = energy absorption coefficient a = tortuosity (for porous materials) (p - porosity 77 = dynamic viscosity (Ns/m ) v= Poisson's ratio p = density (kg/m3) po = density of air (kg/m ) pa = density of the frame in a porous material (kg/m ) p = energy reflection coefficient 9 = angle of incidence ay = stress tensor elements (Pa) CT= flow resistivity (rayls/m) T= energy transmission coefficient co = angular frequency (rad/s) 132 ALp = arithmetic difference in decibel levels A = Viscous dimension (m) A' = Thermal dimension (m) B i b l i o g r a p h y 1 H. Kuttruff. Room Acoustics. Applied Science Publishers, 1976, chapter III. 2 K. R. Fyfe, J. P. Coyette, and P. A. van Vooren. Acoustical and Elasto-Acoustic Analysis Using Finite Element and Boundary Element Methods. Sound and Vibration, December 1991, pp. 16 - 22. 3 V. Easwaran and A. Craggs. On Further Validation and Use of the Finite Element Method to Room Acoustics. Journal of Sound and Vibration, 187(2), 1995, pp. 195-212. 4 J. P. Coyette. Finite Element and Boundary Element Methods for Transient Acoustic Problems. Proceedings of IMAC'95, pp. 1345- 1351. 5 M. Terao and H. Sekine. A Numerical Analysis of Sound Field of a Long Space by a Sub-Region Coupling Approach. Proceedings of Inter-Noise '96, pp. 3007 - 3010. 6 F. Ramiandrasoa, M. A. Hamdi, and M. Haddar. A Sub-Domain Decomposition Method For Large Acoustic Cavities. Proceedings of Inter-Noise'96, pp. 2975 - 2978. 7 J. B. Allen and D. A. Berkley. Image method for efficiently simulating small-room acoustics. Journal of the Acoustical Society of America, 65(4), 1979, pp. 943 -950. 8 J. Borish. Extension of the image model to arbitrary polyhedra. Journal of the Acoustical Society of America, 75(6), 1984, pp. 1827- 1836. 9 A. Kulowski. Algorithmic Representation of the Ray Tracing Technique. Applied Acoustics 18 (1985), pp. 449 -469. 1 0 E. De Geest and H. Patzold. Comparison Between Room Transmission Functions Calculated With a Boundary Element Method and a Ray Tracing Method Including Phase. Proceedings of Inter-Noise '96, pp. 3177 - 3180. 1 1 N. Dadoun, D. G. Kirkpatrick, and J. P. Walsh. The Geometry of Beam Tracing. Association for Computer Machinery, 1985, pp. 55 —61. 1 2 D. van Maercke and J. Martin. The Prediction of Echograms and Impulse Responses within the Epidaure Software. Applied Acoustics, 38(1993), pp. 93 - 114. 13 A. Farina. RAMSETE - A New Pyramid Tracer For Medium and Large Scale Acoustic Problems. Proceedings of Euro-Noise '95, pp. 55 - 60. 14 A. Farina. Pyramid tracing vs. ray tracing for the simulation of sound propagation in large rooms. Computational Acoustics. Computational Mechanics Publications, 1995, pp. 109 - 116. 15 T. Lewers. A Combined Beam Tracing and Radiant Exchange Computer Model of Room Acoustics. Applied Acoustics, 38 (1993), pp. 161 - 178. 1 6 I. A. Drumm and Y. W. Lam. The adaptive beam-tracing algorithm. The Journal of the Acoustical Society of America, 107(3), 2000, pp. 1405 - 1412. 1 7 D. van Maercke. Simulation of Sound Fields in Time and Frequency Domain Using a Geometrical Model. Proceedings of the 12th ICA, Toronto, 1986, paper E l 1-7. 134 1 5 M. Vorlander. Simulation of the transient and steady-state sound propagation in rooms using a new combined ray-tracing/image-source algorithm. Journal of the Acoustical Society of America, 86(1), July 1989, pp. 172 -178. 1 9 B. L. Dalenback. Room acoustic prediction based on a unified treatment of diffuse and specular reflection. Journal of the Acoustical Society ofAmerica, 100(2), August 1996, pp. 899 - 909. 2 0 M. Hodgson. On the prediction of sound fields in large empty rooms. Journal of the Acoustical Society of America, 84(1), 1988, pp. 253 - 261. 2 1 Kinsler. Fundamentals of Acoustics. John Wiley and Sons, 2000. 2 2 D. L. Folds and C. D. Loggins. Transmission and reflection of ultrasonic waves in layered media. Journal of the Acoustical Society of America, 62(5), 1977, pp. 1102 - 1109. 2 3 M. E. Delany and E. N. Bazley. Acoustical properties of fibrous materials. Applied Acoustics, 3 (1970), pp. 105 -16. 2 4 K. Attenborough. Acoustical characteristics of rigid fibrous absorbents and granular media. Journal of the Acoustical Society of America, 73 (1983), pp. 785 - 799. 2 5 D. L. Johnson et. al.. Theory of dynamic permeability and tortuosity in fluid-saturated porous media. Journal of Fluid Mechanics, 176, 1987, pp. 379-402. 2 6 J. F. Allard. Propagation Of Sound In Porous Media. Elsevier Applied Science, 1993, chapter 5. 2 7 M. A. Biot. The theory of propagation of elastic waves in a fluid-saturated porous solid. I Low frequency range. II Higher frequency range. Journal of the Acoustical Society of America, 28 (1956), pp.. 168 - 191. 2 8 R. Panneton and N. Atalla. Numerical prediction of sound transmission through finite multilayer systems with poroelastic materials. Journal of the Acoustical Society of America, 100(1), 1996, pp. 346 - 354. 2 9 J. F. Allard. Propagation Of Sound In Porous Media. Elsevier Applied Science, 1993, chapter 7. 3 0 E. Mariez, S. Sahraoui, and J. F. Allard. Elastic Constants of Polyurethane Foam's Skeleton For Biot Model. Proceedings of Inter-Noise '96, pp. 951 - 954. 3 1 R. L. Willis et. al.. An experimental-numerical technique for evaluating the bulk and shear dynamic moduli of viscoelastic materials. Journal of the Acoustical Society of America, 102(6), 1997, pp. 3549 - 3555. 3 2 M. R. Stinson and Y. Champoux. Propagation of sound and the assignment of shape factors in model porous materials having simple pore geometries. Journal of the Acoustical Society of America, 91(2), 1992, pp. 685 -695. 3 3 Ph. Leclaire et. al.. Determination of the viscous and thermal characteristic lengths of plastic foams by ultrasonic measurements in helium and air. Journal of Applied Physics, 80(4), 1996, pp. 2009 - 2012. 3 4 J. F. Allard. Propagation Of Sound In Porous Media. Elsevier Applied Science, 1993, chapter 6. 3 5 B. Brouard, D. Lafarge, and J. F. Allard. A General Method of Modeling Sound Propagation in Layered Media. Journal of Sound and Vibration, 183(1), 1995, pp. 129- 142. 3 6 K. Attenborough. Acoustical characteristics of porous materials. Physics Reports, 82(3), 1982, pp. 179-227. 3 7 C. Klein and A. Cops. Angle Dependence of the Impedance of a Porous Layer. Acustica, 44, 1980, pp. 258 -264. 135 3 8 J. F. Allard. Propagation of Sound In Porous Media. Elsevier Science Publishers, 1993, Chapter 1. 3 9 J. F. Allard. Propagation of Sound In Porous Media. Elsevier Science Publishers, 1993, p. 146. 4 0 J. F. Allard. Propagation of Sound In Porous Media. Elsevier Science Publishers, 1993, p. 168. 4 1 M. R. Hodgson, A. Wareing, and C. Campbell, "A combined beam-tracing and Biot-transfer-matrix model for predicting sound fields in enclosed spaces", Proceedings of Euro-Noise '98, Munich, Germany, 1023 - 1028 (1998). 4 2 J. F. Allard. Propagation of Sound In Porous Media. Elsevier Applied Science, 1993, pp. 133 - 140. 4 3 B. Brouard et. al.. Measuring and prediction of the surface impedance of a resonant sound absorbing structure. Acta Acustica, 2, 1994, pp. 301 - 306. 4 4 J. Guo and M. Hodgson. Investigation of active noise control in non-diffuse sound fields. Proceedings of Active 99, Fort Lauderdale, FL, USA. Vol. 2, 1999, pp. 621 - 632. 4 5 J. F. Allard. Propagation of Sound in Porous Media. Elsevier Applied Science, 1993, pp 150 - 151. 136 APPENDIX A Interface and Transfer Matrices The transfer-matrix model calculates the surface impedance, reflection coefficient and transmission loss (if applicable), for a composite surface made up of layers of materials. The layers are classified into three categories of materials, fluid, elastic-solid and elastic-porous. Each layer has a finite thickness and extends to infinity in the other two directions. For a given layer shown in Figure A l , the velocities and stresses at a point, a, just on the inside of a boundary are represented as: Fluid layer A l Elastic-solid layer A2 A3 where superscript/denotes fluid material superscript s denotes solid material subscript 1 indicates tangential component (relative to surface) subscript 3 indicates normal component (relative to surface). 137 t xl a b c d Layer 1 Layer 2 Figure A . l : Layer diagram. A transfer matrix, T, relates the velocities and stresses from a point just inside the layer, from the left boundary to the right boundary. For example, in Figure A l , the velocities and stresses at the point a, are related to the point b, as: The boundary conditions at an interface between two adjacent layers are written as a set of homogeneous equations involving transfer matrices. With reference to Figure A l for the interface between layers 1 and 2, the boundary conditions are written as: The transfer matrix for a fluid layer was given in Equation 3.10. The transfer matrices for solid and porous layers are more complicated and given in the next section. The interface matrices I and J depend on the boundary conditions at an interface between two adjacent layers. These boundary conditions depend on the left- and right-layer type combination for the two adjacent layers. The I and J interface-matrices for all combinations are given in Section A.2. V(a) = T-V(b) A4 l],2-V(b) + Jl2-V(c) = 0 A5 138 cu u C H CS X &, JO > cu T3 CO CS cs B >-< C*H CO CU I—I CO =3 o o <D H u II a. H O T—-I m cS C3 cu > CO a CU H o o 1 X f l ci X •j«r o o CO tl I 1 I o o CN I I co O O 1 o o x o o 1 r =3. m <~< ~ 3 --> 1 x 1 CN I n O o cf I ci H CO O o C N I co 5*f CM CO O CJ i 1 CO O . o X o o X =3. -~ q 1 ' c O If o o o o CN CN CO X co O O I cS a o '55 co i) I-I I o o CJ Xi c2 CN <3\ f> cu 1 'E CL) s o 1-4 o & CD on cu -»-» O CQ cu cu cu CU X) CU > CS CU Xi cu CU Xi CU X •o m o II C+H he i d E—1 « &cT CU > CS cu cS c S CU co T3 • »—I CS cu cu CU X) 3 G cu > o cu (3 O CL, S O o "cS CU G O ' c o 3 I >^ The three Biot waves are indicated as 8t and are given as follows: 8\ = co 2{PR-Q2) 2{PR-Q2) [pp22+Rpu-2Qpn-fA] [pp22+Rpn-2QpX2+4l] co 82 = — 3 /V PnPn ~ Pu Pn 2 \ A = (Pp22+Rpxx-2QpX2)2-A{PR-Q2\pxxp22-p22) The following are the ratios of the velocity of the frame to that in the air in the porous material, where i = 1 and 2 represents the compressional waves, and i = 3 represents the shear wave. The velocity ratios indicate the dominant medium in which the wave propagates: p3 coLpX2-Q82 P\2 P22 Dl=(P + QpXk?+k?i)-2Nk? i= 1,2. The variables P, Q, and R are elasticity coefficients in the context of Biot theory: 3 b 0 f R = #Kt 140 The variables pu, pn, and p12 are the inertial coupling terms in the context of Biot theory: Pu =Pl+Pa ~j<4 G(co) CO Pu = 'Pa + J°<P 2 G{co) CO P22 =<PPo+Pa-J(7<ri' Pa = PoAa-1) G(co) co G(co) = 1 + 4ja2r/p0tv Cr2K2(p2 The solid-layer transfer matrix, Ts, is determined similarly to the porous-layer transfer matrix as: T =r'(0)-r'(iy A8 The matrix F ^ ) , in Equation A 8 , is a subset of the rp(x?,) matrix in Equation A 7 with the third and fourth rows, and the third and sixth columns, removed. The Tp(x3) matrix elements removed correspond to the components of fluid pressure and velocity that are absent in the solid layer. The TXxi) matrix is completed with a modification of some of the constants from the Tp(x3) matrix. The velocity ratios, jUj, are set to one, and the wave numbers are defined for the compressional and shear waves in the solid. Several changes are also made to the coefficients D, and Ej. The modifications result in the following for the Vs fa) matrix: r ' ( * 3 ) = cokt cos(kp3x3) -jakp3 sm(kp3x3) -D, cos(^3x3) 2jNktkp3 sin(*p3*3) -jcok,sin(kp3x3) cokp3 cos(kp3x3) jD,sm(kp3x3) jcvks3 sin(*,3*3) cokt cos(^3^3) 2jNk,ks3 sm(ks3x3) - 2Nk,k 3 cos(A: 3^3) N(k*3 - kf)cosO^j^) -coks3 cos( s^3x3) -jo)k,sm(ks3x3) -2M,ks3 cos(ks3x3) -jN(k2s3-k*)sm{ks3x3) A 9 141 The following are the wave numbers for the compressional and shear waves in the solid material, where the subscript p designates a compressional wave, and the subscript s denotes a shear wave. The subscript t indicates the tangential wave that is identical for all layers. The second subscript 3 indicates the X3 component. k, =—sin(0) c kn=— and k, = — ' c. CP The compressional wave speed in an elastic-solid cs = I— The shear wave speed in an elastic-solid V P A = / f e + * , 2 ) + 2 / ^ 3 Ex = Ipk] E2=p{kl-k]) N = p. The shear modulus X and p. are the Lame coefficients: 3 V E A E (l + vX 1 " 2 ^) 2 ( 1 + l /) A.2 Interface Matrices The two sets of interface-matrices, I and J, are unique to the boundary conditions at an interface between adjacent layers. The boundary conditions depend on the combination of left and right layer types that share that interface. The combination of left and right layer types at an 142 interface with corresponding boundary conditions is given in Table A l , for the interface between points b and c in Figure A l . The boundary conditions from Table A l , written in the form of Equation A5, yield the following interface-matrices: Fluid-fluid interface: I/j= 2x2 identiy matrix, and J//= -1/ / Solid-solid interface: h,s = 4x4 identiy matrix, and JSiS = -Porous-porous interface: \,P = 6x6 identiy matrix, and Jp = -1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 </>„ 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 <i>b <f>c Fluid-solid interface: 0 -1 1 0 0 0 and J 0 1 0 0 0 0 1 0 0 0 0 1 Solid-fluid interface: lsJ = J, , ,and f,s-Fluid-porous interface: vf,P * e 0 (We) 0 0 0 and J 0 0 0 0" 0 0 0 0 0 1 l P J 0 0 0 1 0 0 0 0 0 0 1 0 Porous-fluid interface: = J / > p , and 143 Solid-porous interface: Porous-solid interface: 1 0 0 0 0 1 0 0 0 1 0 0 0 0 1 0 0 0 0 1 and J , = 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 1 0 0 0 0 1 0 1 P , S = J.s,P>and 144 Table A.l: Boundary conditions for an interface between adjacent layers in combinations of fluid, solid and porous layer types. Layer-type designations: / = fluid layer type, s = solid layer type, and p = porous layer type. For a porous material, <j>b and <j>c are the porosities at points b and c respectively. Layer-Type Combination Boundary Conditions f-f pib) = p{p) v{ib) = v(ic) s-s 0-3*3(6) = (c) o-;3(6) = o-;3(c) v*(6) = v3!(c) v;ib) = v;ic) p-p a{3 ib) + 0-3*3 ib) = a{3 (c) + o-*3 (c) a;3ib) = a;3ic) <?33(P) °L(P) <l>b </>c <pM$)-vW))=<f>Xvf3ic)-vlic)) v*(6) = v*(c) v[ib) = vlic) fs - pib) = 0-3*3 (c) 0 = o-f3(c) v/(6) = v*(c) f-p v3/(ft) = (l-^e)v,'(c) + ^V3 /(c) - ^ c / ? ( * ) = 0'33(C) o = o-;3(C) - ( l - ^ c ) ^ ) = o-*3(c) s-p v*(6) = v*(c) =v3'(C) v1'W = v1'(c) o-;3(6) = o-;3(c) o-;3(6) = o-3/3(c) + o-3v3(c) 145 APPENDIX B Plate Theory Transfer Matrices B . l Single-plate bounded by air on both sides The reflection coefficient of a single, unbounded plate surrounded by air on both sides, is calculated with the transfer-matrix method for oblique incidence and compared with the reflection coefficient obtained from the well known plate theory. Using plate theory and a matrix approach, the surface impedance is calculated; then the reflection coefficient is determined from the surface impedance via Equation 3.31. Also, for plates, the transmission loss is usually the acoustical property of interest and is shown here to illustrate the acoustical properties of plates. For a single plate loaded by the acoustical pressures p(A) and p(B) on either side, shown in Figure B. 1, the equation of motion for the plate is: D dV j- + p-h-ja>-v! =p(A)-p(B) B.l jco dx,4 where D = Eh3 146 V3(A) 1A P(A) plate air X1 VZplale V3(B) P(B) t X3 Figure B.l: Fluid-loading of a single plate The boundary conditions for the plate/air interface must satisfy continuity of velocity and force balance, and therefore are: vi(A) = v? =v3(B) B.2 and p(A) - p(B) -\j- • + p • h • ja\v> = 0 B.3 Z(A) = v3(A) B.3 Z(B) = P(B) Zc v3(B) cos(8) B.5 where ki = ksm(9), the x/ component of the wave number, same for all layers. The boundary conditions stated in Equations B.2 - B.5 are written in matrix form: 147 c = 0 1 0 Z(A) 1 D \ — -k?+p-h-j6)\ 0 0 0 P(A) 0 -1 v3(A) -1 0 P(B) -1 Z(B)_ = 0 B.6 A solution for Equation B.6 exists if and only if |C| = 0. Therefore, the solution for the surface impedance Z(A) that is required is: Z(A) = - B.7 where Cj = the matrix C with the first row and column removed. C2 = the matrix C with the second row and column removed. The reflection coefficient at A, R(A), is obtained from the surface impedance Z(A) via Equation 3.31. The transmission coefficient - the ratio of transmitted to incident pressure denoted by T(B) - is calculated by re-stating the boundary conditions at the plate/air interface on the transmitted side: p(A)--^^ p(B) = 0 R(A) + l B.8 The boundary condition in Equation B.8 is obtained from the force balance of the plate, p(A) =p(B), and with T(B) = R(A) + 1. The first row of the C matrix in Equation B.6 is replaced with the boundary condition stated in Equation B.8, resulting in the C'matrix: C' = T(B) R(A) + l 0 1 0 0 1 D , 4 kx + p-h- jco ^  0 -1 0 P(A) 0 -1 v3(A) -1 0 P(B) -1 Z(B)_ v3(B) = 0 B.9 148 From Equation B.9, the transmission coefficient T(A) is calculated in a similar manner as the surface impedance Z(A) was calculated from Equation B.6. The determinant of the C matrix must be equal to zero for a non-trivial solution to exist. Therefore, the transmission coefficient is calculated by: T(A) = (R(A) + \) C 3 B.10 where Cy '= the matrix C with the first row and column removed. Cs'— the matrix C with the third row and column removed. The transmission loss of a plate is defined as the reciprocal of the energy or intensity transmission coefficient, r, expressed in decibels (dB). The energy transmission coefficient is related to the pressure transmission coefficient, T, as r = | j j 2 . The transmission loss is expressed as: 7X = -20-log(T) (dB) . B . l l B.2 Double plate with air-gap, bounded by air on both sides The double-plate configuration consists of two drywall plates, separated by a 50 mm air-gap. Both plates are identical with properties defined in Table 4.2, the same material used in the single-plate case. As for the single-plate case, the same acoustical properties were predicted with the transfer-matrix model, and compared to those properties predicted with plate theory. The acoustical properties were predicted with plate theory in a similar manner to the single-plate case via a matrix of the boundary condition equations. This was extended to include both plates and the air-gap. With the aid of Figure B.2, the boundary conditions are written for each side of both plates similar to Equations B.2 - B.5: 149 vi(A) = v 3 l =v3(fl) B.12 v3(b) = v3p2 =v3(fi) B.13 / ? ( ^ ) - /7(a) - <{ — • A:,4 + /?, • /z, • y « ^  = 0 B.14 p{b)-p'E)-• &,4 + p 2 • h2 • jcolvf = 0 B.15 Equations B.4 and B.5 are used again for the boundary conditions for surfaces at points A and B, respectively. The fluid-layer transfer matrix, from Equation 3.10, is used to relate the pressures and velocities at points a and b, in Figure B.2: V3(A) P(A) 1A plate air V3(a) V3fb) plate air V3(B) • • • P(a) P(b) l P(B) V3phle\ y X3 1 VZpbtel Figure B.2: Double-Plate Structure. 150 cos(*3a/r • L) k^-jsm(k3mr.L) ^-^j.sin(k3air.L) K3air _co-pi COSI air KK,r-L) B.16 where ^a(> = &a„.cos(6?), the X3 component of the wave number in the air-gap. Equations B. 12 - B.16 combine to form the matrix equation set of boundary conditions for the double plate as: Z(A) 1 D, —-&,4 +p, •//, jco joy 0 0 0 0 •-^-/sin(* 3£) -cos(k,alrL) 0 1 0 0 -cos(^Z) 0 0 copm -jsm(k3L) 1 jco 0 0 0 0 P(a) 0 -1 v , W P(B) -1 0 MB) -1 Z(B) = 0 B.17 The surface impedance for the double plate is found using Equation B.7, substituting C for Cdouble in Equation B.17. The reflection coefficient is then found via Equation 3.31 with Z = Z(A). The transmission loss is also calculated in a similar manner for the single plate. The top row of the Cdouble matrix in Equation B.17 is replaced with the boundary condition from Equation B.8, in terms of the transmission coefficient T(A). The modified Cd0Ubie matrix becomes C'double, as in the case for C—»C in Equation B.9. The transmission coefficient is calculated from Equation B.10, replacing C with C'double, and the transmission loss is determined from Equation B. l l . The transmission loss, reflection coefficient, and surface impedance were predicted using the transfer-matrix model by modeling each plate as a solid layer, with a fluid layer between. As for the case of the single-plate, the transmission loss predicted by plate theory is displayed, in Figure 4.14, for the same frequencies and angle of incidences. The reflection-coefficient curves, predicted by plate theory for the same frequencies and angles as the single-plate case, are shown 151 in Figure 4.15. In both Figures 4.14 and 4.15, curves are shown only for angles of incidence of 0°, 45°, 75°, and 85°, for clarity. The reflection coefficients predicted with plate theory were compared to transfer-matrix predictions for angles of 0°, 15°, 30, 45°, 60°, 75°, and 85°, shown in Figures 4.16 - 4.22, and used for validation of the transfer matrix because the reflection coefficient is used directly in the beam-tracing model. The surface impedance predicted from the transfer-matrix model is shown in Figure 4.23, included to illustrate the extended-reaction behaviour of the double plate. 152 APPENDIX C Matlab Code C l Matlab m-files for the transfer-matrix model o o %asprop.m % %[R, Z, alpha, t , T, TL] = asprop(f, a, Zc, c, layer_prop, rb) % %This f u n c t i o n c a l c u l a t e s the r e f l e c t i o n c o e f f i c i e n t and surface impedance of %a m u l t i - l a y e r surface that i s f i n i t e i n thickness and i n f i n i t e i n the planar % d i r e c t i o n s . Both a c o u s t i c a l p r o p e r t i e s depend on the frequency and angle of %incidence of an incoming plane wave. The surfaces are thus modelled as %extended r e a c t i o n . Surfaces l a y e r s may be one of three m a t e r i a l types; % f l u i d , e l a s t i c - s o l i d , or p o r o u s - e l a s t i c . The p o r o u s - e l a s t i c l a y e r s are %modelled using B i o t theory. The f o l l o w i n g describes input/output ^parameters: % %INPUTS % % %1) a = angle of incidence, normal = 0 (degrees) %2) f = frequency (Hz) %3) Zc = c h a r a c t e r i s t i c impedance of f l u i d %4) c = speed of sound i n f l u i d (m/s) %5) layer_prop = s t r u c t u r e array with two f i e l d s ; f l u i d and s o l i d . % Each element of the s t r u c t u r e array corresponds to each 153 3 l a y e r i n the surface 5 The f i r s t element represents the f i r s t l a y e r i n the surface 5 - see notes below 56) rb = r i g i d backing, 1 = 'yes', there i s a r i g i d backing ( a f t e r the l a s t 5 l a y e r ) , 0 = 'no', no r i g i d backing, i . e . there i s f l u i d 3 ( a i r ) on the t r a n s m i t t e d side sExtended Reaction Surface Layer P r o p e r t i e s % F l u i d p r o p e r t i e s c o n s i s t of: % 1) thickness (mm) % 2) d e n s i t y (kg/m3) % 3) sound speed (m/s) f o r that f l u i d % A i r s a t u r a t e d p o r e - f l u i d p r o p e r t i e s c o n s i s t of % 1) standard a i r pressure (kPa) % 2) a i r d e n s i t y (kg/m3) % 3) dynamic v i s c o s i t y (Ns/m2) % 4) s p e c i f i c heat r a t i o % 5) P r a n d t l Number % S o l i d p r o p e r t i e s (plates) c o n s i s t of: % 1) t h i c k n e s s (mm) %. 2) d e n s i t y (kg/m3) % 3) Young's modulas (GPa) % 4) Poisson r a t i o % S o l i d p r o p e r t i e s (porous) c o n s i s t of: % 1) thic k n e s s (mm) % 2) frame d e n s i t y (kg/m3) % 3) frame shear modulas (MPa) % 4) frame Poisson r a t i o % 5) flow r e s i s t i v i t y (Ns/m4) % 6) p o r o s i t y % 7) t o r t u o s i t y % 8) viscous dimension (m) % 9) thermal dimension (m) %OUTPUTS 154 %1) R = %2) Z = %3) t = pressure r e f l e c t i o n c o e f f i c i e n t , complex surface impedance ( r a y l , mks u n i t s ) c a l c u l a t i o n time ([hrs, mins, sees]) f u n c t i o n [R, varargout] = asprop(f, a, Zc, c, layer_prop, rb) omega = 2 * p i * f ; %angular frequency (rad/s) % a o f i = a*(pi/180); %convert angles to radians ( d i s a b l e when used with the beam-tracing program) a o f i = a; %keep radian angles (for beam-tracing program) fo r 1 = 1:length(layer_prop) ^construct c e l l array f o r l a y e r p r o p e r t i e s LD{1, 1} = l a y e r _ p r o p ( 1 ) . f l u i d ; LD{1, 2} = l a y e r _ p r o p ( 1 ) . s o l i d ; end t i c f o r i = 1:length(f) f o r j = 1:length(aofi) D = mlyrmat(omega(i), a o f i ( j ) , c, LD, rb) ; [r, z, t c , t l ] = calcprop(D, Zc, a o f i ( j ) , rb) ; Z ( i , j ) = z; %surface impedance R ( i , j ) = r; ^pressure r e f l e c t i o n c o e f f i c i e n t i f nargout > 4 T ( i , j ) = t c ; ^pressure t r a n s m i s s i o n c o e f f i c i e n t T L ( i , j ) = t l ; %transmision l o s s , dB end end end alpha = 1 - abs(R). A2; ^energy absorption c o e f f i c i e n t t = toe; p r o p e r t i e s = {'Z\ 'alpha', ' t ' , " I " , ' TL' } ; nout = max(nargout) - 1; 155 f o r i = l:nout s = ['varargout{', num2str(i), '} = ', p r o p e r t i e s { i } , ' ; ' ] ; e v a l ( s ) ; end Q. O %mlyrmat.m %Function to assemble the i n t e r f a c e and t r a n s f e r matrices f o r a l l l a y e r s . %Inputs are s i n g l e frequency and angle, l a y e r data f u n c t i o n D = mlyrmat(f, a, c, LD, rb) k l = f / c * s i n ( a ) ; Icomponent of wave number p a r a l l e l to surface num_of_layers = size(LD, 1); num_of_inter = num_of_layers + 1; % i n i t i a l i z e D matrix D = [ ] ; f o r b = 1:num_of_inter %determine l a y e r types on each side of i n t e r f a c e switch b case 1 l i t = ' f ' ; r l t = l t y p e ( L D { l , 1}, LD{1, 2}); [I J] = i n t e r f l l t , r l t , LD, 1); T = t m a t r i x ( f , k l , r l t , LD{1, 1}, LD{1, 2}); J = J*T; D = [ I , J] ; case num_of_inter l i t = ltype(LD{b - 1, 1}, LD{b - 1, 2}) ; i f rb r l t = 's ' ; e l s e 156 r l t = ' f' ; end [I J] = i n t e r ( l i t , r l t , LD, b) ; DI = z e r o s ( s i z e ( D , 1), s i z e ( J , 2 ) ) ; D2 = zeros ( s i z e ( I , 1), (size(D, 2) - s i z e ( I , - 2 ) ) ) ; D = [D, DI; D2, I, J] ; otherwise l i t = ltype(LD{b - 1, 1}, LD{b - 1, 2}}; r l t = ltype(LD{b, 1}, LD{b, 2 } ) ; [I J] = i n t e r ( l i t , r l t , LD, b) ; T = t m a t r i x ( f , k l , r l t , LD{b, 1}, LD{b, 2}); J = J*T; DI = z e r o s ( s i z e ( D , 1), s i z e ( J , 2 ) ) ; D2 = z e r o s ( s i z e ( I , 1), (size(D, 2) - s i z e ( I , 2 ) ) ) ; D = [D, DI; D2, I, J] ; end end o o %tmatrix.m % % f u n c t i o n to get the t r a n s f e r matrix f o r a l a y e r % f u n c t i o n T = tmatrix(omega, k l , I t , fp, sp) % % c a l c u l a t e the t r a n s f e r matrix depending on the l a y e r type: % switch I t case ' f T = fluid(omega, k l , f p ) ; case 's' T = solid(omega, k l , sp) ; case 'p' T = porous(omega, k l , fp, sp); end % % inter.m % f u n c t i o n to i n i t i a l i z e i n t e r f a c e matrices f o r a given boundary % l i t stands f o r ' l e f t l a y e r type' 157 % r l t stands f o r ' r i g h t l a y e r type' % l a y e r types are e i t h e r 'f f o r f l u i d , 's' f o r s o l i d , or 'p' f o r porous % int_num i s the i n t e r f a c e number, the number of i n t e r f a c e s i s always the % number % of l a y e r s + 1 f u n c t i o n [ I , J] = i n t e r ( l i t , r l t , LD, b) a l t = s t r c a t ( l l t , r l t ) ; switch a l t case ' f f I = eye(2); J = - I ; case ' f s ' I = [ 0 - 1 ; 1 0; 0 0] ; J = [ 0 1 0 0 ; 0 0 1 0 ; 0 0 0 1 ] ; case 1 f p ' por = LD{b, 2}(6); I = [0, -1; por, 0; (1-por), 0; 0, 0]; J = [0, (1-por), por, 0, 0, 0; 0, 0, 0, 0, 0, 1; 0, 0, 0, 1, 0, 0; 0, 0, 0, 1, 0]; case 'ss' I = eye(4); J = - I ; case ' s f 1 = [ 0 1 0 0 ; 0 0 1 0 ; 0 0 0 1 ] ; J = [0 -1; 1 0; 0 0] ; case 'sp' I = - [1 0 0 0; 0 1 0 0 ; 0 1 0 0 ; 0 0 1 0 ; 0 0 0 1]; J = [ 1 0 0 0 0 0 ; 0 1 0 0 0 0 ; 0 0 1 0 0 0 ; 0 0 0 1 0 1 ; 0 0 0 0 1 0 ] ; case 'pp' p i = LD{b - 1 , 2 } (6) ; p2 = LD{b, 2} (6); I = eye(6); J = - [ 1 0 0 0 0 0 ; 0 1 0 0 0 0 ; 0 (l-p2/pl) p2/pl 0 0 0 ; 0 0 0 1 0 p l / p 2 ) ; 0 0 0 0 1 0 ; 0 0 0 0 0 p l / p 2 ] ; case ' p f por = LD{b - 1, 2}(6); I = [0, (1-por), por, 0, 0, 0; 0, 0, 0, 0, 0, 1; 0, 0, 0, 1, 0, 0; Q, 0, 0, 0, 1, 0]; J = [0, -1; por, 0; (1-por), 0; 0, 0]; case 'ps' 1 = [ 1 0 0 0 0 0 ; 0 1 0 0 0 0 ; 0 0 1 0 0 0 ; 0 0 0 1 0 1 ; 0 0 0 0 1 0 ] ; J = - [ 1 0 0 0; 0 1 0 0 ; 0 1 0 0 ; 0 0 1 0 ; 0 0 0 1]; end a. "5 % f l u i d . m % f u n c t i o n to c a l c u l a t e the t r a n s f e r matrix f o r a f l u i d l a y e r o o f u n c t i o n Tf = f l u i d ( a f r e q , k l , fp) %get m a t e r i a l parameters: h = f p ( l ) * 0 . 0 0 1 ; Iconvert to m d e n s i t y = fp (2) ; c = f p (3) ; % c a l c u l a t e wave number: k = afreq/c; % c a l c u l a t e component of wave number i n the d i r e c t i o n normal to surface: k3 = s q r t ( k A 2 - k l A 2 ) ; % c a l c u l a t e the t r a n s f e r matrix: Tf = zeros(2, 2) ; t f l = a f r e q * d e n s i t y / k 3 ; Tf(1, 1) = cos(k3*h); T f ( l , 2) = t f l * l j * s i n ( k 3 * h ) ; Tf(2, 1) = ( 1 / t f l ) * l j * s i n ( k 3 * h ) ; Tf(2, 2) = cos(k3*h); o, o %solid.m g, o I f u n c t i o n to C a l c u l a t e t r a n s f e r matrix f o r a s o l i d l a y e r 159 f u n c t i o n Ts = s o l i d ( a f r e q , k l , sp) %get user input parameters: h = sp(1)*1.Oe-3; ^convert to m dens i t y = sp(2); E = sp(3)*1.0e9; Iconvert to Pa pr = sp (4 ) ; %C a l c u l a t e Lame c o e f f i c i e n t s : lambda = ( p r * E ) / ( ( l + p r ) * ( l - 2 * p r ) ) ; mu = E/(2* (1 + pr) ); %Set displacement r a t i o vector to empty (for porous l a y e r s o n l y ) : dr = [0 0 0] ; % c a l c u l a t e l o n g i t u d n a l and shear wave speeds (no a t t e n u a t i o n ) : cp = sqrt((lambda + 2*mu)/density); cs = sqrt(mu/density); % C a l c u l a t e wave numbers: kp = afreq/cp; ks = afreq/cs; % C a l c u l a t e component of wave numbers i n the 3 d i r e c t i o n : kp3 = s q r t ( k p A 2 - k l A 2 ) ; ks3 = s q r t ( k s A 2 - k l A 2 ) ; wn = [ k l kp3 ks3 0]; % c a l c u l a t e a d d i t i o n a l parameters: DI = lambda*(kp3 A2 + k l A 2 ) + 2*mu*kp3A2; E l = 2*mu*kl; E2 = mu*(ks3 A2 - k l A 2 ) ; C = [DI E l E2 0 0 0]; %Intermediate matrices (tau) to determine t r a n s f e r matrix (tp) f o r the s o l i d l a y e r 160 t s l = gmffafreq, -h, wn, dr, C, ' s ' ) ; ts2 = gmf(afreq, 0, wn, dr, C, ' s' ) ; % C a l c u l a t e the t r a n s f e r matrix f o r the porous l a y e r : Ts = t s l * p i n v ( t s 2 ) ; % o %porous.m % f u n c t i o n to c a l c u l a t e the t r a n s f e r matrix f o r a porous l a y e r f u n c t i o n Tp = porous(afreq, k l , fp, sp) %get f l u i d parameters f o r the porous l a y e r : % StdPress = f p ( l ) * 1 0 0 0 ; PFDen = fp(2) ; PFVis = f p ( 3 ) ; SpHeatRatio = f p ( 4 ) ; Pr = f p ( 5 ) ; %get s o l i d parameters f o r the porous l a y e r : o. o h = sp(l)*0.001; Iconvert to m frame_den = sp(2); frame_sm = sp(3)*1.0e6; %convert to Pa frame_pr = sp(4); flow_res = sp(5); por = sp(6); t o r t = sp(7); vis_dim = sp(8); therm_dim = sp(9); % C a l c u l a t i o n of intermediate parameters: o % E l a s t i c i t y c o e f f i c i e n t s : % Kf = Bulk modulas of a i r i n s i d e pores % Kb = Bulk modulas of frame % P, Q, R = e l a s t i c i t y c o e f f i c i e n t s r e q u i r e d i n B i o t theory a o 161 % t h e next 4 l i n e s a r e i n t e r m e d i a t e c a l c u l a t i o n s f o r K f K f a = S p H e a t R a t i o * S t d P r e s s ; K f b = a f r e q * P r * t h e r m _ d i m A 2 ; K f c = s q r t ( l + l j * P F D e n * K f b / ( 1 6 * P F V i s ) ) ; K f d = 1 j * t h e r m _ d i m A 2 * P r * a f r e q * P F D e n ; K f = K f a / ( S p H e a t R a t i o - ( S p H e a t R a t i o - 1)*(1 + ( 8 * P F V i s ) / K f d * K f c ) A ( - 1 ) Kb = (2*frame_sm*( frame_pr + 1 ) ) / ( 3 * ( 1 - 2 * f r a m e _ p r ) ) ; P = 4 /3*frame_sm + Kb + (1 - p o r ) A 2 / p o r * K f ; Q = K f * ( 1 - p o r ) ; R = p o r * K f ; % I n e r t i a l c o u p l i n g t e r m s : % r h o _ l l , r h o _ 1 2 , and rho_22 % i n t e r m e d i a t e p a r a m e t e r s g, r , and r h o _ a o o %the next 2 l i n e s a r e i n t e r m e d i a t e c a l c u l a t i o n s f o r g g l = 4 j * t o r t A 2 * P F V i s * P F D e n * a f r e q ; g2 = f l o w _ r e s A 2 * v i s _ d i m A 2 * p o r A 2 ; g = s q r t ( 1 + g l / g 2 ) ; r = 1 j * f l o w _ r e s * p o r A 2 * g / a f r e q ; r h o _ a = P F D e n * p o r * ( t o r t - 1 ) ; r h o _ l l = frame_den + r h o _ a - r ; rho_12 = - r h o _ a + r ; rho_22 = p o r * P F D e n + r h o _ a - r ; %Wave numbers f o r a l l B i o t waves: % k p l and kp2 a r e wave numbers f o r t h e two c o m p r e s s i o n a l waves % ks i s t h e wave number f o r the s h e a r wave %the nex t 4 l i n e s a r e i n t e r m e d i a t e c a l c u l a t i o n s f o r k p l , kp2, and ks . ka = 2*(P*R - Q A 2 ) ; kb = P*rho_22 + R * r h o _ l l - 2*Q*rho_12; kc = r h o _ l l * r h o _ 2 2 - r h o _ 1 2 A 2 ; d e l t a = k b A 2 - 2 * k a * k c ; k p l = a f r e q A 2 / k a * ( k b - s q r t ( d e l t a ) ) ; kp2 = a f r e q A 2 / k a * ( k b + s q r t ( d e l t a ) ) ; ks = afreq A2/frame_sm*(kc/rho_22); %Displacement r a t i o s f o r the frame and a i r : % d i s p _ r a t i o l and d i s p _ r a t i o 2 are f o r the 2 compressional waves % d i s p _ r a t i o 3 i s f o r the shear wave %the next 2 l i n e s are intermediate c a l c u l a t i o n s f o r the displacement % r a t i o s d l = a f r e q A 2 * r h o _ l l ; d2 = afreq A2*rho_12; d i s p _ r a t i o l = (P*kpl - d l ) / ( d 2 - Q*kpl); d i s p _ r a t i o 2 = (P*kp2 - d l ) / ( d 2 - Q*kp2); d i s p _ r a t i o 3 = -rho_12/rho_22; dr = [ d i s p _ r a t i o l d i s p _ r a t i o 2 d i s p _ r a t i o 3 ] ; %Wave number components i n the x3 d i r e c t i o n ( d i r e c t i o n of propagation): % k l 3 and k23 are the x3 components f o r the two compressional wave % numbers % k33 i s the x3 component of the shear wave number kl3 = s q r t ( k p l - k l A 2 ) ; k23 = sqrt(kp2 - k l A 2 ) ; k33 = s q r t ( k s - k l A 2 ) ; wn = [ k l k l 3 k33 k23]; % C o e f f i c i e n t s r e q u i r e d to c a l c u l a t e intermediate matrices f o r the t r a n s f e r %matrix Dl = (P + Q * d i s p _ r a t i o l ) * ( k l A 2 + k l 3 A 2 ) - 2*frame_sm*kl A2; D2 = (P + Q * d i s p _ r a t i o 2 ) * ( k l A 2 + k23 A2)- 2*frame_sm*kl A2; E l = 2*frame_sm*kl; E2 = frame_sm*(k33 A2 - k l A 2 ) ; F l = ( R * d i s p _ r a t i o l + Q ) * ( k l A 2 + k l 3 A 2 ) ; F2 = ( R * d i s p _ r a t i o 2 + Q ) * ( k l A 2 + k23 A2); C = [Dl E l E2 D2 F l F2] ; 163 % Intermediate matrices (tau) to determine t r a n s f e r matrix (tp) f o r the %porous l a y e r t p l = gmf(afreq, -h, wn, dr, C, 'p'); tp2 = gmf(afreq, 0, wn, dr, C, 'p'); %C a l c u l a t e the t r a n s f e r matrix f o r the porous l a y e r : Tp = t p l * p i n v ( t p 2 ) ; o o %gmf.m %Gamma matrix f u n c t i o n to c a l c u l a t e t r a n s f e r matrix %used f o r s o l i d and porous t r a n s f e r matrices o o f u n c t i o n gm = gmf(omega, z, wn, dr, C, I t ) ; % Iomega = angular frequency %z = le n g t h i n the 3 d i r e c t i o n (along thickness of lay e r ) %wn = [ k l kp3 ks3] f o r a s o l i d l a y e r %wn = [ k l k l 3 k33 k23]for a porous l a y e r %dr = [] f o r a s o l i d l a y e r %dr = [ d i s p _ r a t i o l d i s p _ r a t i o 2 disp__ratio3] f or a porous l a y e r %C = [DI E l E2 ] f o r a s o l i d l a y e r %C = [DI E l E2 D2 F l F2] f o r a porous l a y e r % l t = l a y e r type e. o gm = zeros(6, 6); o %elements of gm that are common f o r both s o l i d and porous l a y e r s : gm(l, 1) = omega*wn(1)*cos(wn(2)*z); gm(2, 1) = -lj*omega*wn(2)*sin(wn(2)*z); gm(4, 1) = -C(1)*cos(wn(2)*z); gm(5, 1) = lj*C(2)*wn(2)*sin(wn(2) * z ) ; gm(l, 2) = -lj*omega*wn(l)*sin(wn(2)*z) ; 164 gm(2, 2) = omega*wn(2)*cos(wn(2)*z); gm(4, 2) = 1j * C ( 1 ) * s i n ( w n ( 2 ) * z ) ; gm(5, 2) = -C (2) *wn (2) *cos (wn (2) *z) ; •gm(l, 5) = lj*omega*wn(3)*sin(wn(3)*z); gm(2, 5) = omega*wn(1)*cos(wn(3)*z) ; gm(4, 5) = 1j*C(2)*wn(3)*sin(wn(3)*z); gm(5, 5) = C(3)*cos(wn(3)*z); gm(l, 6) = -omega*wn(3)*cos(wn(3)*z); gm(2, 6) = -1j*omega*wn(1)*sin(wn(3)*z); gm(4, 6) = -C(2)*wn(3)*cos(wn(3)*z); gm(5, 6) = - 1 j * C ( 3 ) * s i n ( w n ( 3 ) * z ) ; i f I t == 'p' % a d d i t i o n a l elements f o r a porous l a y e r : gm (3, 1) = -1j *omega*wn(2)*dr(1)*sin(wn(2) gm (6, 1) = -C(5)*cos(wn(2)*z); gm (3, 2) = omega*dr(1)*wn(2)*cos(wn(2)*z); gm (6, 2) = l j * C ( 5 ) * s i n ( w n ( 2 ) * z ) ; gm(l, 3) = omega*wn(1)*cos(wn(4)*z); gm (2, 3) = -1j *omega*wn(4)*sin(wn(4)*z); gm (3, 3) = - l j *omega*wn(4)*dr(2)*sin(wn(4) gm (4, 3) = -C(4)*cos(wn(4)*z); gm (5, 3) = 1j*C(2)*wn(4)*sin(wn(4)*z); gm (6, 3) = -C(6)*cos(wn(4)*z); gm(l, 4) = - l j *omega*wn(1)*sin(wn(4)*z); gm (2, 4) = omega*wn(4)*cos(wn(4)*z); gm (3, 4) = omega*dr(2)*wn(4)*cos(wn(4)*z); gm(4, 4) = l j * C ( 4 ) * s i n ( w n ( 4 ) * z ) ; gm (5, 4) = -C(2)*wn(4)*cos(wn(4)*z); gm (6, 4) = l j * C ( 6 ) * s i n ( w n ( 4 ) * z ) ; gm(3, 5) = omega*wn(1)*dr(3)*cos(wn(3)*z); 165 gm(3, 6) = -1j*omega*wn(1)*dr(3)*sin(wn(3)*z); e l s e %remove rows and c o l s f o r s o l i d matrix gm(:, 3) = [ ] ; gm (: , 3) = [ ] ; gm(3, :) = [ ] ; gm(5, :) = [ ] ; end 166 C.2 Matlab m-files for the beam-tracing program %inputvar.m l e n t e r a l l input v a r i a b l e s here % g l o b a l s ; % i n i t i a l i z e g l o b a l v a r i a b l e s %1) Common V a r i a b l e s % f r e q = 50:10:4000; f r e q = f r e q ' ; Zc = 415; a i r c = 34 4; a i r _ d e n s i t y = 1.2; %frequency (Hz) % C h a r a c t e r i s t i c Impedance, d e f a u l t value f o r %Sound speed, d e f a u l t value f o r a i r Istandard a i r d e n s i t y omega = 2 * p i * f r e q ; k = omega./c; 52) Room Inputs Lx = 3; Ly = 3; Lz = 3; Coordinates{1, 1} Coordinates{2, 1} Coordinates{3, 1} Coordinates{4 , 1} Coordinates{5, 1} Coordinates{6, 1} [0, 0, 0; 0, Ly, 0; 0,. Ly, Lz; 0, 0, Lz] ; [0, Ly, 0; Lx, Ly, 0; Lx, Ly, Lz; 0, Ly, Lz] ; [Lx, Ly, 0; Lx, 0, 0; Lx, 0, Lz; Lx, Ly, L z ] ; [Lx, 0, 0; 0, 0, 0; 0, 0, Lz; Lx, 0, Lz] ; [0, 0, 0; 0, Ly, 0; Lx, Ly, 0; Lx, 0, 0]; [0, 0, Lz; 0, Ly, Lz; Lx, Ly, Lz; Lx, 0, Lz] ; Source = [0.5, 1.5, 2.0] ; Nb_min = 20; Rec = [2.5, 1.5, 1.8] ; %Source Coordinates %Minimum number of beams %Receiver Coordinates 167 Lw = 80; %Sound Power Lw = repmat(Lw, l e n g t h ( f r e q ) , 1); %sound power of source i n dB, frequency dependent W = (10. A(Lw./10 - 12)); %sound power f o r source i n Watts, f o r each beam Po = s q r t ( ( W * a i r _ d e n s i t y * c ) / ( 2 * p i ) ) ; [ C o e f f i c i e n t s , Normals] = planegeo(Source); Nr = 25; % R e f l e c t i o n Order %Beam v e r t i c e s and centre vectors i c o s a _ f r e q = c e i l ( s q r t ( N b _ m i n / 2 0 ) ) ; Nb = 20*icosa_freq~2; [V, F] = i c o s a ( i c o s a _ f r e q ) ; [V, F] = e l i m v e r t (V, F) ; CV = z e r o s ( s i z e ( F , 1), 3); fo r i = l : s i z e ( F , 1) CV( i , :) = c e n t e r ( V ( F ( i , 1), : ) , V ( F ( i , 2), : ) , V ( F ( i , 3), : ) ) ; end %3) Surface Input Parameters o o extended_reaction = [0, 0, 0, 0, 0, 1] ; LD = ['A', 'A', 'A', 'A', 'A', ' B ' ] ; rb = [0, 0, 0, 0, 0, 1] ; %a) Extended Reaction Surfaces (extended_reaction = 1) o o % F l u i d p r o p e r t i e s c o n s i s t of: % 1) thic k n e s s (mm) % 2) d e n s i t y (kg/m3) % 3) sound speed (m/s) f o r that f l u i d % A i r s a t u r a t e d p o r e - f l u i d p r o p e r t i e s c o n s i s t of : % 1) standard a i r pressure (kPa) 168 % 2) a i r d e n s i t y (kg/m3) % 3) dynamic v i s c o s i t y (Ns/m2) % 4) s p e c i f i c heat r a t i o % 5) P r a n d t l Number % S o l i d p r o p e r t i e s (plates) c o n s i s t of: % 1) thic k n e s s (mm) % 2) d e n s i t y (kg/m3) % 3) Young's modulas (GPa) % 4) Poisson r a t i o % S o l i d p r o p e r t i e s (porous) c o n s i s t of: % 1) thic k n e s s (mm) % 2) frame d e n s i t y (kg/m3) % 3) frame shear modulas (MPa) % 4) frame Poisson r a t i o % 5) flow r e s i s t i v i t y (Ns/m4) % 6) p o r o s i t y % 7) t o r t u o s i t y % 8) viscous dimension (m) % 9) thermal dimension (m) % %Each m u l t i - l a y e r s t r u c t u r e v a r i a b l e c o n s i s t s of two f i e l d s : % 1) f l u i d p r o p e r t i e s % 2) s o l i d p r o p e r t i e s %The vector LD c o n s i s t s of indexes that i n d i c a t e which Layer s t r u c t u r e a p p l i e s to that surface l e x . i f LD(3) = 1, then surface #3 has a m u l t i - l a y e r s t r u c t u r e defined by s t r u c t u r e v a r i a b l e L a y e r l Q. O % p o r e _ f l u i d = [101.3; a i r _ d e n s i t y ; 18.2e-6; 1.4; 0.710]; % f l u i d p r o p e r t i e s f o r a porous m a t e r i a l %rb = [0, 0, 0, 0, 0, 0]; % r i g i d backing, 1 = yes, 0 = no o o % SURFACE TYPE A %b) Local Reaction Surfaces (extended_reaction = 0) %Define SurfaceA: 169 %SurfaceA{l, 1} = [ ] ; %SurfaceA{l, 2} = abs2imp(0.1); %SurfaceA{l, 3} = [ ] ; % %bmtrace.m %beam t r a c i n g main program c l e a r a l l inputvar; % p = z e r o s ( l e n g t h ( f r e q ) , s i z e ( F , 1)); % t i c ; % s t a r t timer f o r b = 1:size(F, 1) u = [V(F(b, 1), : ) ; V(F(b, 2), : ) ; V(F(b, 3), : ) ; CV(b, : ) ] ; S = Source; R_eff = o n e s ( l e n g t h ( f r e q ) , 1); f o r ro = l:Nr r e c e i v e r _ d e t e c t e d = d e t e c t r ( S , Rec, u ( l , : ) , u(2, : ) , u(3, : ) ) ; i f r e c e i v e r _ d e t e c t e d RecPlane = u(4, : ) ; R e c P l a n e ( l , 4) = dot(RecPlane, Rec); %RP i s plane c o e f f i c i e n t s that contain the r e c e i v e r point R r = d i s t ( R e c P l a n e , S, u(4, : ) ) ; % r i s t o t a l path length from l a s t source to R-plane along centre vector z = - l j * k * r ; phase_factor = exp(z); p _ i n s t = ( P o / r ) . * R _ e f f . * ( p h a s e _ f a c t o r ) ; p(:, b) = p ( : , b) + p _ i n s t ; end 170 [S_image, v, r c , su r f ] = p r o j e c t ( S , u); %reset source, d i r e c t i o n v e c t o r s , and c a l c u l a t e e f f e c t i v e r e f l e c t i o n c o e f f i c i e n t S = S_image; u = v; R_eff = R _ e f f . * r c ; end % c a l c u l a t e pressure at r e c e i v e r f o r beam end p b _ t o t a l = sum(p, 2); % t o t a l a c o u s t i c a l pressure (sum pressure f o r each beam a r r i v i n g at r e c e i v e r p o i n t ) . p b _ t o t a l = p b _ t o t a l . * c o n j ( p b _ t o t a l ) ; ^pressure squared, vs. frequency p _ t o t a l = sum(pb_total); %sum of the squared pressures p_ave = mean(pb_total); %average of the squared pressures Calculation_Time = toe; %stop timer % C a l c u l a t i o n time [hrs, mins, sees] = tim e ( C a l c u l a t i o n _ T i m e ) ; t = [hrs, mins, sees]; %output s p l = 1 0 . * l o g l 0 ( p b _ t o t a l . / ( 4 e - 1 0 ) ) ; s p l _ t o t a l = 1 0 . * l o g l 0 ( p _ t o t a l . / ( 4 e - 1 0 ) ) ; spl_ave = 10.*logl0(p_ave./(4e-10)); Q, o % b t _ l o c a l . m % t h i s i s a s p e c i a l m o d i f i c a t i o n to the main beam t r a c i n g program % I t i s used to c a l c u l a t e soun-pressure l e v e l s with m u l t i - l a y e r e d surfaces %modeled as l o c a l - r e a c t i o n , using the t r a n s f e r - m a t r i x model. Q. O %Note that the c l a s s i f i c a t i o n "extended r e a c t i o n " f o r the purpose of t h i s %program i m p l i e s using the 171 % t r a n s f e r matrix model, g l o b a l r c _ l o c a l g l o b a l z _ l o c a l r c _ l o c a l = z e r o s ( l e n g t h ( f r e q ) , 6); z _ l o c a l = z e r o s ( l e n g t h ( f r e q ) , 6); f o r i = 1:6 s u r f a c e _ c l a s s = ['Surface', L D ( i ) ] ; i f e x tended_reaction(i) %NOTE that "extended" here i m p l i e s using t r a n s f e r matrix model e v a l ( [ ' l a y e r _ d a t a = ', s u r f a c e _ c l a s s , ' ; ' ] ) ; [rc l o c a l ( : , i ) , z _ l o c a l ( : , i ) ] = asprop(freq, 0, Zc, c, l a y e r _ d a t a , r b ( i ) ) ; %implement m u l t i - l a y e r surface t r a n s f e r matrix model end end p = z e r o s ( l e n g t h ( f r e q ) , s i z e ( F , 1 ) ) ; t i c ; % s t a r t timer f o r b = l : s i z e ( F , 1) u = [V(F(b, 1), : ) ; V(F(b, 2), : ) ; V(F(b, 3), : ) ; CV(b, : ) ] ; S = Source; R_eff = o n e s ( l e n g t h ( f r e q ) , 1); f o r ro = l:Nr r e c e i v e r _ d e t e c t e d = d e t e c t r ( S , Rec, u ( l , : ) , u(2, : ) , u(3, : ) ) ; i f r e c e i v e r _ d e t e c t e d RecPlane = u(4, : ); (1, 4) = dot(RecPlane, Rec); %RP i s plane c o e f f i c i e n t s that contain the r = d i s t ( R e c P l a n e , S, u(4, r e c e i v e r point R % r i s t o t a l path length from l a s t source to R-plane along centre vector 172 z = - l j * k * r ; phase_factor = exp(z); p _ i n s t = ( P o / r ) . * R _ e f f . * ( p h a s e _ f a c t o r ) ; p(:, b) = p(:, b) + p _ i n s t ; end [S_image, v, r c , surf] = p r o j e c t 2 ( S , u); %reset source, d i r e c t i o n v e c tors, and c a l c u l a t e e f f e c t i v e r e f l e c t i o n c o e f f i c i e n t S = S_image; u = v; R_eff = R _ e f f . * r c ; msg = s p r i n t f ( ' S u r f a c e : L. R.VtBeam No.: %4d / % 4 d \ t R e f l e c t i o n No.: %2d/%2d', b, Nb, ro, Nr) ; disp(msg) end % c a l c u l a t e pressure at r e c e i v e r f o r beam end pb t o t a l = sum(p, 2); % t o t a l a c o u s t i c a l pressure (sum pressure f o r each beam a r r i v i n g at r e c e i v e r p o i n t ) . pb t o t a l = p b _ t o t a l . * c o n j ( p b _ t o t a l ) ; %pressure squared, vs. frequency p t o t a l = sum(pb_total); %sum of the squared pressures p ave = mean(pb_total); %average of the squared pressures Calculation_Time = toe; %stop timer % C a l c u l a t i o n time [hrs, mins, sees] = tim e ( C a l c u l a t i o n _ T i m e ) ; t = [hrs, mins, sees]; %output 173 s p l = 1 0 . * l o g l 0 ( p b _ t o t a l . / ( 4 e - 1 0 ) ) ; s p l _ t o t a l = 1 0 . * l o g l 0 ( p _ t o t a l . / ( 4 e - 1 0 ) ) ; spl_ave = 10.*logl0(p_ave./(4e-10)); Q %icosa.m % f u n c t i o n to construct an icosahedron of frequency, f. There are 20*f A2 v e r t i c e s c a l c u l a t e d . f u n c t i o n [V, F] = i c o s a ( f ) %V = vertex matrix, a 20*f A2x3 matrix of v e r t i c e s that make up the icosahedron %F = face matrix, o p t i o n a l , used to construct patches, deg = pi/180; tao = 0.5*(1+sqrt(5)); %golden r a t i o v = 1; % i n i t i a l i z e vertex index k = 1; %CONSTRUCT THE FIRST FACE IN THE POSITIVE X, Y, Z QUADRANT. %THE OTHER 19 FACES WILL BE CONSTRUCTED FROM THIS ONE USING THE ROTATION FORMULAS fo r i = 0:f fo r j = 0:i xl ( v ) = j ; yl ( v ) = i z l ( v ) = f v = v + 1 ; end end x = x l * s i n ( 7 2 * d e g ) ; y = y l + xl*cos(72*deg); z = repmat(f/2, 1, l e n g t h ( z l ) ) + z l / t a o ; % r o t a t e coordinates to e s t a b l i s h g l o b a l coordinate system: p s i l = 18*deg; t h l = 0; 174 - j ; - i ; ^increment vertex index by 1 phi = 0; D = e u l e r a n g ( p s i l , t h l , p h i ) ; f o r i = l : l e n g t h ( x ) X ( l : 3 , i ) = D * [ x ( i ) ; y ( i ) ; z ( i ) ] ; end % s o r t the coordinates to be c o n s i s t e n t with the F matrix m = 1; h = 1; k = 1; d l = 1; d2 = 1; f o r i = l : f m = m + i ; n = m + i ; X ( l : 3 , m:n) = f l i p l r ( X ( 1 : 3 , m:n)); k = k + i ; d l = [ d l , k ] ; d 2 ( i + 1) = d2(i) + ( i + 1) ; end ^convert to s p h e r i c a l coordinates [ t h e t a l , p h i l , r] = cart2sph(X(1, : ) , X(2, : ) , X(3, : ) ) ; r l = ones(1, l e n g t h ( t h e t a l ) ) ; %Now with r = 1 (unit vectors) f o r a l l v e r t i c e s , convert back to c a r t e s i a n coordinates [ p i , p2, p3] = s p h 2 c a r t ( t h e t a l , p h i l , r l ) ; % C a l c u l a t e F F l = face m a t ( f ) ; %Now construct the other faces v = v - 1; ^remaining top 4 faces: 175 [ t h e t a _ s l , p h i _ s l , F _sl] = c a l c f a c e ( t h e t a l , p h i l , F l , v) %Faces 6 - 1 0 %The 5 t r i a n g l e s around the equator that point downward %Rotation (Euler) angles: p s i = 126*deg; th = acos(-tao/(2 + t a o ) ) ; ph = 54*deg; % c a l c u l a t e r o t a t i o n (coordinate transformation) matrix: D = e u l e r a n g ( p s i , t h , ph) ; f o r i = 1:length(pi) x ( l : 3 , i ) = i n v ( D ) * [ p i ( i ) ; p 2 ( i ) ; p 3 ( i ) ] ; end % f i n d s p h e r i c a l coordinates of new base face [theta2, phi2, r2] = cart2sph(x(1, : ) , x(2, : ) , x(3, :)) F2 = F l + 5*repmat(v, s i z e ( F l ) ) ; [theta_s2, phi_s2, F_s2] = c a l c f a c e ( t h e t a 2 , phi2, F2, v) %Faces 11 - 15 %The 5 t r i a n g l e s around the equator that point downward %Rotation (Euler) angles: p s i = p i / 2 ; th = acos(tao/(2 + t a o ) ) ; ph = -126*deg; % c a l c u l a t e r o t a t i o n (coordinate transformation) matrix: D = e u l e r a n g ( p s i , t h , ph) ; Dinv = i n v ( D ) ; f o r i = 1:length(pi) x ( l : 3 , i ) = D i n v * [ p i ( i ) ; p 2 ( i ) ; p 3 ( i ) ] ; end % f i n d s p h e r i c a l coordinates of new base face [theta3, phi3, r3] = cart2sph(x(1, : ) , x(2, : ) , x(3, :)) F3 = F l + 10*repmat(v, s i z e ( F l ) ) ; [theta_s3, phi_s3, F_s3] = c a l c f a c e ( t h e t a 3 , phi3, F3, v) ; %Faces 16 - 20 % c a l c u l a t e coordinates f o r base t r i a n g l e f o r the bottom 5: theta4 = t h e t a l + repmat(36*deg, s i z e ( t h e t a l ) ) ; phi4 = - p h i l ; F4 = F l + 15*repmat(v, s i z e ( F l ) ) ; % c a l c u l a t e other 4 bottom t r i a n g l e s : [theta_s4, phi_s4, F_s4] = c a l c f a c e ( t h e t a 4 , phi4, F4, v) ; %Now assemble a l l coordinates: theta = [ t h e t a _ s l , theta_s2, theta__s3, t h e t a _ s 4 ] ; phi = [ p h i _ s l , phi_s2, phi_s3, p h i _ s 4 ] ; r = o n e s ( l , l e n g t h ( t h e t a ) ) ; [x, y, z] = sph2cart(theta, p h i , r ) ; V = [x', y', z ' ] ; F = [ F _ s l ; F_s2; F_s3; F_s4]; % %elimvert.m % f u n c t i o n to e l i m i n a t e redundant v e r t i c e s c a l c u l a t e d f o r the source, f u n c t i o n [v2, f2] = e l i m v e r t ( v , f) %set the to l e r a n c e f o r the d i f f e r e n c e between magnitudes of a p a i r of v e r t i c e s to d i s t i n g u i s h i f they are %equal or unique t o l = l e - 6 ; v2 = v; f2 = f; vrows = s i z e ( v , 1); i = 1; c = 1; redundant_vertices = [ ] ; f r e q = s q r t ( s i z e ( f , l ) / 2 0 ) ; fv = vrows/20; v e r t i c e s = 12 + 30*(freq - 1) + 20*(fv - 3*freq); while c <= v e r t i c e s % f i n d the magnitude of the d i f f e r e n c e between each vertex and a l l others v l = repmat(v(i, 1:3), vrows, 1); v d i f f = v l - v; v d i f f = s q r t ( v d i f f ( : , 1). A2 + v d i f f ( : , 2). A2 + v d i f f ( : , 3 ) . A 2 ) ; %set the magnitude of the d i f f e r e n c e between the vertex being t e s t e d and i t s e l f to be 1, % t h i s i s done to avoid t h i s vertex being deleted. % f i n d indeces of redundant v e r t i c e s j = f i n d ( v d i f f < t o l ) ; i f ~isempty(j) % e l i m i n a t e redundant v e r t i c e s redundant_vertices = [redundant_vertices; j ( 2 : l e n g t h ( j ) ) ] ; redundant_vertices = unique(redundant v e r t i c e s ) ; %replace redundant v e r t i c e s with index of unique vertex i n faces matrix f o r k = l : l e n g t h ( j ) [m, n] = f i n d ( f == j ( k ) ) ; f o r p = 1:length(m) f2(m(p), n(p)) = c; end end end %reset number of rows to new s i z e of the v e r t i c e s matrix, increment index i n t o t h i s matrix i = i + 1; while 1 i f - i s e m p t y ( f i n d ( r e d u n d a n t _ v e r t i c e s == i ) ) i = i + 1; e l s e c = c + 1 ; break; end end end 178 v2(redundant_vertices, : ) = [ ] ; [ i , j ] = find(abs(v2) < t o l ) ; fo r k = 1:length(i) v 2 ( i ( k ) , j ( k ) ) = 0; end o Icenter.m % f i n d the center vector ( i n c a r t e s i a n coordinates) f o r a t r i a n g u l a r face. %the coordinates of the three v e r t i c e s f o r the t r i a n g l e are given. %The three v e r t i c e s f o r the t r i a n g u l a r beam face are u n i t vectors that a l l i e on the u n i t sphere %that defines the source. %These p o i n t s l i e on a plane that i s p a r a l l e l to one that i s tangent to the un i t sphere at a point %along the center vector f o r the t r i a n g u l a r face of the beam. This vector i s a u n i t vector that i s %the u n i t normal to the plane defined by the three v e r t i c e s , f u n c t i o n cv = c e n t e r ( v l , v2, v3) le v = center vector % v l , v2, v3 = v e r t i c e s of t r i a n g u l a r face PR = v2 - v l ; PQ = v3 - v l ; N = cross(PR, PQ); N_length = N(1)"2 + N(2) A2 + N(3) A2; N_length = s q r t ( N _ l e n g t h ) ; n = N/N_length; i f dot(n, v l ) < 0 n = -n; end cv = n; % 179 %proj ect.m % t h i s m - f i l e performs a m a j o r i t y of the t r a j e c t o r y - l o o p operations % p r o j e c t s the beam v e r t i c e s on the i n t e r s e c t i n g surface, and c a l c u l a t e s the r e f l e c t i n g d i r e c t i o n v e c t o r s , %and image source, and r e f l e c t i o n c o e f f i c i e n t . f u n c t i o n [ I , v, r c , r s u r f ] = p r o j e c t ( S , u) g l o b a l Coordinates g l o b a l C o e f f i c i e n t s g l o b a l Normals g l o b a l R _ C o e f f i c i e n t s g l o b a l f r e q g l o b a l Zc g l o b a l c g l o b a l LD g l o b a l rb g l o b a l SurfaceA SurfaceB SurfaceC SurfaceD SurfaceE SurfaceF g l o b a l extended_reaction v = zeros (4, 3); r s u r f = i n t s r f ( u ( 4 , : ) , S); [v(4, : ) , a o f i ] = spec_ref(u(4, : ) , r s u r f ) ; f o r i = 1:3 v ( i , :) = s p e c _ r e f ( u ( i , : ) , r s u r f ) ; end I = i s o u r c e ( S , r s u r f ) ; s u r f a c e _ c l a s s = ['Surface', L D ( r s u r f ) ] ; i f e x t e n d e d _ r e a c t i o n ( r s u r f ) e v a l ( [' l a y e r _ d a t a = ', surface_cla.ss, ' ; ' ] ) ; r c = r e f l e c t ( f r e q , a o f i , Zc, c, l a y e r _ d a t a , r b ( r s u r f ) ) ; %implement m u l t i - l a y e r surface t r a n s f e r matrix model el s e %use p r e - a l l o c a t e d t a b l e of r e f l e c t i o n c o e f f i c i e n t s i f e v a l (• [ ' isempty ( ' , sur f ace _ c l a s s , '{1, 2})']) 180 e v a l ( [ 1 r c = ', s u r f a c e _ c l a s s , '{1, 3};']) e l s e i f e v a l ( [ 1 i s e m p t y ( ' , s u r f a c e _ c l a s s , '{1/ 3})']) e v a l ( [ ' Z a = ', s u r f a c e _ c l a s s , '{1, 2};']) rc = l o c r e a c t ( Z a , Zc, a o f i ) ; % R e f l e c t i o n c o e f f i c i e n t s , a l l frequencies end end % %eulerang.m % r o t a t e e u l e r angles p s i , theta, and phi f u n c t i o n D = e u l e r a n g ( p s i , t h e t a , phi) c l = c o s ( p s i ) ; c2 = c o s ( t h e t a ) ; c3 = c o s ( p h i ) ; s i = s i n ( p s i ) ; s2 = s i n ( t h e t a ) ; s3 = s i n ( p h i ) ; D = [ c l * c 3 - s l * c 2 * s 3 , sl*c3+cl*c2*s3, s2*s3; . .. - c l * s 3 - s l * c 2 * c 3 , -sl*s3+cl*c2*c3, s2*c3;... s l * s 2 , - c l * s 2 , c 2 ] ; o o %planegeo.m % C a l c u l a t e s normal vectors f o r planes with c o e f f i c i e n t s f u n c t i o n [C, n] = planegeo(S); g l o b a l Coordinates no_of_surfaces = s i z e ( C o o r d i n a t e s , 1); C = zeros(no_of_surfaces, 4); n = zeros(no_of_surfaces, 3); f o r i = 1:no_of_surfaces P = C o o r d i n a t e s { i , 1}(1, : ) ; 181 Q = C o o r d i n a t e s { i , 1}(2, : ) ; R = C o o r d i n a t e s { i , 1}(3, : ) ; PR = R - P; PQ = Q - P; N = cross(PR, PQ); N_length = N(1) A2 + N(2) A2 + N(3) A2; N_length = s q r t ( N _ l e n g t h ) ; C ( i , 1:3) = N; C ( i , 4) = dot(N, P); n ( i , 1:3) = N/N_length; chkpt = dot(N, S); i f chkpt>C(i, 4) n ( i , 1:3) = - n ( i , 1:3); end end o X) — %contain.m % f u n c t i o n to determine i f the point P, the i n t e r s e c t i o n of a ray and a plane, l i e s w i t h i n the surface %defined by the coordinates f o r that plane f u n c t i o n H = contain(plane, P) g l o b a l Coordinates g l o b a l Normals l e x t r a c t coordinates and normal vector f o r surface # 'plane' coords = Coordinates{plane, 1}; n = Normals(plane, : ) ; %determine e u l e r angles to r o t a t e new coordinate a x i s 182 [azimuth, elev, r] = cart2sph(n(1), n(2), n ( 3 ) ) ; p s i = azimuth - p i / 2 ; theta = elev - p i / 2 ; phi = 0; % c a l c u l a t e transformation matrix and new coordinates E = e u l e r a n g ( p s i , t h e t a , p h i ) ; coords = coords*E'; P = P*E'; %determine i f the point P l i e s i n s i d e the surface defined by the coordinates H = inpolygon(P(1), P(2), coords(:, 1), coords(:, 2 ) ) ; % %Dist.m %Function to c a l c u l a t e the distance from s t a r t point S along d i r e c t i o n vector u, %to the i n t e r s e c t i o n of the ray and a plane with c o e f f i c i e n t s i n the c o e f f i c i e n t matrix. f u n c t i o n s = d i s t ( p l a n e , S, u) g l o b a l C o e f f i c i e n t s i f length(plane) == 1 PI = C o e f f i c i e n t s ( p l a n e , 1:3); P2 = C o e f f i c i e n t s ( p l a n e , 4); el s e PI = p l a n e d , 1:3) ; P2 = plane(1, 4); end s i = d o t ( P I , S); s2 = dot(PI, u); s l = P2 - s i ; i f s2 == 0 183 s = 0; %A "0" i s returned i f the d i r e c t i o n vector i s p e r p i n d i c u l a r to the normal vector f o r the plane e l s e s = s l / s 2 ; end % % i n t s r f . m %determines the surface that i s i n t e r s e c t e d by a ray f u n c t i o n s r f = i n t s r f ( u , S) g l o b a l C o e f f i c i e n t s g l o b a l Normals P = v a l i d s r f ( u , S); surfaces = P(:, 1); dst = P(:, 2); [min_dst, j ] = min(dst); s r f = s u r f a c e s ( j ) ; % lisource.m % f u n c t i o n to c a l c u l a t e the image source of a point r e f l e c t e d about a plane f u n c t i o n S i = i s o u r c e ( S , plane) g l o b a l C o e f f i c i e n t s g l o b a l Normals n = Normals(plane, : ) ; s = d i s t ( p l a n e , S, n); S i = S + 2*s*n; % % v a l i d s r f . m 184 % C a l c u l a t e s p o s s i b l e surfaces that are i n t e r s e c t e d by a ray. f u n c t i o n P = v a l i d s r f ( u , S) g l o b a l Coordinates g l o b a l C o e f f i c i e n t s g l o b a l Normals j = i ; f o r i = 1:size(Normals, 1) p o s s i b l e _ s u r f a c e = dot(u, Normals(i, : ) ) ; i f p o s s i b l e _ s u r f a c e > 0 d = d i s t ( i , S, u); % c a l c u l a t e distance from S to plane, along u Q = S + d*u; I p o i n t of i n t e r s e c t i o n of ray and plane h i t = c o n t a i n ( i , Q) ; %check i f i n t e r s e c t i o n point l i e s w i t h i n coordinates that define surface boundary i f h i t P ( j , 1) = i ; %record plane number P ( j , 2) = d; %record distance j = j + l ; % i ncrement p o s s i b l e surface index end end end 185 

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0089679/manifest

Comment

Related Items