Rayleigh flow

In [1]:

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

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

Theoretical computation of needed power with aerokit¶

In [2]:

import aerokit.aero.Isentropic as Is
import aerokit.aero.MassFlow   as mf
import aerokit.aero.Rayleigh   as ray
#
M0 = .2   # expected upstream Mach number
M1 = .95  # expected downstream Mach number
#
rTt0 = 1.
p1   = 1.
Tt_ratio = ray.Ti_Ticri(M1)/ray.Ti_Ticri(M0)
Pt_ratio = ray.Pi_Picri(M1)/ray.Pi_Picri(M0)
NPR      = Is.PiPs_Mach(M1) / Pt_ratio
Power    = 3.5*(Tt_ratio-1.)*rTt0*(NPR*p1/np.sqrt(rTt0)*mf.WeightMassFlow(M0))
print(Tt_ratio, NPR, Power)

import aerokit.aero.Isentropic as Is import aerokit.aero.MassFlow as mf import aerokit.aero.Rayleigh as ray # M0 = .2 # expected upstream Mach number M1 = .95 # expected downstream Mach number # rTt0 = 1. p1 = 1. Tt_ratio = ray.Ti_Ticri(M1)/ray.Ti_Ticri(M0) Pt_ratio = ray.Pi_Picri(M1)/ray.Pi_Picri(M0) NPR = Is.PiPs_Mach(M1) / Pt_ratio Power = 3.5*(Tt_ratio-1.)*rTt0*(NPR*p1/np.sqrt(rTt0)*mf.WeightMassFlow(M0)) print(Tt_ratio, NPR, Power)

5.751213860870039 2.2040855280805833 8.468636339826547

In [3]:

nx = 100
lx = 6.
meshsim  = unimesh(ncell=nx,  length=10.)

def fenergy(x, q):
    return +Power/lx*(x>2.)*(x<2.+lx)

model = euler.model(source=[None, None, fenergy])
bcL = { 'type': 'insub',  'ptot': NPR*p1, 'rttot': rTt0 }
bcR = { 'type': 'outsub', 'p': p1 }

rhs = modeldisc.fvm(model, meshsim, muscl(vanalbada), 
      bcL=bcL, bcR=bcR)
solver = rk3ssp(meshsim, rhs)

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

finit = rhs.fdata_fromprim([  1., .3, 1. ]) # rho, u, p

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

solver.show_perf()

nx = 100 lx = 6. meshsim = unimesh(ncell=nx, length=10.) def fenergy(x, q): return +Power/lx*(x>2.)*(x<2.+lx) model = euler.model(source=[None, None, fenergy]) bcL = { 'type': 'insub', 'ptot': NPR*p1, 'rttot': rTt0 } bcR = { 'type': 'outsub', 'p': p1 } rhs = modeldisc.fvm(model, meshsim, muscl(vanalbada), bcL=bcL, bcR=bcR) solver = rk3ssp(meshsim, rhs) # computation # nsol = 100 endtime = 100. cfl = 1. finit = rhs.fdata_fromprim([ 1., .3, 1. ]) # rho, u, p fsol = solver.solve(finit, cfl, np.linspace(0., endtime, nsol+1)) solver.show_perf()

cpu time computation (5029 it) : 5.545s
  11.03 µs/cell/it

In [4]:

# Figure / Plot of final results
varlist=['htot', 'mach', 'enthalpy']
nvar = len(varlist)
lines=[None]*nvar # dummy init
fig, ax = plt.subplots(1, nvar, figsize=(5*nvar,4))
for i in range(nvar):
    varname = varlist[i]
    ax[i].set_ylabel(varname) ; ax[i].set_ylim(0., 1.2*np.max(fsol[-1].phydata(varname)))
    ax[i].grid(linestyle='--', color='0.5')
    lines[i], = fsol[-1].plot(varname, 'k-', axes=ax[i])

# Figure / Plot of final results varlist=['htot', 'mach', 'enthalpy'] nvar = len(varlist) lines=[None]*nvar # dummy init fig, ax = plt.subplots(1, nvar, figsize=(5*nvar,4)) for i in range(nvar): varname = varlist[i] ax[i].set_ylabel(varname) ; ax[i].set_ylim(0., 1.2*np.max(fsol[-1].phydata(varname))) ax[i].grid(linestyle='--', color='0.5') lines[i], = fsol[-1].plot(varname, 'k-', axes=ax[i])

In [5]:

import matplotlib.animation as anim
#
def animate(k):
    for i in range(nvar):
        varname = varlist[i]
        fsol[k].set_plotdata(lines[i], varname)
    return lines

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

import matplotlib.animation as anim # def animate(k): for i in range(nvar): varname = varlist[i] fsol[k].set_plotdata(lines[i], varname) return lines ani = anim.FuncAnimation(fig=fig, func=animate, frames=range(nsol+1), interval=100, blit=True) ani

Out[5]:

Once Loop Reflect