Skip to content

Commit

Permalink
Added more docstrings
Browse files Browse the repository at this point in the history
  • Loading branch information
JBertholdt committed Feb 26, 2024
1 parent 5064b77 commit 59357cc
Show file tree
Hide file tree
Showing 4 changed files with 260 additions and 76 deletions.
81 changes: 68 additions & 13 deletions birdpressure/AnalBarber.py
Original file line number Diff line number Diff line change
Expand Up @@ -7,22 +7,63 @@
import matplotlib.pyplot as plt
import numpy as np
class AnalBarber(Timing):

def __init__(self, bird: Bird, scenario: ImpactScenario, grid: GridGenerator, eos: EOS, initial_pressure = 101325):
Timing.__init__(self, bird, scenario)
""" Class for Wilbeck's analytical bird strike solution
Parameters
----------
bird: class
impact_scenario: class
grid: class
eos: class
initial_pressure: int, default: 101325
ambient pressure default at sea-level
Attributes
----------
bird: class
impact_scenario: class
grid: class
eos: class
initial_pressure: int, default: 101325
ambient pressure default at sea-level
U: ndarray
superimposed velocity in x direction [m/s]
V: ndarray
superimposed velocity in y direction [m/s]
W: ndarray
superimposed velocity in z direction [m/s]
c_p: ndarray
array with pressure coefficient over impact area [-]
P: ndarray
array with pressures over impact area [Pa]
f_steady_av: float or ndarray
average force exerted during steady state phase [N]
P_steady_av: float or ndarray
average force exerted during steady state phase [Pa]
P_H: float or ndarray
peak pressure (hugonoit pressure) [MPa]
c_r: float or ndarray
release wave velocity [m/s]
t_b: float or ndarray
peak impact pressure duration [s]
t_c: float or ndarray
duration until shocked area dissipates [s]
"""

def __init__(self, bird: Bird, impact_scenario: ImpactScenario, grid: GridGenerator, eos: EOS, initial_pressure = 101325):
Timing.__init__(self, bird, impact_scenario)


self.initial_pressure = initial_pressure
self.bird = bird
self.scenario = scenario
self.impact_scenario = impact_scenario
self.grid = grid
self.eos = eos

