TaSMET example: a resonant tube driven by a piston

Step 1: import the code models into the Python namespace



In [ ]:

    
from TaSMET import *
import mpld3; mpld3.enable_notebook()

Step 2: Define the global configuration, initialize a TaSystem



In [ ]:

    
Nf=6 # The number of harmonics to solve for
freq=343/4 # The oscillating frequency
omg=2*pi*freq # Radian frequency
# Define a global configuration
gc=Globalconf_airSTP(Nf,freq)

# And: create a TaSystem object
sys=TaSystem(gc)

Generate a Piston object



In [ ]:

    
Km=1e-4 # Stiffness of the spring
Cm=0    # Mechanical damping
Mp=0.1   # Mass of the spring
R=0.02 # Radius of the tube of 1 cm
Sl=Sr=pi*R**2 # Left and right surface area of the piston

Ma_low=0.001
Ma_high=0.02
Ma=Ma_low
u1=Ma*gc.c0()
x1=u1/omg

Vr=1.1*x1*Sr #
Vl=10*Vr # Very compliant left volume
# Generate a piston configuration
pc=PistonConfiguration(Sl,Sr,Vl,Vr,Mp,Km,Cm)

# Finally, build the piston
piston=Piston(pc)
# To refer to the piston, we need an ID:
piston.setID('piston')

Generate a tube model

First step: generate 1D grid

Example: linear grid



In [ ]:

    
lineargrid=LinearGrid(50,1)
plot(lineargrid.getx(),'o'); xlabel('N'); ylabel('x');

Better: boundary layer grid



In [ ]:

    
blgrid=BlGrid(1,gc.deltanu0()/30,0.01)
plot(blgrid.getx(),'o'); xlabel('N'); ylabel('x');

Generate a geometry



In [ ]:

    
# Choices: ConeTube, CylindricalTube, VertPlates
geom=CylindricalTube(blgrid.getx(),R,True)

Finally: generate the LaminarDuct



In [ ]:

    
tube=LaminarDuct(geom)
# Allow for exchanging heat with the wall
tube.setInsulated(False)
# Set an ID to refer to it
tube.setID('tube')

Connect the piston to the tube, apply boundary conditions

Connect the tube to the piston



In [ ]:

    
# A connector to connect the tube on the left side to the piston
# on the right side
ductpistoncon=DuctPistonConnector('tube',left,'piston',right)

Apply a wall boundary condition on the right side of the tube



In [ ]:

    
wallbc=IsoTWall('tube',right,var(gc,293.15),True)

Apply a mechanical boundary condition on the piston



In [ ]:

    
forcevar=var(gc)
forcevar.setadata(1,1000) # Apply a force of 10N at the first harmonic: F_1=10 N
#forcevar.setadata(3,8) # Apply a force of 8N at tha second harmonic: F_2=8 N
mechbc=MechBc('piston',Varnr_F,forcevar)

# Or, if we want a kinematic constraint on the piston
xvar=var(gc)
xvar.setadata(1,x1)

mechbc=MechBc('piston',Varnr_x,xvar)

Add all segments and boundary conditions to the TaSystem



In [ ]:

    
setTaSMETTracer(40)
# Segments
sys+=tube
sys+=piston
# Boundary conditions
sys+=mechbc
sys+=ductpistoncon
sys+=wallbc

Show if everything went well so far



In [ ]:

    
#sys.show(4)

Create a solver object and solve the system!



In [ ]:

    
sol=Solver()
sol.setMaxiter(100)
sol.setMaxdampfac(0.5)



In [ ]:

    
sol.solve(sys)
semilogy(sol.getSp().Funer())

Do the postprocessing!



In [ ]:

    
piston=sys.getPiston('piston')
tube=sys.getDuct('tube')
x=tube.getx()

The tube

Plot the pressure in the tube at different time instances



In [ ]:

    
T=1/freq # Obtain the period
Np=15 # Number of time instance to plot
figure(figsize=(12,4))
hold('on')
for i in range(Np):
    plot(x,tube.getValueT(Varnr_p,i*T/(Np+1))) # Obtain the value of the pressure in the tube at a certain time instance
xlabel('x')
ylabel('p_t')



In [ ]:

    
figure(figsize=(12,4))
hold('on')
for i in range(Np):
    plot(x,tube.getValueT(Varnr_U,i*T/(Np+1))) # Obtain the value of the pressure in the tube at a certain time instance
xlabel('x')
ylabel('U_t')

Looking at phasors



In [ ]:

    
figure(figsize=(12,3))
T1=tube.getValueC(Varnr_T,1)
plot(x,T1.real)
xlabel('x')
ylabel('T_1')

Looking at phasors



In [ ]:

    
figure(figsize=(12,3))
T2=tube.getValueC(Varnr_T,2)
plot(x,T2.real)



In [ ]:

    
mH=tube.getValue(Varnr_mH,0)
Q=tube.getValue(Varnr_Q,0)
#figure(figsize=(12,3))
#plot(x,mH+Q)

Variables in the piston



In [ ]:

    
figure(figsize=(12,3))
subplot(121)
plot(piston.Vr().timeResponse())
subplot(122)
plot(piston.pr().timeResponse())