Acoustics/Interface

In [1]:

import matplotlib.pyplot as plt
%matplotlib inline
import numpy as np
import flowdyn.mesh           as mesh
import flowdyn.xnum           as Xnum
import flowdyn.integration    as Tnum
import flowdyn.modelphy.euler as euler
import flowdyn.modeldisc      as modeldisc
#
plt.rcParams['figure.dpi'] = 120
plt.rcParams['savefig.dpi'] = 120
plt.rcParams["animation.html"] = "jshtml"  # for matplotlib 2.1 and above, uses JavaScript

import matplotlib.pyplot as plt %matplotlib inline import numpy as np import flowdyn.mesh as mesh import flowdyn.xnum as Xnum import flowdyn.integration as Tnum import flowdyn.modelphy.euler as euler import flowdyn.modeldisc as modeldisc # plt.rcParams['figure.dpi'] = 120 plt.rcParams['savefig.dpi'] = 120 plt.rcParams["animation.html"] = "jshtml" # for matplotlib 2.1 and above, uses JavaScript

Set up a acoustic wave propagating through an interface¶

A planar acoustic wave is defined by a (gaussian) profile of u. Using (non linear) characteristic invariant, the associated perturbation of p is computed to ensure a forward propagating wave. Density is then computed to ensure a isentropic flow. But an interface (discontinuity of density) is defined at x=6.

In [2]:

meshsim  = mesh.unimesh(ncell=200,  length=1.)
#meshref  = unimesh(ncell=1000, length=1.)

model = euler.model()

bcL  = { 'type': 'sym' } # not physical but can work
bcR  = { 'type': 'sym' } # for wall
xnum = Xnum.muscl(Xnum.vanalbada) ; flux = 'hllc'
#xnum = Xnum.extrapol1() ; flux = 'centered'

rhs = modeldisc.fvm(model, meshsim, numflux=flux, num=xnum, bcL=bcL, bcR=bcR)
solver = Tnum.lsrk26bb(meshsim, rhs)

# computation
#
nsol    = 100
endtime = .8
cfl     = .8

# initial functions
def fu(x):
    vmag = .01 #; k = 10.
    return vmag*np.exp(-500*(x-.2)**2) #*np.sin(2*np.pi*k*x)
def fp(x): # gamma = 1.4
    return (1. + .2*fu(x))**7.  # satisfies C- invariant to make only C+ wave
def frho(x):
    rhoratio = 10.
    return 1.4 * ( fp(x)**(1./1.4)*(x<.6) + rhoratio*(x>.6) )

xc    = meshsim.centers()
finit = rhs.fdata_fromprim([ frho(xc), fu(xc), fp(xc) ]) # rho, u, p

fsol = solver.solve(finit, cfl, np.linspace(0., endtime, nsol+1))
solver.show_perf()

# Figure / Plot

varname='pressure' # mach, pressure, entropy
ttime = [ fsol[i].time for i in range(nsol+1) ]
xx, xt = np.meshgrid(xc, ttime)
solgrid = [ fsol[i].phydata(varname) for i in range(nsol+1) ]
vmin, vmax = np.min(solgrid), np.max(solgrid)
#
# Figure / Plot of final results
fig, (ax1,ax2) = plt.subplots(1, 2, figsize=(12,4))
ax1.set_ylabel(varname) ; ax1.set_ylim(vmin, vmax)
ax1.grid(linestyle='--', color='0.5')
finit.plot(varname, 'k-', axes=ax1)
line1, = fsol[-1].plot(varname, 'b-', axes=ax1)
ax2.set_ylabel('t') ; ax2.set_xlim(0., 1.)
#ax2.grid(linestyle='--', color='0.5')
flood  = ax2.contour(xx, xt, solgrid, np.linspace(vmin, vmax, 50))
line2, = ax2.plot([0., 10.], [ttime[-1], ttime[-1]], 'k--')
plt.show()

meshsim = mesh.unimesh(ncell=200, length=1.) #meshref = unimesh(ncell=1000, length=1.) model = euler.model() bcL = { 'type': 'sym' } # not physical but can work bcR = { 'type': 'sym' } # for wall xnum = Xnum.muscl(Xnum.vanalbada) ; flux = 'hllc' #xnum = Xnum.extrapol1() ; flux = 'centered' rhs = modeldisc.fvm(model, meshsim, numflux=flux, num=xnum, bcL=bcL, bcR=bcR) solver = Tnum.lsrk26bb(meshsim, rhs) # computation # nsol = 100 endtime = .8 cfl = .8 # initial functions def fu(x): vmag = .01 #; k = 10. return vmag*np.exp(-500*(x-.2)**2) #*np.sin(2*np.pi*k*x) def fp(x): # gamma = 1.4 return (1. + .2*fu(x))**7. # satisfies C- invariant to make only C+ wave def frho(x): rhoratio = 10. return 1.4 * ( fp(x)**(1./1.4)*(x<.6) + rhoratio*(x>.6) ) xc = meshsim.centers() finit = rhs.fdata_fromprim([ frho(xc), fu(xc), fp(xc) ]) # rho, u, p fsol = solver.solve(finit, cfl, np.linspace(0., endtime, nsol+1)) solver.show_perf() # Figure / Plot varname='pressure' # mach, pressure, entropy ttime = [ fsol[i].time for i in range(nsol+1) ] xx, xt = np.meshgrid(xc, ttime) solgrid = [ fsol[i].phydata(varname) for i in range(nsol+1) ] vmin, vmax = np.min(solgrid), np.max(solgrid) # # Figure / Plot of final results fig, (ax1,ax2) = plt.subplots(1, 2, figsize=(12,4)) ax1.set_ylabel(varname) ; ax1.set_ylim(vmin, vmax) ax1.grid(linestyle='--', color='0.5') finit.plot(varname, 'k-', axes=ax1) line1, = fsol[-1].plot(varname, 'b-', axes=ax1) ax2.set_ylabel('t') ; ax2.set_xlim(0., 1.) #ax2.grid(linestyle='--', color='0.5') flood = ax2.contour(xx, xt, solgrid, np.linspace(vmin, vmax, 50)) line2, = ax2.plot([0., 10.], [ttime[-1], ttime[-1]], 'k--') plt.show()

cpu time computation (202 it) : 0.646s
  16.00 µs/cell/it

In [3]:

import matplotlib.animation as anim
#
def animate(k):
    fsol[k].set_plotdata(line1, varname)
    line2.set_data([0., 10.], [ttime[k], ttime[k]])
    return line1, line2

ani = anim.FuncAnimation(fig=fig, func=animate, frames=range(nsol+1), interval=100)#, blit=True)
ani

import matplotlib.animation as anim # def animate(k): fsol[k].set_plotdata(line1, varname) line2.set_data([0., 10.], [ttime[k], ttime[k]]) return line1, line2 ani = anim.FuncAnimation(fig=fig, func=animate, frames=range(nsol+1), interval=100)#, blit=True) ani

Out[3]:

Once Loop Reflect

Comments¶

The interaction of the acoustic wave with the interface produces a reflected wave and a transmitted one. Waves velocities are well captured. However, the propagation suffers (a non expected) dissipation.