self.min_axis = self.grid.min_axis
self.maj_axis = self.grid.maj_axis
#self.min_axis = self.grid.min_axis
#self.maj_axis = self.grid.maj_axis
self.U, self.V, self.W = self.get_UVW()
self.c_p = 1 - (self.V ** 2 + self.W ** 2) / (self.scenario.get_impact_velocity() ** 2)
self.P = self.c_p * 1 / 2 * self.bird.get_density() * self.scenario.get_impact_velocity() ** 2
self.c_p = 1 - (self.V ** 2 + self.W ** 2) / (self.impact_scenario.get_impact_velocity() ** 2)
self.P = self.c_p * 1 / 2 * self.bird.get_density() * self.impact_scenario.get_impact_velocity() ** 2
self.f_steady_av, self.P_steady_av = self.get_averages()
self.P_H = self.det_peak_P()
self.c_r = self.eos.get_c_r(self.P_H)
Expand All @@ -32,6 +73,7 @@ def __init__(self, bird: Bird, scenario: ImpactScenario, grid: GridGenerator, eo


def get_UVW(self):
"""Function to compute superimposd velocity in x,y and z direction."""
uvw = np.zeros([self.grid.ids.shape[0], 3])
for idx_1, id_cent in enumerate(self.grid.ids):

Expand All @@ -48,16 +90,19 @@ def get_UVW(self):
u, v, w, r_1, r_2, r_3, r_4 = self.get_velocity(x, y, z, eta_1, eta_2, zeta_1, zeta_2)
uvw[idx_1, :] += [u, v, w]

U = self.scenario.get_impact_velocity() * np.sin(self.scenario.get_impact_angle()) + self.scenario.get_impact_velocity() * np.sin(self.scenario.get_impact_angle()) / (
2 * np.pi) * uvw[:, 0]
V = self.scenario.get_impact_velocity() * np.cos(self.scenario.get_impact_angle()) + self.scenario.get_impact_velocity() * np.sin(self.scenario.get_impact_angle()) / (
2 * np.pi) * uvw[:, 1]
W = self.scenario.get_impact_velocity() * np.sin(self.scenario.get_impact_angle()) / (2 * np.pi) * uvw[:, 2]
U = self.impact_scenario.get_impact_velocity() * np.sin(self.impact_scenario.get_impact_angle()) + \
self.impact_scenario.get_impact_velocity() * np.sin(self.impact_scenario.get_impact_angle()) / (2 * np.pi) * uvw[:, 0]

V = self.impact_scenario.get_impact_velocity() * np.cos(self.impact_scenario.get_impact_angle()) + \
self.impact_scenario.get_impact_velocity() * np.sin(self.impact_scenario.get_impact_angle()) / (2 * np.pi) * uvw[:, 1]

W = self.impact_scenario.get_impact_velocity() * np.sin(self.impact_scenario.get_impact_angle()) / (2 * np.pi) * uvw[:, 2]

return U, V, W


def get_velocity(self, x, y, z, eta_1, eta_2, zeta_1, zeta_2):
"""Function to compute velocity field induced by uniformly distributed sources."""

r_1 = np.sqrt(x ** 2 + (y - eta_1) ** 2 + (z - zeta_1) ** 2)
r_2 = np.sqrt(x ** 2 + (y - eta_2) ** 2 + (z - zeta_1) ** 2)
Expand All @@ -80,6 +125,7 @@ def get_velocity(self, x, y, z, eta_1, eta_2, zeta_1, zeta_2):
return u, v, w, r_1, r_2, r_3, r_4

def mesh_results(self, result):
"""Function to save results in mesh format."""

result_cells = np.zeros_like(self.grid.mesh_X)

Expand All @@ -89,6 +135,7 @@ def mesh_results(self, result):
return result_cells

def plot_result(self, results, results_name: str, result_unit: str):
"""Function to plot results"""

result_cells = self.mesh_results(results)

Expand All @@ -104,6 +151,12 @@ def plot_result(self, results, results_name: str, result_unit: str):
plt.show()

def get_axis_data(self,result):
"""Function to plot results along major and minor axis
Notes
-----
In order to use this function a grid with uneven number of elements needs to be used.
"""
if (self.grid.n_elements%2) == 0:
print('Error: n_elements needs to be odd to generate axis data')
else:
Expand Down Expand Up @@ -142,6 +195,7 @@ def get_axis_data(self,result):
return maj_result, maj_ax_val, min_result, min_ax_val

def get_averages(self):
"""Function to compute steady state average values."""

cell_area = self.grid.resolution**2
impact_area = cell_area*len(self.P)
Expand All @@ -155,6 +209,7 @@ def get_averages(self):
return force, average_P

def show_averaged_force_history(self):
"""Function to plot force history"""
self.f_peak_av = self.P_H*self.grid.impact_area
self.f_hist, self.t_hist, self.impulse_tot = self.find_force_history(self.f_peak_av, self.f_steady_av,self.c_r)
self.plot_f_history(self.f_hist,self.t_hist)
Expand Down
78 changes: 60 additions & 18 deletions birdpressure/AnalWilbeck.py
Original file line number Diff line number Diff line change
Expand Up @@ -8,22 +8,38 @@
import numpy as np

class AnalWilbeck(Timing):

def __init__(self, bird: Bird, impact_scenario: ImpactScenario, grid: GridGenerator, eos: EOS, initial_pressure = 101325):
"""
""" Class for Wilbeck's analytical bird strike solution
Parameters
----------
bird:
impact_scenario
grid
eos
initial_pressure
"""
bird: class
impact_scenario: class
grid: class
eos: class
initial_pressure: int, default: 101325
ambient pressure default at sea-level
Attributes
----------
initial_pressure: int, default: 101325
ambient pressure default at sea-level
eos: class
bird: class
impact_scenario: class
grid: class
u_s: float or ndarray
shock wave velocity [m/s]
P_H: float or ndarray
peak (hugonoit) pressure [N/m^2]
c_r: float or ndarray
release wave velocity [m/s]
P_s: float ndarray
steady state stagnation point pressure [N/m^2]
"""

def __init__(self, bird: Bird, impact_scenario: ImpactScenario, grid: GridGenerator, eos: EOS, initial_pressure = 101325):
Timing.__init__(self, bird, impact_scenario)

self.initial_pressure = initial_pressure

self.eos = eos
self.bird = bird
self.impact_scenario = impact_scenario
Expand All @@ -32,16 +48,12 @@ def __init__(self, bird: Bird, impact_scenario: ImpactScenario, grid: GridGenera
self.u_s = self.find_u_s()
self.P_H = self.det_peak_P()
self.c_r = self.eos.get_c_r(self.P_H)

self.P_s = self.get_P_s()

# def det_peak_P(self, use_us = True):
# if use_us:
# return self.bird.get_density()*self.u_s*self.impact_scenario.get_impact_velocity()
# return self.bird.get_density()*self.bird.get_sound_speed()*self.impact_scenario.get_impact_velocity()
def get_averages(self, circular_cross = True):
"""Function to compute average pressure and force during steady state"""

def get_averages(self, Circular_cross = True):
if Circular_cross:
if circular_cross:
impact_A = np.pi * (self.bird.diameter/2)**2
else:
impact_A = np.pi * self.grid.min_axis * self.grid.maj_axis
Expand All @@ -52,17 +64,45 @@ def get_averages(self, Circular_cross = True):
return P_avg, F_avg

def get_P_s(self):
"""Funtion to compute steady state stagnation point pressure"""
return 1/2 * self.bird.get_density() * self.impact_scenario.get_impact_velocity()**2 * np.sin(self.impact_scenario.get_impact_angle())
def find_leach(self,r,xi_2 =2.58 ):
"""Function for Leach & Walker steady state pressure distribution.
Parameters
----------
r: float or ndarray
distance from center [m]
xi_2: float, default: 2.58
material constant [-]
Returns
-------
P: float or ndarray
Pressure distribution [Pa]
"""
P = self.P_s * (1 - 3 * (r / (xi_2 * (self.bird.get_diameter()/2))) ** 2 + 2 * (r / (xi_2 * (self.bird.get_diameter()/2))) ** 3)
return P

def find_banks(self,r,xi_1 = 0.5):
"""Function for Banks & Chandrasekhara steady state pressure distribution.
Parameters
----------
r: float or ndarray
distance from center [m]
xi_1: float, default: 0.5
material constant [-]
Returns
-------
P: float or ndarray
Pressure distribution [Pa]
"""
P = self.P_s * np.e ** (-xi_1 * (r / (self.bird.get_diameter()/2)) ** 2)
return P


def show_leach(self):
"""Funtion to plot Leach & Walker steady state pressure distribution """
r = np.linspace(0,4*self.bird.get_length()/2,50)
P_dist = self.find_leach(r)

Expand All @@ -84,6 +124,7 @@ def show_leach(self):
self.P_leach = P_dist

def show_banks(self):
"""Funtion to plot Banks & Chandrasekhara steady state pressure distribution """
r = np.linspace(0,4*self.bird.get_length()/2,50)
P_dist = self.find_banks(r)

Expand All @@ -106,6 +147,7 @@ def show_banks(self):


def show_averaged_force_history(self):
"""Function to plot force history"""
self.p_steady_av, self.f_steady_av = self.get_averages()
self.f_peak_av = self.P_H*self.grid.impact_area
self.f_hist, self.t_hist, self.impulse_tot = self.find_force_history(self.f_peak_av, self.f_steady_av, self.c_r)
Expand Down
60 changes: 50 additions & 10 deletions birdpressure/EOS.py
Original file line number Diff line number Diff line change
Expand Up @@ -5,23 +5,48 @@


class EOS():
""" Equation of state function to determine c_r
Parameters
----------
bird: class
name: str
name of eos
Attributes
----------
bird: class
name: str
name of eos
name_lst: list
list of available eos names
rho: sympy.core.symbol.Symbol
density sympu symbol
P_0: int
reference pressure Murnaghan [Pa]
B: int
Murnaghan material constant [Pa]
gamma: float
Murnaghan material constant [-]
P_func: sympy.core.mul.Mul
sympy function of chosen eos
K: int
linear eos bulk modulus [Pa]
mu: sympy.core.mul.Mul
change of density during impact
k : int
polynomial eos experimental constant
"""

def __init__(self, bird: Bird, name: str):
"""
Equation of state Class uses sympy to use different equations
:param bird: Bird Class
:param name: Name EOS you want to use ['Murnaghan', 'Linear_EOS', 'Polynomial_EOS']
:param peak_pressur: Set pressure around which c_r is determined
"""


self.bird = bird
self.name = name
self.name_lst = ['Murnaghan', 'Linear_EOS', 'Polynomial_EOS']



self.rho = sy.symbols('rho', real = True, positive=True)

if self.name == "Murnaghan":
self.P_0 = 0
self.B = 128e6
Expand Down Expand Up @@ -52,23 +77,38 @@ def __init__(self, bird: Bird, name: str):


def get_c_r(self,peak_pressure):
rho = self.get_density(peak_pressure)
"""Compute release wave velocity by determining eos slope at peak pressure
Parameters
----------
peak_pressure: float
hugonoit pressure [Pa]
Returns
-------
c_r: float
release wave velocity [m/s]
"""
rho_peak = self.get_density(peak_pressure)

c_r_squared = sy.diff(self.P_func)
c_r = np.sqrt(float(c_r_squared.subs(self.rho,rho)))
c_r = np.sqrt(float(c_r_squared.subs(self.rho,rho_peak)))
return c_r

def get_density(self,P):
"""Compute density from pressure"""
eq = sy.Eq(self.P_func,P)

return sy.nsolve(eq,self.rho, self.bird.get_density())

def get_pressure(self,rho):
"""Compute pressure from density"""
return self.P_func.subs(self.rho,rho)



def eos_list(self):
"""Print string of eos options"""
print('Current options of EOS:')
for eos in self.name_lst:
print(eos)
Loading

0 comments on commit 59357cc

Please sign in to comment.