TMmodel infinite layer examples
Absorber
# Example with 3 layups for sound absorption
import numpy as np
import matplotlib.pyplot as plt
import pyva.models as mds
import pyva.systems.infiniteLayers as iL
# my packages
import pyva.properties.materialClasses as matC
plt.close('all')
omega = 2*np.pi*np.logspace(1,4,100)
# default fluid
air = matC.Fluid(eta = 0.0)
fibre1 = matC.EquivalentFluid(porosity = 0.98, \
flow_res = 25000.,\
tortuosity = 1.02, \
length_visc = 90.e-6, \
length_therm = 180.e-6,\
rho_bulk = 0.98*1.20 + 30. , \
rho0 = 1.208, \
dynamic_viscosity = 1.81e-5 )
# Thickness
h1 = 0.1
h2 = 0.2
angles = np.linspace(0,np.pi/2)
angle45 = np.pi/4
kx45 = air.wavenumber(omega)*np.sin(angle45)
#kx = mC.DataAxis(kx,typestr='wavenumber')
z0 = air.z0
fibre_10cm = iL.FluidLayer(h1,fibre1)
fibre_20cm = iL.FluidLayer(h2,fibre1)
perforate = iL.PerforatedLayer(0.005, 0.001, distance = 0.02)
TMM_fibre_10 = mds.TMmodel((fibre_10cm,))
TMM_fibre_20 = mds.TMmodel((fibre_20cm,))
TMM_perf_fibre_20 = mds.TMmodel((perforate, fibre_10cm,))
alpha_fibre_10 = TMM_fibre_10.absorption_diffuse(omega,theta_max=np.pi/2,signal = False)
alpha_fibre_20 = TMM_fibre_20.absorption_diffuse(omega,theta_max=np.pi/2,signal = False)
alpha_perf_fibre_10 = TMM_perf_fibre_20.absorption_diffuse(omega,theta_max=np.pi/2,signal = False)
# %% plot1
# plot impedance of fibre results for publishing
plt.figure(1)
plt.plot(omega,alpha_fibre_20,label='20cm fibre')
plt.plot(omega,alpha_fibre_10,label='10cm fibre')
plt.plot(omega,alpha_perf_fibre_10,label='perforate + 10cm fibre' )
plt.xscale('log')
plt.xlabel('$\omega/$s$^{-1}$')
plt.ylabel('absorption')
plt.legend()
plt.savefig('../source/images/TMM_absorber_abs_diffuse.png')
Double Walls
# -*- coding: utf-8 -*-
#
# Example with 3 tubes showing the capabilities of DynamicMatrix
import numpy as np
import matplotlib.pyplot as plt
import pyva.models as mds
import pyva.systems.infiniteLayers as iL
import pyva.properties.structuralPropertyClasses as stPC
# my packages
import pyva.properties.materialClasses as matC
plt.close('all')
omega = 2*np.pi*np.logspace(1.5,4,1000)
# default fluid
air = matC.Fluid(eta = 0.01)
rho_bulk = 0.98*1.20 + 30.
fibre1 = matC.EquivalentFluid(porosity = 0.98, \
flow_res = 25000.,\
tortuosity = 1.02, \
length_visc = 90.e-6, \
length_therm = 180.e-6,\
rho_bulk = rho_bulk , \
rho0 = 1.208, \
dynamic_viscosity = 1.81e-5 )
# Create props
alu = matC.IsoMat()
alu1mm = stPC.PlateProp(0.001,alu)
# Thickness
omega1 = 6000
k1 = np.real(air.wavenumber(omega1))
angles = np.linspace(0,np.pi/2)
#kx = k1*np.sin(angles)
#kx = mC.DataAxis(kx,typestr='wavenumber')
theta = 50/180*np.pi
kx = air.wavenumber(omega)*np.sin(theta)
z0 = air.z0
# Fluid layer
air_5cm = iL.FluidLayer(0.05,air)
fibre_5cm = iL.FluidLayer(0.05,fibre1)
# Limp mass layer
heavy_1kg = iL.MassLayer(0.001, 1000)
heavy_2kg7 = iL.MassLayer(0.001, 2700)
# Plate layer
iL_alu1mm = iL.PlateLayer(alu1mm,)
#heavy_1kg = iL.PlateLayer(alu03mm)
T_alu = mds.TMmodel((iL_alu1mm,))
T_alu_air_mass = mds.TMmodel((iL_alu1mm,air_5cm,heavy_1kg))
T_alu_fibre_mass10 = mds.TMmodel((iL_alu1mm,fibre_5cm,heavy_1kg))
T_alu_fibre_mass27 = mds.TMmodel((iL_alu1mm,fibre_5cm,heavy_2kg7))
tau_alu = T_alu.transmission_diffuse(omega,signal = False)
tau_alu_air_mass = T_alu_air_mass.transmission_diffuse(omega,theta_step=np.pi/1000,signal = False)
tau_alu_fibre_mass10 = T_alu_fibre_mass10.transmission_diffuse(omega,signal = False)
tau_alu_fibre_mass27 = T_alu_fibre_mass27.transmission_diffuse(omega,signal = False)
#%% fig1
plt.figure(1)
plt.plot(omega,-10*np.log10(tau_alu),label='alu 1mm')
plt.plot(omega,-10*np.log10(tau_alu_air_mass),label='alu-air5cm-1kg')
plt.plot(omega,-10*np.log10(tau_alu_fibre_mass10),label='alu-fibre5cm-1kg')
plt.plot(omega,-10*np.log10(tau_alu_fibre_mass27),label='alu-fibre5cm-2.7kg')
plt.xscale('log')
plt.xlabel('$\omega/$s$^{-1}$')
plt.ylabel('$TL$')
plt.legend()
plt.savefig('../source/images/TMM_DW_transmission.png')
#%%fig 2
Material Parameter Determination
# Example of JCA-parameter identification
#
# See Atalla Y, Panneton, R.; Inverse Acoustical Characterizatin of Open Cell Porous Media Using Impedance Tube Measurements
#
# Taken from figure 7 - not 6
# Numerics
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import dual_annealing
# pyva packages
import pyva.properties.materialClasses as matC
import pyva.models as mds
import pyva.systems.infiniteLayers as iL
plt.close('all')
omega = 2*np.pi*np.geomspace(200,2000,100)
freq = omega/2/np.pi
# default fluid
# air = matC.Fluid(eta = 0.0)
# Frequency range
f_min = 500
# select test
test = input('Enter test selector 1-4!')
test = int(test)
# switch if fibres are calculated based on limp model
fibre_sw = True # False means rigid
# name test
test_str = ('foam1','foam2','fibrous1','fibrous2')
# Given test data
if test == 1:
# Foam 1
flow_res = 4971.
porosity = 0.97
rho_bulk = 21.6
h = 49.92e-3
#h = 98.7e-3
limp = False
# reference values
tor_ref = 1.31
Lam_visc_ref = 123.19e-6
Lam_term_ref = 289.54e-6
# Environment
T = 21.5+273.15
P0 = 0.992
elif test == 2:
# Foam 2
flow_res = 8197.
porosity = 0.95
rho_bulk = 23.9
h = 34.35e-3
#h = 68.77e-3
limp = False
# reference values
tor_ref = 1.42
Lam_visc_ref = 133.10e-6
Lam_term_ref = 212.65e-6
# Environment
T = 24.6+273.15
P0 = 0.997
elif test == 3:
# Fibrous 1
flow_res = 21235.
porosity = 0.94
rho_bulk = 89.6
h = 23.37e-3
#h = 46.83e-3
limp = fibre_sw
# reference values
tor_ref = 1.00
Lam_visc_ref = 48.62e-6
Lam_term_ref = 114.39e-6
# Environment
T = 22.6+273.15
P0 = 1.002
elif test == 4:
# Fibrous 2
flow_res = 50470.
porosity = 0.89
rho_bulk = 150.0
h = 37.38e-3
#h = 18.98e-3
limp = fibre_sw
# reference values
tor_ref = 1.00
Lam_visc_ref = 41.51e-6
Lam_term_ref = 114.27e-6
# Environment
T = 24.6+273.15
P0 = 0.977
# Set air under test conditions
air_test = matC.Fluid.air(T,P0)
# Create function from model
def fibre_fit(tor,Lam_visc,Lam_term):
"""
Function of fibre material with 3 parameters for optimisation
Parameters
----------
tor : float
tortuosity.
Lam_visc : float
characteristic viscous length.
Lam_term : float
characteristic termal length.
Returns
-------
EquivalentFluid object with 3 paramters
"""
return matC.EquivalentFluid(flow_res, porosity , tor, rho_bulk, \
Lam_visc, Lam_term,\
rho0 = air_test.rho0, \
c0 = air_test.c0,\
dynamic_viscosity = air_test.dynamic_viscosity,\
Cp = air_test.Cp, heat_conductivity=air_test.heat_conductivity, \
limp = limp )
def layer_fit(tor,Lam_visc,Lam_term):
"""
Function of fibre material layer with 3 parameters for optimisation
Parameters
----------
tor : float
tortuosity.
Lam_visc : float
characteristic viscous length.
Lam_term : float
characteristic termal length.
Returns
-------
TMmodel object with 3 paramters
"""
return mds.TMmodel((iL.FluidLayer(h,fibre_fit(tor,Lam_visc,Lam_term)),))
def impedance_fit(f,tor,Lam_visc,Lam_term):
"""
Function of fibre material layer with 3 parameters for optimisation
Parameters
----------
f : ndarray of float
frequnecy.
tor : float
tortuosity.
Lam_visc : float
characteristic viscous length.
Lam_term : float
characteristic termal length.
Returns
-------
Signal
surface impedance.
"""
return layer_fit(tor,Lam_visc,Lam_term).impedance(2*np.pi*f,0.).ydata.flatten()
z0 = air_test.z0
# Import test data
f_test,Zs_re,Zs_im = np.loadtxt ('.//data//'+test_str[test-1]+'.csv',
unpack = True,
usecols = (0,1,2), skiprows = 1,
delimiter = ',')
# create absolute values
Zs = (Zs_re+1j*Zs_im)*z0
# index for frequency range selection
i_freq = f_test >= f_min
# Create cost function with given test data
def cost_function(x):
"""
Cost function of mean square residual
Parameters
----------
x : array of 3 parameters
(tortuosity, visc_length, term_length).
Returns
-------
float
sum of squared difference between test and model impedance.
"""
return np.sum(np.abs(impedance_fit(f_test[i_freq],x[0],x[1],x[2])-Zs[i_freq])**2)
# set bounds
lw = [1.,1.e-6,1.e-6] # lower bounds
up = [4.,4.e-4,4.e-4] # upper bounds
bounds=list(zip(lw, up))
# set start value
x0 = [1., 30.1e-6, 60.1e-6]
#res = least_squares(cost_function, x0,bounds = bounds)
res = dual_annealing(cost_function, bounds = bounds)
#%% plot1
# plot impedance of fibre results for publishing
plt.figure(1)
plt.plot(f_test,Zs_re,'d:',label='Re')
plt.plot(f_test,Zs_im,'d:',label='Im')
1
# check fit function f,por,sigma,Lam_visc,Lam_term
plt.plot(freq,np.real(impedance_fit(freq, *res.x)/z0), label = 'pyva_Re')
plt.plot(freq,np.imag(impedance_fit(freq, *res.x)/z0), label = 'pyva_Im')
# check Atalla solution ,por,sigma,Lam_visc,Lam_term
plt.plot(freq,np.real(impedance_fit(freq, tor_ref, Lam_visc_ref,Lam_term_ref)/z0), label = 'Atalla_Re')
plt.plot(freq,np.imag(impedance_fit(freq, tor_ref, Lam_visc_ref,Lam_term_ref)/z0), label = 'Attala_Im')
X0 = 550; X1 = 900; Y0 = -3.6; DY = -0.6
plt.text(X0, Y0, "Fit pyva:", fontsize=12)
plt.text(X1, Y0, r"$\alpha_\infty = {0:.2f}$".format(res.x[0]), fontsize=12)
plt.text(X1, Y0+DY, r"$\Lambda = {0:.2f}\mu$m".format(res.x[1]*1e6), fontsize=12)
plt.text(X1, Y0+2*DY, r"$\Lambda' = {0:.2f}\mu$m".format(res.x[2]*1e6), fontsize=12)
plt.text(X0, Y0+3*DY, "Atalla:", fontsize=12)
plt.text(X1, Y0+3*DY, r"$\alpha_\infty = {0:.2f}$".format(tor_ref), fontsize=12)
plt.text(X1, Y0+4*DY, r"$\Lambda = {0:.2f}\mu$m".format(Lam_visc_ref*1e6), fontsize=12)
plt.text(X1, Y0+5*DY, r"$\Lambda' = {0:.2f}\mu$m".format(Lam_term_ref*1e6), fontsize=12)
#plt.xscale('log')
plt.xlim((200,2000))
plt.ylim((-7,3))
plt.xlabel('$f/$Hz')
plt.ylabel('$Z_s$')
plt.title(test_str[test-1])
plt.legend()
plt.savefig('../source/images/atalla_JCA_parameter_'+test_str[test-1]+'.png')
Absorber Design with Poroelastic Materials
# -*- coding: utf-8 -*-
"""
Allard example of section 11.7.2 of [All2009]_
Example with lay-up given in figure 11.16 and table 11.7.
Table 11.7 is not correct - the fibrous layer thickness is 12.5mm noz 1.25mm.
Figure 11.17 shows the resuls but in Hz not in kHz as in [All2009]_
Alexander Peiffer
"""
import numpy as np
import matplotlib.pyplot as plt
import pyva.models as mds
import pyva.systems.infiniteLayers as iL
import pyva.properties.structuralPropertyClasses as stPC
import pyva.properties.materialClasses as matC
plt.close('all')
# Frequency Range
freq = np.linspace(200,1400,100) # kHz
omega = 2*np.pi*freq
# Define Materials
air = matC.Fluid(eta = 0.0)
# Carpet as Poroelastic
carpet_solid = matC.IsoMat(E=20000.,rho0=60,nu=0.,eta=0.5)
carpet_mat = matC.PoroElasticMat(carpet_solid, \
flow_res = 5000., \
porosity = 0.99, \
tortuosity = 1., \
length_visc = 23.E-6, length_therm = 28.E-6)
# Impervious scree as PlateProp
screen_mat = matC.IsoMat(E=30000.,rho0 = 2000, nu=0.49)
screen_prop = stPC.PlateProp(0.003, screen_mat)
# Fibre as Poroelastic
fibre_solid = matC.IsoMat(E=100000.,rho0=60,nu=0.,eta=0.88)
fibre_mat = matC.PoroElasticMat(fibre_solid, \
flow_res = 33000., \
porosity = 0.98, \
tortuosity = 1.1, \
length_visc = 50.E-6, length_therm = 110.E-6)
# Air impedance for normalisation
z0 = air.z0
# Define infiniteLayers
carpet1 = iL.PoroElasticLayer(carpet_mat, 0.0035)
carpet2 = iL.PoroElasticLayer(carpet_mat, 0.0035)
screen = iL.ImperviousScreenLayer(screen_prop)
fibre = iL.PoroElasticLayer(fibre_mat, 0.0125)
# Create lay-up as TMmodel
TMM_layup = mds.TMmodel((carpet1,carpet2,screen,fibre))
# Calculate normal surface impedance
Z = TMM_layup.impedance_allard(omega,kx=0,signal = False)
# Calculate normal absorption
alpha0 = TMM_layup.absorption(omega,kx=0.0,signal = False,allard=True)
# Calculate diffure absorption
alpha_diff = TMM_layup.absorption_diffuse(omega,theta_max=np.pi*78/180,theta_step = np.pi/180,signal = False,allard=True)
#%% plot impedance of according to fig 11.17 of [All2009]_
# The frequency is in Hz not in kHz as given in fig 11.17
plt.figure(1)
plt.plot(freq,np.real(Z)/z0 ,label='Re')
plt.plot(freq,np.imag(Z)/z0 ,label='Im')
plt.xlabel('$f$/Hz')
plt.ylabel('$z/z0$')
plt.legend()
plt.savefig('../source/images/carpet_fibre_impedance.png')
#%% plot absorption
plt.figure(2)
plt.plot(freq,alpha0 ,label='normal')
plt.plot(freq,alpha_diff, label='diffuse')
plt.xlabel('$f$/Hz')
plt.ylabel( r'$\alpha$')
plt.legend()
plt.savefig('../source/images/carpet_fibre_absorption.png')
Acoustic transmission design with poroelastic foam and rubber
# pyva double wall exampple
# using Allards transfer matrix method
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import pyva.models as mds
import pyva.systems.infiniteLayers as iL
import pyva.properties.structuralPropertyClasses as stPC
import pyva.properties.materialClasses as matC
plt.close('all')
# default fluid
air = matC.Fluid(rho0=1.21,eta = 0.0)
z0 = air.z0
# Angles and frequncies
freq = np.geomspace(100, 16000,100)
omega = np.pi*2*freq
# Isotropic materials
alu = matC.IsoMat(eta = 0.1)
rubber = matC.IsoMat(E=2.6e6,rho0=1200,nu=0.49,eta=0.00)
# Melamin foam paramters
melamin_vac = matC.IsoMat(E=3.0e5,rho0=12.0,nu=0.4,eta=0.1) # Frame in vaccuum E6
melamin = matC.PoroElasticMat(melamin_vac, \
flow_res = 30000., \
porosity = 0.99, \
tortuosity = 1.01, \
length_visc = 250.E-6, length_therm = 550.E-6)
# plate properties
alu_1mm = stPC.PlateProp(0.001, alu)
rubber_2mm = stPC.PlateProp(0.002, rubber)
# test the foam as solid
foam_3cm_solid = stPC.PlateProp(0.03, melamin_vac)
# Foam Layers
iL_foam_3cm = iL.PoroElasticLayer(melamin, 0.03)
iL_foam_3cm_solid = iL.SolidLayer(foam_3cm_solid)
# rubber and alu as solid- and screen layer
iL_rubber_solid_2mm = iL.SolidLayer(rubber_2mm)
iL_rubber_imper_2mm = iL.ImperviousScreenLayer(rubber_2mm)
iL_alu_solid_1mm = iL.SolidLayer(alu_1mm)
iL_alu_imper_1mm = iL.ImperviousScreenLayer(alu_1mm)
# Mass of Fluid as gap
# iL_nothing = iL.MassLayer(1e-6,1.)
iL_nothing = iL.FluidLayer(0.00001,air)
# TMmpodel using the solid or impervious screen formulation
alu_melamin_rubber_solid = mds.TMmodel((iL_alu_solid_1mm,iL_foam_3cm,iL_rubber_solid_2mm))
alu_melamin_rubber_imper = mds.TMmodel((iL_alu_imper_1mm,iL_foam_3cm,iL_rubber_imper_2mm))
# Try other variations
alu_melamin_rubber_decoup = mds.TMmodel((iL_alu_imper_1mm,iL_nothing,iL_foam_3cm,iL_rubber_imper_2mm))
alu_melamin_rubber_foam_as_solid = mds.TMmodel((iL_alu_imper_1mm,iL_foam_3cm_solid,iL_rubber_imper_2mm))
tau_imper = alu_melamin_rubber_imper.transmission_diffuse(omega,theta_max=78/180*np.pi,allard=True,signal=False)
tau_solid = alu_melamin_rubber_solid.transmission_diffuse(omega,theta_max=78/180*np.pi,allard=True,signal=False)
tau_decoup = alu_melamin_rubber_decoup.transmission_diffuse(omega,theta_max=78/180*np.pi,allard=True,signal=False)
tau_foam_as_solid = alu_melamin_rubber_foam_as_solid.transmission_diffuse(omega,theta_max=78/180*np.pi,allard=True,signal=False)
#%% Plot1
plt.figure(1)
plt.plot(freq,-10*np.log10(tau_imper),label='impervious screen')
plt.plot(freq,-10*np.log10(tau_solid),label='solid')
plt.plot(freq,-10*np.log10(tau_decoup),label='decoup')
plt.plot(freq,-10*np.log10(tau_foam_as_solid),label='foam as solid')
plt.xlabel('f/Hz')
plt.ylabel('TL/dB')
plt.xscale("log")
plt.legend()
plt.savefig('../source/images/allard_DW_TL.png')
Perforated Absorber created with plates with bending wave motion
# Example with plate fibre layups and perforation
import numpy as np
import matplotlib.pyplot as plt
import pyva.models as mds
import pyva.systems.infiniteLayers as iL
import pyva.properties.structuralPropertyClasses as stPC
import pyva.properties.materialClasses as matC
plt.close('all')
omega = 2*np.pi*np.logspace(1,4,100)
# default fluid
air = matC.Fluid(eta = 0.0)
# Fibre material
fibre1 = matC.EquivalentFluid(porosity = 0.98, \
flow_res = 25000.,\
tortuosity = 1.02, \
length_visc = 90.e-6, \
length_therm = 180.e-6,\
rho_bulk = 0.98*1.20 + 30. , \
rho0 = 1.208, \
dynamic_viscosity = 1.81e-5 )
# plate material and properties
steel = matC.IsoMat(E=15.E10,rho0=7800.,nu=0.27)
steel5mm = stPC.PlateProp(0.005, steel)
nothing = stPC.PlateProp(0.000001, steel)
# infinite layers
il_fibre_10cm = iL.FluidLayer(0.1,fibre1)
il_perforate = iL.PerforatedLayer(0.005, 0.001, distance = 0.02)
il_steel_5mm = iL.PlateLayer(steel5mm)
il_nothing_perf = iL.PlateLayer(nothing,perforation = il_perforate)
il_steel_5mm_perf = iL.PlateLayer(steel5mm,perforation = il_perforate)
il_steel_5mm_s_perf = iL.ImperviousScreenLayer(steel5mm,perforation = il_perforate)
# layups definition
TMM_fibre_10 = mds.TMmodel((il_fibre_10cm,))
TMM_perf_fibre_10 = mds.TMmodel((il_perforate, il_fibre_10cm,))
TMM_steel_fibre_10 = mds.TMmodel((il_steel_5mm, il_fibre_10cm,))
TMM_steel_perf_fibre_10 = mds.TMmodel((il_steel_5mm_perf, il_fibre_10cm,))
TMM_steel_s_perf_fibre_10 = mds.TMmodel((il_steel_5mm_s_perf, il_fibre_10cm,))
TMM_nothing_perf_fibre_10 = mds.TMmodel((il_nothing_perf, il_fibre_10cm,))
# absorption calculation
alpha_steel_perf_fibre_10 = TMM_steel_perf_fibre_10.absorption_diffuse(omega,theta_max=np.pi/2,signal = False)
alpha_steel_s_perf_fibre_10 = TMM_steel_s_perf_fibre_10.absorption_diffuse(omega,theta_max=np.pi/2,signal = False,allard=True)
alpha_nothing_perf_fibre_10 = TMM_nothing_perf_fibre_10.absorption_diffuse(omega,theta_max=np.pi/2,signal = False)
alpha_fibre_10 = TMM_fibre_10.absorption_diffuse(omega,theta_max=np.pi/2,signal = False)
alpha_perf_fibre_10 = TMM_perf_fibre_10.absorption_diffuse(omega,theta_max=np.pi/2,signal = False)
alpha_steel_fibre_10 = TMM_steel_fibre_10.absorption_diffuse(omega,theta_max=np.pi/2,signal = False)
#%% plot1
plt.figure(1)
plt.plot(omega,alpha_fibre_10,label='10cm fibre')
plt.plot(omega,alpha_nothing_perf_fibre_10,':',label='1µm perforated nothing + 10cm fibre' )
plt.plot(omega,alpha_perf_fibre_10,label='perforate + 10cm fibre' )
#plt.plot(omega,alpha_steel_fibre_10,label='5mm steel + 10cm fibre' )
plt.plot(omega,alpha_steel_perf_fibre_10,'<',label='5mm perforated steel + 10cm fibre' )
plt.plot(omega,alpha_steel_s_perf_fibre_10,'.',label='5mm perforated solid steel + ...' )
plt.xscale('log')
plt.ylim((0,1.4))
plt.xlabel('$\omega/$s$^{-1}$')
plt.ylabel('absorption')
plt.legend(loc=2)
plt.savefig('../source/images/TMM_perforate_test_abs_diffuse.png')