-
Notifications
You must be signed in to change notification settings - Fork 19
/
spectra.py
93 lines (86 loc) · 2.8 KB
/
spectra.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
# -*- coding: utf-8 -*-
"""
Created on Tuesday May 29 13:08:01 2018
@author: Tony Saad
"""
# !/usr/bin/env python
from scipy import interpolate
import numpy as np
from numpy import pi
class file_spectrum:
def __init__(self,specfname):
filespec = np.loadtxt(specfname)
kfile=filespec[:,0]
efile=filespec[:,1]
self.especf = interpolate.interp1d(kfile, efile,'cubic')
self.kmin = kfile[0]
self.kmax = kfile[len(kfile) - 1]
def evaluate(self,k):
return self.especf(k)
class cbc_spectrum:
def __init__(self):
cbcspec = np.loadtxt('cbc_spectrum.txt')
kcbc=cbcspec[:,0]*100
ecbc=cbcspec[:,1]*1e-6
self.especf = interpolate.interp1d(kcbc, ecbc,'cubic')
self.kmin = kcbc[0]
self.kmax = kcbc[len(kcbc) - 1]
def evaluate(self,k):
return self.especf(k)
class vkp_spectrum:
def __init__(self,ke=40.0,nu=1.0e-5,urms=0.25):
self.nu = nu
self.urms = urms
self.ke = ke
self.kmin = 0
self.kmax = 1e6
def evaluate(self, k):
ke = self.ke
Nu = self.nu
urms = self.urms
# computed from input to satisfy homogeneous turbulence properties
Kappae = np.sqrt(5.0/12.0)*ke
L = 0.746834/Kappae #integral length scale
Alpha = 1.452762113;
# L = 0.05 # integral length scale
# Kappae = 0.746834/L
Epsilon = urms*urms*urms/L;
KappaEta = pow(Epsilon,0.25)*pow(Nu,-3.0/4.0);
r1 = k/Kappae
r2 = k/KappaEta
espec = Alpha*urms*urms/ Kappae * pow(r1,4)/pow(1.0 + r1*r1,17.0/6.0)*np.exp(-2.0*r2*r2)
return espec
class kcm_spectrum:
def __init__(self,station=0):
self.station_ = station
self.epst_ = [22.8, 9.13, 4.72, 3.41]
self.lt_ = [0.25, 0.288, 0.321, 0.332]
self.etat_= [0.11e-3, 0.14e-3, 0.16e-3, 0.18e-3]
self.ck_ = 1.613;
self.alfa_ = [0.39, 1.2, 4.0, 2.1, 0.522, 10.0, 12.58]
self.eps_ = self.epst_[station]
self.ckEpsTwo3rd_ = self.ck_*pow(self.eps_, 0.666666666667)
self.l_ = self.lt_[station]
self.eta_ = self.etat_[station]
def evaluate(self, k):
kl = k*self.l_;
keta = k*self.eta_;
term1 = kl / pow( (pow(kl, self.alfa_[1]) + self.alfa_[0]), 0.83333333333333)
term1 = pow(term1,5.66666666666667)
term2 = 0.5 + 0.31830988618379*np.arctan(self.alfa_[5]*np.log10(keta) + self.alfa_[6])
espec = self.ckEpsTwo3rd_*pow(k,-1.66666666666667)*term1 * np.exp(-self.alfa_[3]*keta) * (1.0 + self.alfa_[4]*term2)
return espec
class pq_spectrum:
# Implements the Passot-Pouquet Spectrum
# http://www.wseas.us/e-library/conferences/2011/Corfu/ASM/ASM-19.pdf
def __init__(self,ke=506, uavg=1.5):
self.uavg = uavg
self.ke = ke
self.kmin = 0
self.kmax = 1e6
def evaluate(self, k):
ke = self.ke
uavg = self.uavg
kke = k/ke
espec = 16.0*uavg*uavg/ke * np.sqrt(2.0/np.pi) * pow(kke,4) * np.exp(-2.0*(kke)*(kke))
return espec