Towards Impedance Characterization of Carbon-Carbon Ultrasonically Absorptive Coatings via the Inverse Helmholtz Problem Danish Patel , Prateek Gupta , and Carlo Scalo ∗ ∗ † 7 1 Purdue University, West Lafayette, Indiana 47907, USA 0 2 Thomas Rothermel, and Markus Kuhn ‡ ‡ n German Aerospace Center (DLR), Stuttgart, Germany a J 9 We present a numerical method to determine the complex acoustic impedance at the 2 open surface of an arbitrarily shaped cavity, associated to an impinging planar acoustic wave with a given wavenumber vector and frequency. We have achieved this by developing ] n thefirstinverseHelmholtzSolver(iHS),whichimplicitlyreconstructsthecomplexacoustic y waveform–at a given frequency–up to the unknown impedance boundary, hereby providing d the spatial distribution of impedance as a result of the calculation for that given frequency. - We show that the algebraic closure conditions required by the inverse Helmholtz problem u l are physically related to the assignment of the spatial distribution of the pressure phase f over the unknown impedance boundary. The iHS is embarassingly parallelizable over the . s frequency domain allowing for the straightforward determination of the full broadband c impedance at every point of the target boundary. In this paper, we restrict our analysis i s to two-dimensions only. We first validate our results against Rott’s quasi one-dimensional y thermoacoustic theory for viscid and inviscid constant-area rectangular ducts, test our iHS h in a fully unstructured fashion with a geometrically complex cavity, and finally, present a p [ simplified, two-dimensional analysis of a sample of carbon-carbon ultrasonically absorptive coatings (C/C UACs) manufactured in DLR-Stuttgart, and used in the hypersonic transi- 1 tion delay experiments by Wagner et al. in AIAA 2012-5865. The final goal is to model v C/C UACs with multi-oscillator Time Domain Impedance Boundary Conditions (TDIBC) 1 developed by Lin et al. in JFM (2016) to be applied in direct numerical simulations (DNS) 9 focused on the overlying boundary layer, eliminating the need to simultaneously resolve 3 the microscopic porous structure of the C/C UACs. 8 0 . 1 I. Introduction 0 7 Boundary layers undergo laminar-to-turbulent transition as a result of the development and growth of 1 perturbations due to a variety of factors including, but not limited to, free-stream turbulence, free-stream : v Mach number, surface-to-free-stream temperature ratio, surface roughness, surface curvature and sound i radiation.1,2 In the case of hypersonic boundary layers with small free-stream disturbances evolving over X smooth surfaces and slender geometries, Mack modes govern the transition dynamics.3,4 Among them, r a second-modewavesareofparticularimportance: theyaremanifestasconvectivelyamplifiedacousticmodes resonating inside the boundary layer with frequencies typically in the ultrasonic range.5 As a result, several passivecontroltechniquesforhypersonicboundarylayertransitiondelayhavebeenfocusingonthedesignof ultrasonically absorptive coatings, or UAC(figure 1).6 They have been shown to attenuate resonant modes7 viathermo-viscousattenuationofwavestrappedinthecoatings’cavities,and,therefore,todelaytransition. ∗GraduateResearchAssistant,SchoolofMechanicalEngineering †AssistantProfessor,SchoolofMechanicalEngineering,[email protected] ‡ResearchScientist,Dep. SpaceSystemIntegration AssubmittedtoAIAASciTech2017,availableathttp://arc.aiaa.org/doi/pdf/10.2514/6.2017-0460 1of14 AmericanInstituteofAeronauticsandAstronautics Acoustic energy absorption by UACs can be conveniently modeled by assigning equivalent impedance boundary conditions. In general, the (normal) acoustic impedance at a surface, , is defined as, S pˆ(x;ω) Z(x;ω)= x Uˆ(x;ω).nˆ(x) ∀ ∈S where Uˆ and pˆ are Fourier transforms of the velocity and pressure at the surface, respectively, nˆ is the surface normal directed away from the flow side and into the cavity, x is the position vector on the surface, and ω is the angular frequency. Hereafter, vectorial quantities and/or arrays will be indicated in bold. For a quiescent fluid, the acoustic power transmitted through a surface—due to interaction with the constituent frequencies contained in an incident wave packet—with given impedance is, W˙ = ∞ 1Re[Z(ω)] Uˆ (ω) nˆ 2dA dω. (1) ac 2 · Z−∞(cid:26)ZS (cid:12) (cid:12) (cid:27) (cid:12) (cid:12) The acoustic impedance at a surface retains all the neces(cid:12)sary infor(cid:12)mation regarding the two-way acoustic couplingbetweenwavesoneithersideofthesurface. Assuch,iftheacousticimpedanceweretobeimposed, it would control the direction and magnitude of acoustic energy transfer. Figure1: Illustrationofa7o half-angleround-tipconemodelequippedwithC/CUACemployedbyWagner4 to damp second mode instabilities. The microporous structure of such coatings is highly irregular and only statistically predictable (see inset figure). Inthiswork,weproposeastrategytomodelopencavitiessuchastheonespresentinUACsbycalculating theirequivalentacousticimpedance(figure2). Thegoalisthentofitappropriatemulti-oscillatorTDIBCsto thecalculatedimpedancefollowingtheapproachbyLinetal. (2016),8 andapplytheminfullycompressible Navier-Stokes solvers focused only on resolving the overlying transitional boundary layer, hence saving significant computational resources at every time-step compared to pore-resolved cases. (a) (b) Figure 2: Illustration of a plane wave impinging on a porous surface (a) with closed cavities, which can be modeled by assigning equivalent impedance boundary conditions (IBCs) (b). We note that this cannot be achieved with a classic acoustic analysis based on Helmholtz solvers. Such solvers rely on the formulation of an eigenvalue problem which yields a set of eigenmodes and corresponding complex eigenvalues containing frequency and growth/decay rate information. These solvers require con- ditions, and especially impedance, to be specified at all boundaries (figure 3a). We rather need a tool to analyze a cavity with a single open boundary where impedance is to be determined (figure 3b) and not specified or known a priori. 2of14 AmericanInstituteofAeronauticsandAstronautics (a) (b) Figure 3: Sample two-dimensional computational setup for (a) direct and (b) inverse Helmholtz solvers. The direct solver requires impedance at all boundaries to be specified, whereas the inverse solver is able to provide the spatial distribution of the impedance at an open boundary as an outcome of the calculation. In the current study, we achieve this by constructing the first inverse Helmholtz Solver (iHS). We begin with the linearized governing equations for a compressible fluid, coupled with a generalized equation of state, which are then cast into the inverse Helmholtz formulation in the frequency domain where part of the governing equations are extended to the unknown impedance boundary (section II). This problem is, however, rank deficient,9 i.e. unclosed. This issue is resolved by applying appropriate—and physically realizable—algebraic closure conditions which are also discussed below. We then show validation of the iHS based on the acoustic impedance prediction for a constant-area rectangular duct against Rott’s linear thermoacousticwavetheoryforbothviscidandinviscidtestcases(sectionIII).Wethenpresenttheanalysis of an idealized complex toy geometry (section IV), and finally, perform a simplified two-dimensional study of the acoustic absorption characteristics of a real C/C UAC geometry provided by DLR-Stuttgart applying the proposed inverse Helmholtz Solver methodology to a single cavity (section V). II. Methodology: Inverse Helmholtz Solver (iHS) In this section, we outline the inverse Helmholtz problem formulation starting with the linearization of the compressible flow governing equations for a fluid with a general equation of state. Although a more comprehensive derivation with mean flow and gradients is possible, we restrict the scope of this paper to a uniform, quiescent base flow and two-dimensional geometries without loss of generality for the iHS methodology. II.A. Governing Equations The fully non-linear governing equations for a general compressible fluid are, ∂ ∂ (ρ)+ (ρu )=0, (2a) j ∂t ∂x j ∂ ∂ ∂ ∂ (ρu )+ (ρu u )= p+ τ , (2b) i i j ij ∂t ∂x −∂x ∂x j i j ∂ ∂ ∂ (ρe)+ [u (ρe+p)]= (u τ q ), (2c) j i ij j ∂t ∂x ∂x − j j where x , x , and x (or x, y and z) are Cartesian coordinates, u are the velocity components in each of 1 2 3 i these directions, and p, ρ, and e are, respectively, the pressure, density, and internal energy per unit mass. 3of14 AmericanInstituteofAeronauticsandAstronautics The shear stress τ , and heat flux q , are defined as, ij j ∂u ∂u j k τ =2µ +λ δ , ij ij ∂x ∂x i k ∂T q = κ , j − ∂x j where µ and λ are the first and second viscosity coefficients respectively and κ is the coefficient of thermal conductivity. The field variables are decomposed as, p(x;t)=p +p(x;t) 0 ′ T(x;t)=T +T (x;t) 0 ′ ρ(x;t)=ρ0+ρ′(x;t) (3) e(x;t)=e +e(x;t) 0 ′ u (x;t)=u (x;t) i ′i where the subscript 0 denotes the base (quiescent) state and the primes denote fluctuations. Upon substi- tution of the equations (3) into the system of equations (2) and dropping terms higher than first order in primed variables, we obtain the linearized equations, ∂ρ ∂u ′ +ρ ′k =0, (4a) 0 ∂t ∂x k ∂e ∂ρ ∂u ∂u ∂ ∂ ρ ′ +e ′ +ρ e ′k +p ′k κ T =0, (4b) 0 0 0 0 0 ′ ∂t ∂t ∂x ∂x − ∂x ∂x k k k (cid:18) k (cid:19) ∂u ∂p ∂ ρ ′i + ′ = τ . (4c) 0 ∂t ∂x ∂x i′j i j Inthegeneraltheoryforlinearacoustics,primarythermoviscouslossesarecapturedviatermscontaining the dynamic viscosity, µ in τ and the thermal conductivity, κ. In addition to these losses, in the ultrasonic ij regime, bulk viscous losses become significant and are accounted for in λ, the coefficient of second viscosity, which can be calculated from the bulk viscosity, µ λ+(2/3)µ. For the latter, we refer to the findings by b ≡ Lin et al.10 where the bulk viscosity, and hence the second viscosity, λ, can be expressed as a function of the frequency, base pressure and temperature. The total differentials of pressure and internal energy for a general equation of state read, 1 α v dp= dρ+ dT (5) ρβ β T T de=c dT B dρ (6) v 0 − ρ∂p 1 with, α = ∂T v , β = v∂p − , B = C0 and, C = T ∂p0 p , v −(cid:16) ∂p (cid:12)(cid:12) (cid:17) T −(cid:18) ∂v(cid:12)T(cid:19) 0 (ρ0)2 0 (cid:18) 0 ∂T0(cid:12)v− 0(cid:19) ∂v T(cid:12) (cid:12) (cid:12) (cid:16) (cid:12) (cid:17) (cid:12) (cid:12) where the coefficient of th(cid:12)ermal expansion, αv,(cid:12)the isothermal compressibility, βT, and cons(cid:12)tants B0 and (cid:12) C can be obtained analytically from any given equation of state. For a calorically perfect ideal gas, these 0 variables revert to, 1 1 α = β = B =C =0. v T 0 0 T p 0 0 Substituting (5) and (6) in equation (4) and eliminating density perturbations yields, ∂p ∂T ∂u ρ β ′ ρ α ′ +ρ ′k =0, (7a) 0 T 0 v 0 ∂t − ∂t ∂x k ∂T ∂p ∂u ∂ ∂ (ρ c +C α ) ′ C β ′ +p ′k κ T =0, (7b) 0 v 0 v 0 T 0 ′ ∂t − ∂t ∂x − ∂x ∂x k k (cid:18) k (cid:19) ∂u ∂p ∂ ∂ ∂ ∂u ρ ′i + ′ =µ u +(µ+λ) ′k , (7c) 0 ∂t ∂x ∂x ∂x ′i ∂x ∂x i k (cid:18) k (cid:19) i (cid:18) k(cid:19) which is the final form that will be discretized. 4of14 AmericanInstituteofAeronauticsandAstronautics II.B. Numerics and Algebraic Formulation Convertingthesystemofequations(7)fromtimedomaintofrequencydomain,collectingalltermscontaining the frequency, ω, and applying a spatial discretization yields the following system of linear equations, iωX +A X =b (8) 0 eig 0 where, T X0 = pˆ Tˆ uˆ1 uˆ2 uˆ3 and, b=0. (cid:16) (cid:17) We adopt a second order discretization wise staggered variable arrangement on an unstructured grid (shown in figure 4), where thermodynamic field variables (pressure p, specific internal energy e, and hence temperature T and density ρ) are defined at the cell centers and velocity field components u at the face i centers. (a) (b) Figure4: Staggeredvariablearrangementovera(a)QUAD,and(b)TRIcell. Allthermodynamicvariables, p,T,ρ,eandc ,arecollocatedoncellcenters( ),whilevelocitycomponents,u ,arelocatedonfacecenters v i • Equation(8)withappropriateboundaryconditionswouldyieldaconventionalHelmholtzproblem,which is essentially an eigenvalue problem that is typically considered in linear acoustic analysis.11 Solving this problem returns a finite set of frequencies (eigenvalues) and wave forms (eigenvectors). Note that homoge- neous boundary conditions (b=0) are required by an eigenvalue solver. On the other hand, in the inverse Helmholtz solver approach, the vector b in equation (8) is no longer null as it contains closure conditions (explained below) to be applied at the unknown impedance boundary (termed IBC). The value of ω is now real and a given, and the system (9) is solved for an extended set of unknowns, which include velocities and pressures at the boundary (uˆ , uˆ , uˆ and pˆ ), 1,IBC 2,IBC 3,IBC IBC in addition to the previous unknowns in X . The momentum equation is now written at the IBC and is 0 symbolically represented by blocks I and II in the overall inverse Helmholtz problem formulation, which reads, b pˆ pˆ p  Tˆ   Aeig II  Tˆ   b..T  uˆ uˆ . iω uuuˆˆˆpˆ321uuˆˆI,,,BIII123BBBCCCC + IIII  uuuˆˆˆpˆ321uuˆˆI,,,BIII123BBBCCCC = bIB............ C . (9)        5of14 AmericanInstituteofAeronauticsandAstronautics In order to close the system, we introduce a pressure variable, pˆ at the IBC. This later allows us IBC to evaluate directly the value of impedance at the IBC, without having to rely on extrapolation techniques that may be inaccurate for complex geometries. Closure conditions are formulated as, pˆ IBC,m+1 =Ψ , (10) m pˆ IBC,m symbolically represented by block III of the system (9), and can be physically related to the shape of the impinging wave. For example, in the case of a two-dimensional planar wave of the shape, p = wave p exp( ik(xcos(θ )+ysin(θ )))—where, p is the amplitude, k the wave number and, θ the angle of n n n ∞ − ∞ incidence of the incident wave on the unknown impedance boundary—equation (10) becomes, pˆ IBC,m+1 =exp i kcos(θ ) x x +ksin(θ ) y y =Ψ (11) n IBC,m+1 IBC,m n IBC,m+1 IBC,m m pˆ − − − IBC,m (cid:16) (cid:8) (cid:0) (cid:1) (cid:0) (cid:1)(cid:9)(cid:17) Finally, to achieve full rank in the system we assign an arbitrary reference pressure, pˆ in an arbitrarily ref chosen cell corresponding to one row of block III. We note that the final value of impedance is independent fromthespecificvalueofthereferencepressuresincetheformerisaratiooftwoquantities,bothproportional to pˆ . ref These two closure conditions are finally incorporated into the system of equations (9), restoring its full rank, allowing to solve for the discrete pressure and velocity components at the IBC as well as inside the domain. Once we have obtained the complete solution, the impedance can directly be calculated as, pˆ IBC Z(x;ω)= (12) ~u .nˆ IBC IBC where, ~u = (uˆ ,uˆ ,uˆ ) is the velocity vector and nˆ the surface normal facing inside, IBC 1,IBC 2,IBC 3,IBC IBC towards the domain. III. Validation Against Rott’s Thermoacoustic Theory The basic working principles of the inverse Helmholtz Solver (iHS), introduced in a general form in the previous section, are here explained with a simple 1-D numerical example. Using a staggered grid approach with a set frequency, ω, Rott’s wave equations,12,13 1 duˆ iωpˆ= ρ a2 (13a) −1+(γ 1)f 0 0dx κ − 1 dpˆ iωUˆ = (1 f )A (13b) ν − − ρ dx 0 with, f = tanh(cid:18)(√1+2νi)/yω0(cid:19), f = tanh δ(ν1+/√i)Py0r , and y = h ν (1+i)y0 κ ((cid:16)1+i)y0 (cid:17) 0 2 √2ν/ω δν/√Pr (cid:18) (cid:19) (cid:16) (cid:17) canbespatiallyintegratedfromonehardend(ontheleftinfigure5),startingwithagivenarbitrarypressure value, all the way to the unknown impedance boundary at x = H. Here, pˆ represents discrete pressure perturbations and Uˆ represents the volume flow rate through a constant cross-sectional area A = h 1, × whereunitdepthinz axishasbeenassumed. Theseequationsaccountforviscouslossesatthewall(viathe thermoviscous functions f and f ) and represent the conservation of mass and energy combined (13a), and κ ν momentum (13b). Following this procedure, the impedance at the end of the duct for any given frequency is a result of the calculation and not an input as it would be for a classic acoustic eigenvalue problem. 6of14 AmericanInstituteofAeronauticsandAstronautics Figure 5: Inverse wave reconstruction via one dimensional spatial integration of Rott’s wave equations. The impedance at the end of the tube is calculated directly via Z(x;ω)= 1(pˆ +pˆ )A/Uˆ . 2 N−1 N N The cross-sectional averaged velocity at the end of the duct, Uˆ(x=H)=Uˆ can be used to recover the N full velocity profile in the wall-normal direction using the equation: Uˆ cosh(η) N uˆ(x=H,y)= 1 (14) A − cosh(η ) (cid:18) 0 (cid:19) iω iω η =y , and η(y)=y with 0 y <y . 0 0 0 ν ν ≤ r r The wall-normal profile of acoustic impedance at the end of the duct for each frequency, ω, can then be calculated as, 1(pˆ +pˆ ) Z(x,ω)= 2 N−1 N . (15) uˆ(x=H,y) Note that in the absence of viscosity, f and f should be explicitly set to zero and the equations (13) ν κ revert to those for isentropic acoustics. Figure6: Computational setupto analyzeasimple ductusingthe inverse HelmholtzSolver (iHS). Thegrid shown here is just meant for illustrative purposes. Figure 6 shows the two-dimensional computational set up used in the iHS runs, compared against Rott’s theory in the following subsections. The dimensions of the duct are H =4 mm and h=1 mm. Results are first obtained for an inviscid duct, with lateral slip walls for a range of frequencies in the ultrasonic regime, and then for a viscous duct with a viscosity of ν =0.05 m2/s. III.A. Inviscid Duct Figure 7 compares solutions obtained from the iHS against semi-analytical solutions obtained from Rott’s wave equations over several decades of frequencies for an inviscid duct. In the absence of shear stress, no wave attenuation should be present (if not numerical) and the acoustic velocity profile is flat everywhere. In figure 7, the periodic admittance peaks represent resonant frequencies at which theduct length is a multiple integeroftheacousticwavelength,andwherethesurfacenormaladmittanceattheendoftheductdiverges. 7of14 AmericanInstituteofAeronauticsandAstronautics Thequalityoftheagreementdecreasesathigherfrequencieswherethenumericalsolutioneventuallybecomes under-resolved in the axial direction. 3 2 k C B 1 I Y k 0 0 a 0 ρ 1 2 0 10 20 30 40 50 60 70 80 fH/a 0 Figure 7: Admittance magnitude obtained from the iHS (—) and from Rott’s wave equations ( ) for an ◦ inviscid two-dimensional duct. III.B. Viscous Duct 0.05 0.09 0.08 0.04 0.07 0.06 0.03 0.05 0.04 0.02 0.03 0.01 0.02 0.01 0.00 0.00 -0.5 -0.3 -0.1 0.1 0.3 0.5 -0.5 -0.3 -0.1 0.1 0.3 0.5 0.35 3.0 0.30 2.5 0.25 2.0 0.20 1.5 0.15 1.0 0.10 0.05 0.5 0.00 0.0 -0.5 -0.3 -0.1 0.1 0.3 0.5 -0.5 -0.3 -0.1 0.1 0.3 0.5 Figure 8: Comparisons of results obtained from the iHS (—) and Rott’s wave equations ( ) for a two- ◦ dimensional viscid duct at different frequencies. Figure 8 shows how solutions from the iHS compare with Rott’s thermoacoustic theory for viscous cases. It is observed that at lower frequencies there is poor matching with the semi-analytical solution, primarily due to the fact that Rott’s thermoacoustic theory chooses to treat the duct aspect ratio as very large, not accounting for edge effects, whereas in our case we have H/h = 4 . At high frequencies, edge effects are confined near the hard end and away from the unknown impedance boundary, returning perfect matching with Rott’s solutions. In general, it was found that there is good agreement with Rott’s solutions at frequencies where the Stokes boundary layer thickness was less than 20 % of the duct width for this specific geometry. 8of14 AmericanInstituteofAeronauticsandAstronautics IV. Two-Dimensional Toy Cavity Figure9: Sampleunstructuredcomputationalmeshdiscretizingidealizedtoycavitygeometryanalyzedusing the iHS. Inthissection,wepresenttheanalysisofanidealizedUACgeometry(figure9)usingtheinverseHelmholtz solver developed in the previous sections. The impinging wave is taken to be a planar wave that is first oriented normally (figure 10a), and then at an angle of 135o measured with respect to the positive x axis (figure 10b). (a) (b) Figure 10: Pressure waveform and magnitude of normal admittance for a 700 kHz planar wave for θ = 90o (a) and θ =135o (b), where θ is the wave’s incidence angle measured with respect to the positive x axis. Figure10showsthepressurecontoursandsurfacenormalspecificadmittancemagnitudeattheimpedance boundaryoftheidealizedgeometry,inresponsetoaplanarwave. Inbothcases,theviscousboundarylayers 9of14 AmericanInstituteofAeronauticsandAstronautics responsible for wave attenuation are resolved. While the admittance at the edges–where the IBC intersects with the hard, no-slip wall–is numerically almost the same for both cases, the intermediate profiles are different. In the normal case, the admittance is lowest near the center obstacle, due to interference from waves reflected off the obstacles on either side of it. For the angled case, the minima occurs approximately in the region bounded by the three obstacles on the right. The peak admittance in the angled case is significantly higher than that in the case of normal incidence, due to the fact that the angled wave is able to traverse more distance before it encounters the obstacles on the right. V. Two-Dimensional Analysis of C/C UAC In this section, we employ the iHS to analyze the impedance of C/C UACs. Figure11: Two-dimensionalslicesofasubvolumeofthethreedimensionalC/CUACinsertshowninfigure1, each 1.5 mm 1.5 mm in size. × Carbon fiber reinforced carbon ceramic (C/C) belongs to a class of materials commonly known as ce- ramic matrix composites (CMCs). C/C CMCs have excellent thermal resistance and a low coefficient of thermal expansion, rendering them ideal for structural uses at high temperatures typically encountered in re-entry vehicles.14,15 Moreover, porous C/C surfaces have been shown to damp second-mode instabilities in hypersonic boundary layers, delaying the laminar-to-turbulent trasition.16–18 C/CoriginatesasasemifinishedproductfromtheC/C-SiCmanufacturingprocessviatheliquidsilicone infiltration route as described by Krenkel.19 The entire process is depicted in figure 13. The manufacturing comprisesthefabricationofagreenbodyofcarbonfiberreinforcedplastic(CFRP),inthiscaseconsistingof carbon fibers impregnated with a phenolic resin processed via autoclave technique. The green body is cured and subsequently pyrolized at temperatures up to 1650C. This process step converts the matrix component from phenolic resin to amorphous carbon. During the following cool-down the shrinkage of the pyrolized matrix is hindered by the carbon fibers resulting in a pattern of micro cracks. The micro crack system of the C/C CMC used here is shown in figure 11. The tomographic images were taken with the aid of a GE Phoenix Nanotomc 180nF computer tomograph, using a voxel size of 1µm and (cid:13) a volume of interest (VOI) of 15003 voxels. It can be seen that the C/C is based on a CRFP with a 0—90 lay-up of carbon fiber fabric. Therefore, the material exhibits an orthotropic layout and behaviour, resulting in different thermo-mechanical and 10of14 AmericanInstituteofAeronauticsandAstronautics