Abschlussarbeit Strömungsakustik

Wenn du dir nicht sicher bist, in welchem der anderen Foren du die Frage stellen sollst, dann bist du hier im Forum für allgemeine Fragen sicher richtig.
Antworten
marteand
User
Beiträge: 7
Registriert: Montag 8. Juli 2013, 18:29

Hi @all

Ich habe ein Programm geschrieben mit dem mann die Schallausbreitung durch inhomogene Medien berechnen/beschreiben kann. Nun will ich dieses Programm noch weiter entwickeln. Die 3d Darstellung ist noch nicht perfekt und es ist noch nicht parametrisiert. Außerdem hätte ich gerne mehrere "ziele" (Mikrofone)
zu denen man den weg berechnen kann bzw. darstellen.
Würde mich über jede Hilfe sehr freuen. Das ganze ist im Rahmen einer Abschlussarbeit.
Ich habe das auch schon in einem anderen Forum gepostet, bin mir aber nicht sicher ob es da richtig ist, darum hier nochmal....

Programmtext:

from numpy import* #importiert alle Dateien vom Modul numpy; ermoeglicht Vektor- und Matrizenrechnung mit arrays

from scipy.integrate import odeint #importiert aus dem Modul Scipy die Funktion odeint zum loesen der DGL's
from scipy.optimize import fmin #importiert Funktionen von Scipy zur Minimierung des Abstands zwischen Mikrofon und Strahl

import matplotlib.pyplot as plt #importiert Dateien zum ploten
from mayavi import mlab #zur 3D-Darstellung der Loesung

from time import time #zur Ermittlung der Programmlaufzeit

c=343 #globale Konstante c Schallgeschwindigkeit [m/s]; die Schallgeschwindigkeit wird wie die Lufttemperatur und Luftfeuchte als konstant betrachtet
v0 = 50. #globale Konstante Stroemungsgeschwindigkeit [m/s]

n = 10 # Exponent der Geschwindigkeitsfunktion;globale Konstante Dimensionslos;es muss mindestens n = 2 sein bei groesser werden dem n wird das geschwindigkeitsprofil rechteckfoermiger
ko= 1 # Koeffizient der Geschwindigkeitsfunktion; globale Konstante, Dimensionslos;veraendert den Durchmesser des Geschwindigkeitsprofils
zx = 2
zy = 2
zz = 2
#n0 = array((nx, ny, nz))/sqrt(nx**2+ny**2+nz**2)
#x0 = array((sx, sy, sz))
ziel = array((zx, zy, zz))

n0 = array((1, 1, 1))/sqrt(3)
x0 = array((0, 0, 0.))
#ziel = array((1, 1, 1))
#s0 = n0 / (c+dot(v(x0)[0],n0))
#y0 = hstack((x0,s0))

t = linspace(0, 0.08, 2000)

def v(xx):
'Funktion zur Beschreibung Geschwindigkeitsprofils'
x,y,z = xx
v = array((v0*(1-tanh(ko*((y**2+z**2)**n))), 0., 0.)) #Geschwindigkeitsprofil
dv = array( ((0,0,0),((-n*v0*y**(n-1)*(1-tanh(ko*((y**2+z**2)**n)))**2),0,0),((-n*v0*z**(n-1)*(1-tanh(ko*((y**2+z**2)**n)))**2),0,0)) ) #Ableitung des Geschwindigkeitsprofils
return v,dv #

def f(y, t):
'Funktion zur Beschreibung der DGLs'
x = y[0:3]
s = y[3:6]
vv, dv = v(x)
sa = sqrt((s*s).sum()) # Betrag von s-Vektor
f0 = c*s/sa + vv
f1 = dot(s, dv)
return hstack((f0, f1))

def abstd(n):
'Funktion zur Minimierung des Abstands'
n0 = n/sqrt((n*n).sum())
s0 = n0 / (c+dot(v(x0)[0], n0))
y0 = hstack((x0, s0))

soln = odeint(f, y0, t)
ort = soln[:,0:3]
abstand = sqrt((abs(ziel-ort)**2).sum(1))
indmin = abstand.argmin()
tmin = t[indmin]
a = ort[indmin]
u = ort[indmin+1] - a
d = cross(ziel-a,u)
d = sqrt((d*d).sum())/sqrt((u*u).sum())
return abstand[indmin], tmin, d

def abstd2(n):
return abstd(n)[0]

ti = time()
#while d>0.001: ;versuch Hauptprogramm in schleife zu packen'

def prog():
'Funktion für Hauptprogramm'

ti = time()

s0 = n0 / (c+dot(v(x0)[0],n0))
y0 = hstack((x0,s0))
soln = odeint(f, y0, t)
ort = soln[:,0:3]
return prog

res = fmin(abstd2, n0,maxfun=300,full_output=0)

print time()-ti
print res
print abstd(res)
print sqrt(((x0-ziel)**2).sum())/c

s0 = n0 / (c+dot(v(x0)[0],n0))
y0 = hstack((x0,s0))
soln = odeint(f, y0, t)
ort = soln[:,0:3]
plt.plot(t,ort[:,0],label='x')
plt.plot(t,ort[:,1],label='y')
plt.plot(t,ort[:,2],label='z')
plt.legend()
plt.show()
pl=mlab.plot3d(ort[:,0],ort[:,1],ort[:,2],color=(1,0,0))
mlab.axes(pl)
mlab.show()


#for y in arange(-1,1,0.8):
#print v( array((0,y,0)) )
#y = arange(-2,2,0.01)
#vx = 1-tanh(ko*((y**2+z**2)**n
#vx = array((1-tanh(ko*(y**n)), 0., 1-tanh(ko*(y**n)))
#vxy = (-n*v0*y**(n-1)*(1-tanh(ko*((y**2+z**2)**n)))**2)
#plt.plot (y,vxy)
#plt.plot (y,vx)
#plt.show()
#pl=mlab.plot3d(vx[:,0],vx[:,1],vx[:,2],color=(0,1,0))
#mlab.axes(pl)
#mlab.show()
Antworten