To receive notifications about scheduled maintenance, please subscribe to the mailing-list gitlab-operations@sympa.ethz.ch. You can subscribe to the mailing-list at https://sympa.ethz.ch

Commit 681815f5 authored by mastc's avatar mastc 🎯
Browse files

Revert "Merge branch 'patch/yosioka' into 'patch/Plotting'"

This reverts commit 96a41b84
parent 96a41b84
Pipeline #88972 passed with stages
in 4 minutes and 46 seconds
Yosioka (1955)
==============
`Link to paper
<https://www.ingentaconnect.com/content/dav/aaua/1955/00000005/00000003/art00004>`_
.. inheritance-diagram::
gorkov.solutions.Yosioka1955.arf
:top-classes: gorkov.solutions.base_arf.BaseARF
:parts: 1
.. currentmodule:: gorkov.solutions.Yosioka1955
.. autosummary::
:toctree: generated/
:template: class_property.rst
:nosignatures:
~baseclass.BaseYosioka
~arf.YosiokaARF
......@@ -73,4 +73,3 @@ from gorkov.core.functions import (
import gorkov.solutions.King1934 as King1934
import gorkov.solutions.Gorkov1962 as Gorkov1962
import gorkov.solutions.Yosioka1955 as Yosioka1955
......@@ -42,16 +42,12 @@ class BackgroundField:
self._position = PassiveVariable(position, "particle position d")
# Dependent variables
self._E_bar = ActiveVariable(self._compute_E_bar,
"mean total energy")
self._A = ActiveVariable(self._compute_A,
"velocity potential amplitude A")
self._A_in = ActiveVariable(self._reset_A_in,
"velocity potential amplitude A_n")
# Dependencies
self._E_bar.is_computed_by(self.fluid._rho_f, self.fluid._k_f)
self._A.is_computed_by(self.fluid.frequency._omega,
self.fluid._rho_f,
self._p_0)
......@@ -115,14 +111,6 @@ class BackgroundField:
# Getters for Dependent Variables
# -------------------------------------------------------------------------
@property
def E_bar(self) -> float:
"""mean total energy
:rtype: float
"""
return self._E_bar.value
@property
def A(self) -> complex:
"""Amplitude of the pressure wave
......@@ -180,9 +168,6 @@ class BackgroundField:
# Methods
# -------------------------------------------------------------------------
def _compute_E_bar(self) -> float:
return self.rho_f * self.k_f ** 2
def _compute_A_in(self, N: int) -> None:
N_old = len(self._A_in.value)
......
......@@ -26,7 +26,7 @@ class InviscidFluid(BaseFrequencyComposite):
super().__init__(frequency)
# Independent Variables
self._rho_f = PassiveVariable(rho, "density rho of the fluid")
self._rho_f = PassiveVariable(rho, "density rho")
self._c_f = PassiveVariable(c, "speed of sound c")
# Dependent variables
......@@ -67,7 +67,7 @@ class InviscidFluid(BaseFrequencyComposite):
def rho_f(self) -> float:
"""Density of the fluid.
:getter: returns the value for the density of the fluid
:getter: returns the value for the density
:rtype: float
:setter: automatically invokes :func:`~gorkov.core.BaseVariable.notify`
:type: float
......
......@@ -4,6 +4,7 @@ import numpy as np
from scipy.integrate import romb, quad
from scipy.special import hankel1, lpmv, jv, yv, factorial
from gorkov import log
# Often used NumPy and SciPy Functions and Classes
......
......@@ -26,7 +26,7 @@ class RigidSolid(BaseFrequencyComposite):
# Logging
# Independent variables
self._rho_s = PassiveVariable(rho, "density rho of the solid")
self._rho_s = PassiveVariable(rho, "density rho")
if type(self) is RigidSolid:
log.debug(f'Creating {self}')
......@@ -42,7 +42,7 @@ class RigidSolid(BaseFrequencyComposite):
def rho_s(self) -> float:
"""Density of the solid.
:getter: returns the value for the density of the solid
:getter: returns the value for the density
:rtype: float
:setter: automatically invokes
:func:`gorkov.core.variable.BaseVariable.notify`
......
from gorkov.solutions.Yosioka1955.baseclass import BaseYosioka
from gorkov.solutions.Yosioka1955.scatteringfield import ScatteringField
from gorkov.solutions.Yosioka1955.arf import ARF
from __future__ import annotations
from typing import Union, Optional
from gorkov import (
log, pi, sin, cos, conj,
BaseSphereFrequencyComposite,
Sphere, Frequency, integrate)
from gorkov.solutions.base_arf import BaseARF
from .baseclass import BaseYosioka
from gorkov.solutions.Yosioka1955.scatteringfield import ScatteringField
class ARF(BaseYosioka, BaseARF, BaseSphereFrequencyComposite):
"""ARF according to Yosioka's theory of 1955
.. warning::
This model is based on the following assumptions:
- Inviscid fluid
- Spheres of finite compressibility
- plane progressive / standing wave field
"""
def __init__(self, frequency: Union[Frequency, float, int],
radius: Union[Sphere, float, int],
rho_s: float, rho_f: float, c_f: float, c_s: float,
p_0: float, position: Optional[float],
wave_type: str = "standing",
particle_type: str = "bubble", n_max: int = None) -> None:
""" Initializes the class and sets all parameter needed for the
computation.
:param frequency: Frequency [Hz]
:type frequency: Union[Frequency, float, int]
:param radius: Radius of the solid [m]
:type radius: Union[Sphere, float, int]
:param rho_s: Density of the fluid-like solid [kg/m^3]
:type rho_s: float
:param rho_f: Density of the fluid [kg/m^3]
:type rho_f: float
:param c_f: Speed of sound of the fluid [m/s]
:type c_f: float
:param p_0: Pressure amplitude of the field [Pa]
:type p_0: float
:param position: Position within the standing wave field [m]
:type position: float
:param wave_type: either standing or progressive wavefield
:type wave_type: str
:rtype: None
"""
# init of parent class
BaseSphereFrequencyComposite.__init__(self, frequency, radius)
BaseYosioka.__init__(self,
frequency=frequency,
radius=radius,
rho_s=rho_s,
rho_f=rho_f,
c_f=c_f,
c_s=c_s,
p_0=p_0,
position=position,
wave_type=wave_type,
particle_type=particle_type,
n_max=n_max)
self.scatter = ScatteringField(frequency=frequency, radius=radius,
rho_s=rho_s, rho_f=rho_f,
c_f=c_f, c_s=c_s,
p_0=p_0, wave_type=wave_type,
particle_type=particle_type,
position=position,
n_max=n_max)
def acoustic_radiation_force(self) -> float:
"""Computes the ARF and returns the force in Newton. Depending on the
user-input assumptions are made for either small spheres or bubbles.
:rtype: float
"""
if self.particle_type == "analytical":
out = self._compute_analytical_arf()
elif self.particle_type == "numerical":
out = self._compute_numerical_arf()
elif self.particle_type == "special":
out = self._compute_special_arf()
elif self.particle_type == "bubble":
out = self._compute_bubble_arf()
elif self.particle_type == "sphere":
out = self._compute_sphere_arf()
else:
print(self.particle_type)
raise NotImplementedError('Solution is no implemented (yet).')
surface = 4 * pi * self.R_0 ** 2
return out * surface
def _compute_analytical_arf(self) -> float:
"""Computes general analytical solution for the ARF according to
eq. (44)
:rtype: float
"""
if self.n_max is None:
self.n_max = 10
n = 0
out = [0 for _ in range(5)]
while n < self.n_max:
# (45)
s1 = (self.K_n(n) * conj(self.K_n(n + 1))).real
s1 *= n + 1
out[0] += s1
# (46)
s2 = (self.M_n(n) * conj(self.M_n(n + 1))).real
s2 *= n * (n + 1) * (n + 2)
out[1] += s2
# (47)
s3 = (self.K_n(n) * conj(self.M_n(n + 1))).real
s3 *= (n + 1) * (n + 2)
s4 = (self.M_n(n) * conj(self.K_n(n + 1))).real
s4 *= n * (n + 1)
out[2] += s3
out[3] += s4
# (48)
s5 = (self.M_n(n) * conj(self.M_n(n + 1))).real
s5 *= n + 1
out[4] += s5
n += 1
out[0] *= -2 * pi * self.R_0 ** 2 * self.rho_f
out[1] *= 2 * pi * self.xlambda ** 2 * self.rho_f
out[2] *= -2 * pi * self.R_0 * self.xlambda * self.rho_f
out[3] *= 2 * pi * self.R_0 * self.xlambda * self.rho_f
out[4] *= (-2 * pi * (self.k_f ** 2) *
(self.R_0 ** 2) * (self.xlambda ** 2) * self.rho_f)
print(out)
# (44)
return sum(out)
def _compute_numerical_arf(self) -> float:
"""Computes general numerical solution for the ARF according to
eq. (11)
:rtype: float
"""
p_uv = self.compute_p_numerical("uv")
p_ut = self.compute_p_numerical("ut")
p_uv_ut = self.compute_p_numerical("uv_ut")
p_phi = self.compute_p_numerical("phi")
return sum([p_uv, p_ut, p_uv_ut, p_phi])
def compute_p_numerical(self, dir: str) -> float:
out = self.scatter.omega / (2 * pi)
if dir == "uv":
out *= integrate(
self.time_averaging_uv, 0.0,
(2.0 * pi / self.scatter.omega))
elif dir == "ut":
out *= integrate(
self.time_averaging_ut, 0.0,
(2.0 * pi / self.scatter.omega))
elif dir == "uv_ut":
out *= integrate(
self.time_averaging_uv_ut, 0.0,
(2.0 * pi / self.scatter.omega))
elif dir == "phi":
out *= integrate(
self.time_averaging_phi, 0.0,
(2.0 * pi / self.scatter.omega))
return out
def time_averaging_uv(self, t):
out = []
for dt in t:
self.scatter.t = dt
out.append(-pi * self.R_0 ** 2 * self.rho_f *
integrate(self.integrand_uv, 0, pi))
return out
def integrand_uv(self, x):
out = self.scatter.radial_acoustic_fluid_velocity(
r=self.scatter.R_0, theta=x, t=self.scatter.t) ** 2
out *= sin(x) * cos(x)
return out
def time_averaging_ut(self, t):
out = []
for dt in t:
self.scatter.t = dt
out.append(pi * self.scatter.rho_f *
integrate(self.integrand_ut, 0, pi))
return out
def integrand_ut(self, x):
out = self.scatter.tangential_acoustic_fluid_velocity(
r=self.scatter.R_0, theta=x, t=self.scatter.t)
out *= self.scatter.R_0
out **= 2
out *= sin(x) * cos(x)
return out
def time_averaging_uv_ut(self, t):
out = []
for dt in t:
self.scatter.t = dt
out.append(
2 * pi * self.scatter.R_0 *
self.scatter.xlambda * self.scatter.rho_f *
integrate(self.integrand_uv_ut, 0, pi))
return out
def integrand_uv_ut(self, x):
out = self.scatter.radial_acoustic_fluid_velocity(
r=self.scatter.R_0, theta=x, t=self.scatter.t)
out *= self.scatter.tangential_acoustic_fluid_velocity(
r=self.scatter.R_0, theta=x, t=self.scatter.t)
out *= self.scatter.R_0
out *= sin(x) ** 2
return out
def time_averaging_phi(self, t):
out = []
for dt in t:
self.scatter.t = dt
out.append(
(-pi * self.scatter.R_0 ** 2
* self.scatter.rho_f / self.scatter.c_f ** 2) *
integrate(self.integrand_uv, 0, pi))
return out
def integrand_phi(self, x):
out = self.scatter.d_Phi1_dt(
r=self.scatter.R_0, theta=x, t=self.scatter.t) ** 2
out *= sin(x) * cos(x)
return out
def _compute_special_arf(self) -> float:
if self.wave_type == 'traveling':
return self._traveling_wave_solution_special()
else:
raise NotImplementedError('Solution is no implemented (yet).')
def _traveling_wave_solution_special(self) -> float:
"""based on eq. (56)
:rtype: float
"""
out = 0
n = 0
while True:
num = ((n + 1) *
((self._D(1, n) * self._D(2, n + 1)
- self._D(2, n) * self._D(1, n + 1)))
** 2)
denom = ((self._D(1, n) ** 2 + self._D(2, n) ** 2) *
(self._D(1, n + 1) ** 2 + self._D(2, n + 1) ** 2))
out_prev = out
try:
out += num / denom
except ZeroDivisionError:
print("Early Termination due to zero division")
break
# condition when to terminate
# (change is less than 0.1% and at least 1000 iterations)
try:
error = abs((out - out_prev) / out_prev)
# at the beginning
except ZeroDivisionError:
error = 1
if n >= 100 or error < 0.001:
# print("Termination after {} iterations".format(n))
break
n += 1
out *= 2 * pi * self.rho_f
return out
def _compute_sphere_arf(self):
order_sphere = 1.0
mag = 10
conditions = [self.xlambda < (1.0 / mag) * order_sphere,
self.xlambda > (1.0 * mag) * order_sphere,
(self.k_s * self.R_0) ** 2 > 1.0,
(self.k_f * self.R_0) ** 2 > 1.0]
for i, c in enumerate(conditions):
if c:
text = 'Criteria for small spheres might not be fulfilled.'
text += 'Consider using another model approach'
log.info(text)
break
if self.wave_type == 'standing':
return self._standing_wave_solution_sphere()
else:
return self._traveling_wave_solution_sphere()
def _traveling_wave_solution_sphere(self) -> float:
"""based on eq. (58)
:rtype: float
"""
out = (2 * pi * self.rho_f *
(self.k_f * self.R_0) ** 6 * self._F.value)
return out
def _standing_wave_solution_sphere(self) -> float:
"""based on eq. (62)
:rtype: float
"""
out = (4 * pi * self.R_0 ** 3 * self.k_f * self.field.E_bar
* sin(2 * self.k_f * self.position) * self._F.value)
return out
def _compute_bubble_arf(self):
order_bubble = (self.k_f * self.R_0) ** 2
mag = 10
conditions = [self.xlambda < (1.0 / mag) * order_bubble,
self.xlambda > (1.0 * mag) * order_bubble,
(self.k_s * self.R_0) ** 2 > 1.0,
(self.k_f * self.R_0) ** 2 > 1.0]
for i, c in enumerate(conditions):
if c:
text = 'Criteria for small bubbles might not be fulfilled.'
text += 'Consider using another model approach'
log.info(text)
break
if self.wave_type == 'standing':
return self._standing_wave_solution_bubble()
else:
return self._traveling_wave_solution_bubble()
def _traveling_wave_solution_bubble(self) -> float:
"""based on eq. (67)
:rtype: float
"""
out = (2 * pi * self.rho_f * self.sigma ** 2
* (self.k_f * self.R_0) ** 6)
out /= (self.sigma ** 2 * (self.k_f * self.R_0) ** 6
+ (3 * self.xlambda - (self.k_f * self.R_0) ** 2) ** 2)
return out
def _standing_wave_solution_bubble(self) -> float:
"""based on eq. (73)
:rtype: float
"""
out = -(4 * pi * self.rho_f *
sin(2 * self.k_f * self.position) * self._F.value)
return out
if __name__ == '__main__':
pass
from __future__ import annotations
from typing import Union
from gorkov import (
log, pi, full_range,
BaseSphereFrequencyComposite,
Sphere, Frequency,
BackgroundField,
InviscidFluid,
ActiveVariable, PassiveVariable,
StringFormatter as SF,
SpecialFunctions as sp
)
class BaseYosioka(BaseSphereFrequencyComposite):
"""Baseclass for Yosioka's (1955) ARF solution
"""
def __init__(self, frequency: Union[Frequency, float, int],
radius: Union[Sphere, float, int],
rho_s: float,
rho_f: float, c_f: float, c_s: float,
p_0: float, position: float,
particle_type: str = "bubble",
wave_type: str = "standing", n_max=None) -> None:
"""Initializes the class and sets all parameter needed for the
computation.
:param frequency: Frequency [Hz]
:type frequency: Union[Frequency, float, int]k
:param radius: Radius of the solid [m]
:type radius: Union[Sphere, float, int]
:param rho_s: Density of the fluid-like solid [kg/m^3]
:type rho_s: float
:param rho_f: Density of the fluid [kg/m^3]
:type rho_f: float
:param c_f: Soundspeed of the fluid [m/s]
:type c_f: float
:param p_0: Pressure amplitude of the field [Pa]
:type p_0: float
:param position: Position within the standing wave field [m]
:type position: Optional[float]
:param wave_type: either standing or progressive wavefield
:type wave_type: str
:rtype: None
"""
# init of parent class
BaseSphereFrequencyComposite.__init__(self, frequency, radius)
# Initialize Components
self.solid = InviscidFluid(self.frequency, rho_s, c_s)
self.fluid = InviscidFluid(self.frequency, rho_f, c_f)
self.field = BackgroundField(self.fluid, p_0, wave_type, position)
# Independent variables
self.n_max = n_max
self._particle_type = PassiveVariable(
particle_type,
"Type of the particle")
# Independent variables:
self._omega = PassiveVariable(2 * pi * frequency, "Angular Frequency")
self._R_0 = PassiveVariable(radius, "Particle-Radius")
# Dependent variables
self._xlambda = ActiveVariable(self._compute_xlambda,
"ratio of densities")
self._sigma = ActiveVariable(self._compute_sigma,
"ratio of speed of sounds")
self._F = ActiveVariable(self._compute_F,
"density-compressibility factor")
self._A_n = ActiveVariable(self._reset_coeff,
"Coefficient A_n"
" for scattered potential")
self._B_n = ActiveVariable(self._reset_coeff,
"Coefficient B_n"
" for potential inside particle")
self._M_n = ActiveVariable(self._reset_coeff,
"Coefficient M_n")
self._K_n = ActiveVariable(self._reset_coeff,
"Coefficient K_n")
# Dependencies
self._xlambda.is_computed_by(
self.solid._rho_f, self.fluid._rho_f
)
self._sigma.is_computed_by(
self.fluid._k_f, self.solid._k_f
)
self._F.is_computed_by(
self._sigma,
self._xlambda,
self.field._wave_type,
self._particle_type
)
self._A_n.is_computed_by(
self._xlambda, self.fluid._k_f,
self.solid._k_f, self._R_0
)
self._B_n.is_computed_by(
self._xlambda, self.fluid._k_f,
self.solid._k_f, self._R_0
)
self._M_n.is_computed_by(self._omega, self.solid._k_f,
self._R_0, self._B_n
)
self._K_n.is_computed_by(self._omega, self.solid._k_f,
self._R_0, self._B_n
)
log.info(str(self))
log.debug(repr(self))
def __repr__(self):
return (
f'YosiokaARF({self.f}, {self.R_0}, {self.rho_s}, '
f'{self.rho_f}, {self.c_f}, '
f'{self.p_0}, {self.position}'
)