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 3af45b6e authored by Cyrill Mast's avatar Cyrill Mast
Browse files

use E_ac and I_ac

add plots
minor changes
parent c5ac6b96
......@@ -76,7 +76,9 @@ class ARF(ScatteringField, BaseARF):
self._sigma,
self._xlambda,
self.field._wave_type,
self._particle_type)
self._particle_type,
self.solid._k_f,
self._R_0)
self._K_n.is_computed_by(
self.solid._k_f,
......@@ -157,8 +159,18 @@ class ARF(ScatteringField, BaseARF):
out[4] *= (-2 * pi * (self.k_f ** 2) *
(self.R_0 ** 2) * (self.xlambda ** 2) * self.rho_f)
# (44)
return sum(out)
# (eq. 44)
out = sum(out)
"""
Yosioka has a different definition for the standing wave,
therefore we need to introduce an opposite sign in the case
of the standing wave field
"""
if self.wave_type == "standing":
out *= -1
return out
def _compute_special_arf(self) -> float:
if self.wave_type == 'traveling':
......@@ -219,27 +231,27 @@ class ARF(ScatteringField, BaseARF):
out = self._standing_wave_solution_sphere()
else:
out = self._traveling_wave_solution_sphere()
return out * (abs(self.field.A / 2) ** 2).real
return out.real
def _traveling_wave_solution_sphere(self) -> float:
"""
based on eq. (58)
based on eq. (59)
:rtype: float
"""
out = (2 * pi * self.rho_f *
(self.k_f * self.R_0) ** 6
out = (pi * self.R_0 ** 2 * 4 * (self.k_f * self.R_0) ** 4
* self.I_ac / self.c_f
* self.density_compressibility_factor)
return out
def _standing_wave_solution_sphere(self) -> float:
"""
based on eq. (61)
based on eq. (62)
:rtype: float
"""
out = (4 * pi * self.rho_f *
(self.k_f * self.R_0) ** 3
out = (pi * self.R_0 ** 2 * 4 * self.k_f * self.R_0
* self.E_ac
* sin(2 * self.k_f * self.position)
* self.density_compressibility_factor)
return out
......@@ -262,30 +274,30 @@ class ARF(ScatteringField, BaseARF):
out = self._standing_wave_solution_bubble()
else:
out = self._traveling_wave_solution_bubble()
return out * (abs(self.field.A / 2) ** 2).real
return out.real
def _traveling_wave_solution_bubble(self) -> float:
"""
based on eq. (67)
based on eq. (68)
:rtype: float
"""
out = (2 * pi * self.rho_f * self.sigma ** 2
* (self.k_s * self.R_0) ** 6)
out /= (self.sigma ** 2 * (self.k_s * self.R_0) ** 6
+ (3 * self.xlambda -
(self.k_s * self.R_0) ** 2) ** 2)
out = (4 * pi * self.R_0 ** 2
* (self.k_s * self.R_0) ** 4
* self.I_ac)
out /= (self.c_f * (self.sigma ** 2 * (self.k_s * self.R_0) ** 6
+ (3 * self.xlambda - (self.k_s * self.R_0) ** 2) ** 2))
return out
def _standing_wave_solution_bubble(self) -> float:
"""
based on eq. (73)
based on eq. (74)
:rtype: float
"""
out = -(4 * pi * self.rho_f *
sin(2 * self.k_f * self.position)
* self.density_compressibility_factor)
out = -4 * pi / self.k_f ** 2
out *= self.E_ac
out *= sin(2 * self.k_f * self.position)
out *= self.density_compressibility_factor
return out
# -------------------------------------------------------------------------
......@@ -347,10 +359,10 @@ class ARF(ScatteringField, BaseARF):
F -= 1 / (3 * self.xlambda * self.sigma ** 2)
# (75)
elif self.particle_type == "bubble":
F = (self.sigma * self.k_s * self.R_0 *
(3 * self.xlambda - (self.k_s * self.R_0) ** 2))
temp1 = 3 * self.xlambda - (self.k_s * self.R_0) ** 2
F = self.sigma * self.k_s * self.R_0 * temp1
F /= (self.sigma ** 2 * (self.k_s * self.R_0) ** 6
+ (3 * self.xlambda - (self.k_s * self.R_0) ** 2) ** 2)
+ temp1 ** 2)
else:
raise NotImplementedError
return F
......
......@@ -290,6 +290,7 @@ class BaseYosioka(BaseSphereFrequencyComposite):
# -------------------------------------------------------------------------
# Wrappers for Dependent Solid Attributes
# -------------------------------------------------------------------------
@property
def k_s(self) -> float:
"""
......@@ -299,6 +300,26 @@ class BaseYosioka(BaseSphereFrequencyComposite):
"""
return self.solid.k_f
# -------------------------------------------------------------------------
# Wrappers for Dependent Field Attributes
# -------------------------------------------------------------------------
@property
def I_ac(self) -> float:
"""
Wraps to
:attr:`gorkov.core.backgroundfields.BackgroundField.acoustic_intensity`
"""
return self.field.I_ac
@property
def E_ac(self) -> float:
"""
Wraps to
:attr:`gorkov.core.backgroundfields.BackgroundField.E_ac`
"""
return self.field.E_ac
# ----------------------------------
# Wrappers for independent variables
# ----------------------------------
......
......@@ -81,14 +81,6 @@ class ScatteringField(BaseYosioka, BaseScattering):
self.solid._k_f, self._R_0, self.field._A_in
)
def __repr__(self):
return (
f'ScatteringField({self.f}, {self.R_0}, {self.rho_s}, '
f'{self.rho_f}), {self.c_f}. {self.c_s}, {self.p_0}, '
f'{self.wave_type}, {self.particle_type}, '
f'{self.position}'
)
# -------------------------------------------------------------------------
# Velocity Potentials
# -------------------------------------------------------------------------
......
......@@ -4,15 +4,24 @@ from gorkov import (
pi,
Yosioka1955, King1934, Gorkov1962)
from plotting.plot_Acoustofluidics import PlotAcoustofluidics
import matplotlib.pyplot as plt
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import rc
rc('font', **{'family': 'sans-serif', 'sans-serif': ['Helvetica']})
# for Palatino and other serif fonts use:
# rc('font',**{'family':'serif','serif':['Palatino']})
rc('text', usetex=True)
plt.rc('text', usetex=True)
plt.rc('font', family='serif')
class TestYosiokaARF(BaseTest):
def setUp(self) -> None:
self.line_palette = ["-", "-.", ":", "--", "dotted", "dashdot"]
# Geometry
self.R_0 = 50e-6
......@@ -20,24 +29,16 @@ class TestYosiokaARF(BaseTest):
self.f = 1e5
# Solid parameters
self.E_s = 0.32e10
self.nu_s = 0.35
self.rho_s = 1.05e3
self.c_s = 1.5e3
# Fluid parameters
self.rho_f = 1e3
self.c_f = 1.5e3
self.eta_f = 1e-3
self.eta_p = 9e-3
self.zeta_f = 0
self.zeta_p = 0
self.lambda_M = 1e-6
# Incident field
self.p0 = 10e3
self.wave_type = 'traveling'
self.inf_type = 'radius'
self.pos = 1e-4
n_max = 10
......@@ -113,10 +114,10 @@ class TestYosiokaARF(BaseTest):
self.R_0 = 50e-6
self.f = 1e5
self.rho_s = 1.225
self.c_s = 1.5e3
self.rho_f = 1e3
self.c_s = 343
self.rho_f = 997
self.c_f = 1.5e3
self.p0 = 10e3
self.p0 = 1e4
self.pos = 1e-4
order_bubble = (self.compute_k_f() * self.R_0) ** 2
......@@ -163,11 +164,13 @@ class TestYosiokaARF(BaseTest):
for key in self.cls_dict:
self.cls_dict[key].f = self.f
self.cls_dict[key].c_f = self.c_f
self.cls_dict[key].c_s = self.c_s
self.cls_dict[key].rho_f = self.rho_f
self.cls_dict[key].rho_s = self.rho_s
self.cls_dict[key].p_0 = self.p0
self.cls_dict[key].position = self.pos
self.cls_dict[key].R_0 = self.R_0
self.cls_dict[key].wave_type = self.wave_type
def test_xlambda(self) -> None:
for _ in range(10):
......@@ -189,6 +192,16 @@ class TestYosiokaARF(BaseTest):
except AttributeError:
pass
def test_F(self) -> None:
self.cls_bubble.wave_type = "traveling"
print(self.cls_bubble.density_compressibility_factor)
self.cls_bubble.wave_type = "standing"
print(self.cls_bubble.density_compressibility_factor)
self.cls_sphere.wave_type = "traveling"
print(self.cls_sphere.density_compressibility_factor)
self.cls_sphere.wave_type = "standing"
print(self.cls_sphere.density_compressibility_factor)
def test_compute_coefficients(self) -> None:
"""
check if coefficients are computed correctly
......@@ -215,30 +228,30 @@ class TestYosiokaARF(BaseTest):
def test_plot_coefficients(self) -> None:
n = np.arange(10)
A_n = []
B_n = []
K_n = []
M_n = []
self.cls_general.wave_type = "standing"
self.cls_general.position = 0
coeffs = [K_n, M_n]
self.cls_general.wave_type = "traveling"
self.cls_general.position = 1e-4
for N in n:
A_n.append(self.cls_general.A_n(N).real)
B_n.append(self.cls_general.B_n(N).real)
K_n.append(self.cls_general.K_n(N).real)
M_n.append(self.cls_general.M_n(N).real)
A_n /= A_n[0]
B_n /= B_n[0]
K_n /= K_n[0]
M_n /= M_n[0]
K_n.append(abs(self.cls_general.K_n(N)))
M_n.append(abs(self.cls_general.M_n(N)))
for c in coeffs:
c /= c[0]
fig, ax = plt.subplots()
ax.set_xlabel("n")
ax.set_ylabel("Imaginary part of Coefficient")
ax.set_title("Standing wave field at position 0")
ax.semilogy(n, abs(A_n))
ax.semilogy(n, abs(B_n))
ax.semilogy(n, abs(K_n))
ax.semilogy(n, abs(M_n))
ax.legend(["A_n", "B_n", "K_n", "M_n"])
ax.set_xlabel(r"Order $n$")
ax.set_ylabel(r"Absolute value of the complex coefficients")
ax.set_title(r"The coefficients $K_{n}$ and $M_{n}$")
for i, c in enumerate(coeffs):
print(c)
ax.scatter(n, c, linestyle=self.line_palette[i])
ax.set_yscale('log')
# ax.set_xscale('log')
ax.legend([r"$K_n$", r"$M_n$"])
# ax.semilogy(n, abs(B_n_real), linestyle=self.line_palette[1])
# ax.semilogy(n, abs(K_n_real), linestyle=self.line_palette[2])
# ax.semilogy(n, abs(M_n_real), linestyle=self.line_palette[3])
plt.show()
def test_arf(self) -> None:
......@@ -275,73 +288,199 @@ class TestYosiokaARF(BaseTest):
def test_special_solutions(self) -> None:
wave_types = ["traveling", "standing"]
for w in wave_types:
# standing wave solutions are off somehow
if w == "traveling":
print(w)
for key in self.cls_dict:
self.cls_dict[key].wave_type = w
self.sphere_change_parameters()
self.assign_parameters()
# self.assertAlmostEqual(
# self.cls_gorkov.acoustic_radiation_force(),
# self.cls_sphere.acoustic_radiation_force(), 3e-1)
# self.assertAlmostEqual(
# self.cls_general.acoustic_radiation_force(),
# self.cls_sphere.acoustic_radiation_force(), 3e-1)
self.bubble_change_parameters()
self.assign_parameters()
# self.assertAlmostEqual(
# self.cls_gorkov.acoustic_radiation_force(),
# self.cls_bubble.acoustic_radiation_force(), 3e-1)
# self.assertAlmostEqual(
# self.cls_general.acoustic_radiation_force(),
# self.cls_bubble.acoustic_radiation_force(), 3e-1)
def test_arf_plots(self) -> None:
for key in self.cls_dict:
self.cls_dict[key].wave_type = w
self.sphere_change_parameters()
self.assign_parameters()
self.assertAlmostEqual(
self.cls_gorkov.acoustic_radiation_force(),
self.cls_sphere.acoustic_radiation_force(), 3e-1)
self.assertAlmostEqual(
self.cls_general.acoustic_radiation_force(),
self.cls_sphere.acoustic_radiation_force(), 3e-1)
self.bubble_change_parameters()
self.assign_parameters()
self.assertAlmostEqual(
self.cls_gorkov.acoustic_radiation_force(),
self.cls_bubble.acoustic_radiation_force(), 3e-1)
self.assertAlmostEqual(
self.cls_general.acoustic_radiation_force(),
self.cls_bubble.acoustic_radiation_force(), 3e-1)
def test_small_spheres(self) -> None:
plt.style.use("ggplot")
fig, ax = plt.subplots()
ax.set_xlabel("Position")
ax.set_ylabel("Amplitude of the ARF")
ax.set_title("ARF on different locations")
ax.set_xlabel(r"$(k_{f} a)^{2}$")
ax.set_ylabel(r"Deviation to general solution in \%")
ax.set_title(r"Small Sphere Approximation")
self.sphere_change_parameters()
self.wave_type = "traveling"
self.assign_parameters()
radius = np.linspace(4e-3, 1e-9, int(1e2))
sol_general = []
sol_sphere = []
ka_squared = []
for r in radius:
self.R_0 = r
self.assign_parameters()
ka_squared.append((self.R_0 * self.compute_k_f()) ** 2)
sol_general.append(self.cls_general.acoustic_radiation_force())
sol_sphere.append((self.cls_sphere.acoustic_radiation_force()))
sol_sphere_normalized = [
100 * abs(sol_sphere[i] - sol_general[i])
/ sol_general[i] for i in range(len(sol_sphere))
]
ax.semilogx(
ka_squared,
sol_sphere_normalized,
label="small spheres",
color="b")
ax.set_ylim(-10, 100)
ax.set_xlim(1e-3, 1.05)
plt.legend()
plt.show()
def test_bubbles(self) -> None:
plt.style.use("ggplot")
fig, ax = plt.subplots()
ax.set_xlabel(r"$(k_{f} a)^{2}$")
ax.set_ylabel(r"Deviation to general solution in \%")
ax.set_title(r"Bubble Approximation")
self.bubble_change_parameters()
self.wave_type = "standing"
self.assign_parameters()
numbers = np.linspace(1, 100, 100)
radi = self.compute_xlambda() / self.f
radius = [n * radi for n in numbers]
sol_general = []
sol_bubble = []
ka_squared = []
for r in radius:
self.R_0 = r
self.assign_parameters()
ka_squared.append((self.R_0 * self.compute_k_f()) ** 2)
sol_general.append(self.cls_general.acoustic_radiation_force())
sol_bubble.append((self.cls_bubble.acoustic_radiation_force()))
sol_bubble_normalized = [
abs((sol_bubble[i] - sol_general[i])
/ sol_general[i]) for i in range(len(sol_bubble))
]
ax.plot(ka_squared,
sol_bubble_normalized,
label="bubbles approximation - normalized",
color="b")
# ax.set_xlim(0)
ax.relim()
plt.legend()
plt.show()
print(self.compute_xlambda())
def test_frequency_sweep(self) -> None:
"""
Recreates Fig. 2 from the Yosioka paper
"""
plt.style.use("ggplot")
fig, ax = plt.subplots()
ax.set_xlabel(r"$\frac{k^{*} a}{\sqrt{3 \lambda}}$")
ax.set_ylabel(r"Amplitude of the ARF")
ax.set_title("Bubble Approximation")
self.bubble_change_parameters()
self.wave_type = "traveling"
self.assign_parameters()
frequency = np.linspace(1e3, 2e5, int(1e2))
sol_general = []
sol_bubble = []
x_data = []
for f in frequency:
self.f = f
self.assign_parameters()
x_data.append(
(self.R_0 * self.compute_k_s())
/ (np.sqrt(3 * self.compute_xlambda())))
sol_general.append(self.cls_general.acoustic_radiation_force())
sol_bubble.append(self.cls_bubble.acoustic_radiation_force())
ax.semilogy(
x_data,
sol_bubble,
label="bubble approximation",
color="b",
linestyle="-.")
ax.semilogy(
x_data,
sol_general,
label="general solution",
color="r",
linestyle=":")
plt.legend()
plt.show()
def test_arf_traveling(self) -> None:
self.sphere_change_parameters()
self.assign_parameters()
plt.style.use("ggplot")
fig, ax = plt.subplots()
ax.set_xlabel(r"Position")
ax.set_ylabel(r"Amplitude of the ARF (N)")
ax.set_title(r"ARF on different locations")
pos = np.linspace(-1e-2, 1e-2, int(1e3))
legend = []
for key in self.cls_dict:
line_palette = ["-", "-.", ":", "--", "dotted", "dashdot"]
for i, key in enumerate(self.cls_dict):
try:
self.cls_dict[key].wave_type = "traveling"
if key == "bubble" or key == "king" or key == "special":
raise NotImplementedError
arf = []
for p in pos:
self.cls_dict[key].position = p
arf.append(
abs(self.cls_dict[key].acoustic_radiation_force()))
ax.semilogy(pos, arf)
self.cls_dict[key].acoustic_radiation_force())
ax.plot(pos, arf, linestyle=line_palette[i])
legend.append(key)
except NotImplementedError:
pass
ax.legend(legend)
plt.show()
for key in self.cls_dict:
self.cls_dict[key].wave_type = "standing"
try:
a = PlotAcoustofluidics(ARF=self.cls_dict[key])
print(f' This is : {key} for the standing wavefield')
a.plot_ARF()
except NotImplementedError:
pass
for key in self.cls_dict:
self.cls_dict[key].wave_type = "traveling"
def test_arf_standing(self) -> None:
self.sphere_change_parameters()
self.assign_parameters()
plt.style.use("ggplot")
fig, ax = plt.subplots()
ax.set_xlabel(r"Position (m)")
ax.set_ylabel(r"Amplitude of the ARF (N)")
ax.set_title(r"ARF on different locations")
pos = np.linspace(-1e-2, 1e-2, int(1e3))
legend = []
line_palette = ["-", "-.", ":", "--", "dotted", "dashdot"]
for i, key in enumerate(self.cls_dict):
try:
a = PlotAcoustofluidics(ARF=self.cls_dict[key])
print(f' This is : {key} for the traveling wavefield')
a.plot_ARF()
if key == "bubble" or key == "king":
raise NotImplementedError
self.cls_dict[key].wave_type = "standing"
arf = []
for p in pos:
self.cls_dict[key].position = p
arf.append(
(self.cls_dict[key].acoustic_radiation_force()))
ax.plot(pos, arf, linestyle=line_palette[i])
legend.append(key)
except NotImplementedError:
pass
ax.legend(legend)
plt.show()
def test_str(self) -> None:
print(self.cls_general.__str__())
def test_repr(self) -> None:
print(self.cls_general.__repr__())
if __name__ == '__main__':
unittest.main()
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment