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

analytical potentials and velocities

parent 4cb20500
......@@ -2,10 +2,9 @@ from __future__ import annotations
from typing import Union, Optional
from gorkov import (
log, pi, sin, cos,
log, pi, sin, cos, conj,
BaseSphereFrequencyComposite,
Sphere, Frequency, integrate)
import numpy as np
from gorkov.solutions.base_arf import BaseARF
from .baseclass import BaseYosioka
from gorkov.solutions.Yosioka1955.scatteringfield import ScatteringField
......@@ -107,20 +106,20 @@ class ARF(BaseYosioka, BaseARF, BaseSphereFrequencyComposite):
self.n += 1
# (45)
s1 = 0.5 * (K_n_prev * np.conj(self.K_n_0)).real
s1 = 0.5 * (K_n_prev * conj(self.K_n_0)).real
s1 *= 2 * self.n
s1 *= (-2 * pi * self.R_0 ** 2 * self.rho_f)
out[0] += s1
# (46)
s2 = 0.5 * (M_n_prev * np.conj(self.M_n_0)).real
s2 = 0.5 * (M_n_prev * conj(self.M_n_0)).real
s2 *= 2 * (self.n - 1) * self.n * (self.n + 1)
s2 *= (2 * pi * self.xlambda ** 2 * self.rho_f)
out[1] += s2
# (47)
s3 = 0.5 * (K_n_prev * np.conj(self.M_n_0)).real
s4 = 0.5 * (M_n_prev * np.conj(self.K_n_0)).real
s3 = 0.5 * (K_n_prev * conj(self.M_n_0)).real
s4 = 0.5 * (M_n_prev * conj(self.K_n_0)).real
s3 *= 2 * self.n * (self.n + 1)
s3 *= (-2 * pi * self.R_0 * self.xlambda * self.rho_f)
s4 *= 2 * (self.n - 1) * self.n
......@@ -129,7 +128,7 @@ class ARF(BaseYosioka, BaseARF, BaseSphereFrequencyComposite):
out[3] += s4
# (48)
s5 = 0.5 * (M_n_prev * np.conj(self.M_n_0)).real
s5 = 0.5 * (M_n_prev * conj(self.M_n_0)).real
s5 *= 2 * self.n
s5 *= (-2 * pi * (self.k_f ** 2) *
(self.R_0 ** 2) * (self.xlambda ** 2) * self.rho_f)
......
......@@ -3,8 +3,8 @@ from gorkov import (ActiveVariable,
Frequency, Sphere,
LegendreFunctions as leg,
SpecialFunctions as sp,
exp, log)
exp, log, cos, sin)
from scipy.special import lpmv
from gorkov.solutions.Yosioka1955 import BaseYosioka
from gorkov.solutions.base_scattering import BaseScattering
......@@ -35,6 +35,30 @@ class ScatteringField(BaseYosioka, BaseScattering):
"velocity potential of scattered waves")
self._Phi_star = ActiveVariable(self._compute_Phi_star,
"velocity potential inside particle")
self._d_Phi_1_dr = ActiveVariable(
self._compute_d_Phi_1_dr,
"Derivative of Phi_1 w.r.t. r")
self._d_Phi_1_dtheta = ActiveVariable(
self._compute_d_Phi_1_dtheta,
"Derivative of Phi_1 w.r.t. theta")
self._d_Phi_i_dr = ActiveVariable(
self._compute_d_Phi_i_dr,
"Derivative of Phi_i w.r.t. r")
self._d_Phi_i_dtheta = ActiveVariable(
self._compute_d_Phi_i_dtheta,
"Derivative of Phi_i w.r.t. theta")
self._d_Phi_s_dr = ActiveVariable(
self._compute_d_Phi_s_dr,
"Derivative of Phi_s w.r.t. r")
self._d_Phi_s_dtheta = ActiveVariable(
self._compute_d_Phi_s_dtheta,
"Derivative of Phi_s w.r.t. theta")
self._d_Phi_star_dr = ActiveVariable(
self._compute_d_Phi_star_dr,
"Derivative of Phi_star w.r.t. r")
self._d_Phi_star_dtheta = ActiveVariable(
self._compute_d_Phi_star_dtheta,
"Derivative of Phi_star w.r.t. theta")
self._radial_acoustic_fluid_velocity = ActiveVariable(
self.radial_acoustic_fluid_velocity, "radial acoustic velocity")
self._radial_particle_velocity = ActiveVariable(
......@@ -48,40 +72,95 @@ class ScatteringField(BaseYosioka, BaseScattering):
# Dependencies:
self._Phi_1.is_computed_by(
self._Phi_i, self._Phi_s)
self._d_Phi_1_dr.is_computed_by(
self._d_Phi_i_dr, self._d_Phi_s_dr)
self._d_Phi_1_dtheta.is_computed_by(
self._d_Phi_i_dtheta, self._d_Phi_s_dtheta)
self._Phi_i.is_computed_by(
self._delta_n, self.fluid._k_f,
self._r, self._theta, self._t)
self._d_Phi_i_dtheta.is_computed_by(
self._delta_n, self.fluid._k_f,
self._r, self._theta, self._t)
self._d_Phi_i_dr.is_computed_by(
self._delta_n, self.fluid._k_f,
self._r, self._theta, self._t)
self._Phi_s.is_computed_by(
self._A_n, self._C_n, self.fluid._k_f,
self._r, self._theta, self._t)
self._d_Phi_s_dr.is_computed_by(
self._A_n, self._C_n, self.fluid._k_f,
self._r, self._theta, self._t)
self._d_Phi_s_dtheta.is_computed_by(
self._A_n, self._C_n, self.fluid._k_f,
self._r, self._theta, self._t)
self._Phi_star.is_computed_by(
self._A_n, self._D_n, self.fluid._k_f,
self._r, self._theta, self._t)
self._d_Phi_star_dr.is_computed_by(
self._A_n, self._D_n, self.fluid._k_f,
self._r, self._theta, self._t)
self._d_Phi_star_dtheta.is_computed_by(
self._A_n, self._D_n, self.fluid._k_f,
self._r, self._theta, self._t)
self._radial_particle_velocity.is_computed_by(
self._Phi_star, self._theta, self._t, self._r)
self._d_Phi_star_dr)
self._tangential_particle_velocity.is_computed_by(
self._Phi_star, self._theta, self._t, self._r)
self._d_Phi_star_dtheta, self._r)
self._radial_acoustic_fluid_velocity.is_computed_by(
self._Phi_1, self._theta, self._t, self._r)
self._d_Phi_1_dr)
self._tangential_acoustic_fluid_velocity.is_computed_by(
self._Phi_1, self._theta, self._t, self._r)
self._d_Phi_1_dtheta, self._r)
@property
def Phi_1(self) -> complex:
return self._Phi_1.value
@property
def d_Phi_1_dr(self) -> complex:
return self._d_Phi_s_dr.value
@property
def d_Phi_1_dtheta(self) -> complex:
return self._d_Phi_1_dtheta.value
@property
def Phi_i(self) -> complex:
return self._Phi_i.value
@property
def d_Phi_i_dr(self) -> complex:
return self._d_Phi_i_dr.value
@property
def d_Phi_i_dtheta(self) -> complex:
return self._d_Phi_i_dtheta.value
@property
def Phi_s(self) -> complex:
return self._Phi_s.value
@property
def d_Phi_s_dr(self) -> complex:
return self._d_Phi_s_dr.value
@property
def d_Phi_s_dtheta(self) -> complex:
return self._d_Phi_s_dtheta.value
@property
def Phi_star(self) -> complex:
return self._Phi_star.value
@property
def d_Phi_star_dr(self) -> complex:
return self._d_Phi_star_dr.value
@property
def d_Phi_star_dtheta(self) -> complex:
return self._d_Phi_star_dtheta.value
@property
def v_ac_r(self) -> float:
return self._radial_acoustic_fluid_velocity.value
......@@ -115,21 +194,21 @@ class ScatteringField(BaseYosioka, BaseScattering):
Compute velocity potential outside the sphere
according to (17)
"""
if self.r < self.R_0:
print("Warning: Phi_1 is only defined outside the particle!")
# (17)
return self.Phi_i + self.Phi_s
def _compute_d_Phi_1_dr(self) -> complex:
return self.d_Phi_i_dr + self.d_Phi_s_dr
def _compute_d_Phi_1_dtheta(self) -> complex:
return self.d_Phi_i_dtheta + self.d_Phi_s_dtheta
def _compute_Phi_i(self) -> complex:
"""
Compute the incident velocity potential
according to (16) for travelling waves
and (27) for standing waves
"""
if self.r < self.R_0:
print("Warning: Phi_i is only defined outside the particle!")
# (16)
if self.wave_type == "traveling":
coeff = self.legendre_coeffs(lambda n:
......@@ -149,15 +228,69 @@ class ScatteringField(BaseYosioka, BaseScattering):
coeff)
return out
def _compute_d_Phi_i_dr(self) -> complex:
"""
Compute the incident velocity potential
according to (16) for travelling waves
and (27) for standing waves
"""
# (16)
if self.wave_type == "traveling":
coeff = self.legendre_coeffs(lambda n:
(2 * n + 1) * (-1j) ** n
* self.k_f *
sp.d_jn(n, self.k_f * self.r))
# (27)
elif self.wave_type == "standing":
coeff = self.legendre_coeffs(lambda n:
(2 * n + 1) * (-1j) ** n
* self.delta_n
* self.k_f *
sp.d_jn(n, self.k_f * self.r))
else:
# No incident potential field
coeff = 0
raise NotImplementedError
out = exp(1j * (self.omega * self.t))
out *= leg.first_cos_poly(self.theta,
coeff)
return out
def _compute_d_Phi_i_dtheta(self) -> complex:
"""
Compute the incident velocity potential
according to (16) for travelling waves
and (27) for standing waves
"""
# (16)
if self.wave_type == "traveling":
coeff = self.legendre_coeffs(lambda n:
(2 * n + 1) * (-1j) ** n
* sp.jn(n, self.k_f * self.r))
# (27)
elif self.wave_type == "standing":
coeff = self.legendre_coeffs(lambda n:
(2 * n + 1) * (-1j) ** n
* self.delta_n
* sp.jn(n, self.k_f * self.r))
else:
# No incident potential field
coeff = 0
raise NotImplementedError
out = exp(1j * (self.omega * self.t))
degrees = len(coeff)
out *= sum([coeff[n] * (lpmv(2, n, cos(self.theta)))
for n in range(degrees)])
out *= -sin(self.theta)
out /= self.r
return out
def _compute_Phi_s(self) -> complex:
"""
Compute the scattered velocity potential
according to (18) for travelling waves
and (29) for standing waves
"""
if self.r < self.R_0:
print("Warning: Phi_s is only defined outside the particle!")
# (18)
if self.wave_type == "traveling":
coeff = self.legendre_coeffs(lambda n:
......@@ -171,19 +304,71 @@ class ScatteringField(BaseYosioka, BaseScattering):
else:
# No incident potential field
coeff = 0
raise NotImplementedError
out = exp(1j * (self.omega * self.t))
out *= leg.first_cos_poly(self.theta, coeff)
return out
def _compute_d_Phi_s_dr(self) -> complex:
"""
Compute the scattered velocity potential
according to (18) for travelling waves
and (29) for standing waves
"""
# (18)
if self.wave_type == "traveling":
coeff = self.legendre_coeffs(lambda n:
(2 * n + 1) * (-1j) ** n * self.A_n
* self.k_f *
sp.d_hn(self.n, self.k_f * self.r))
# (29)
elif self.wave_type == "standing":
coeff = self.legendre_coeffs(lambda n:
(2 * n + 1) * (-1j) ** n * self.C_n
* self.k_f *
sp.d_hn(self.n, self.k_f * self.r))
else:
# No incident potential field
coeff = 0
raise NotImplementedError
out = exp(1j * (self.omega * self.t))
out *= leg.first_cos_poly(self.theta, coeff)
return out
def _compute_d_Phi_s_dtheta(self) -> complex:
"""
Compute the scattered velocity potential
according to (18) for travelling waves
and (29) for standing waves
"""
# (18)
if self.wave_type == "traveling":
coeff = self.legendre_coeffs(lambda n:
(2 * n + 1) * (-1j) ** n * self.A_n
* sp.hn(self.n, self.k_f * self.r))
# (29)
elif self.wave_type == "standing":
coeff = self.legendre_coeffs(lambda n:
(2 * n + 1) * (-1j) ** n * self.C_n
* sp.hn(self.n, self.k_f * self.r))
else:
# No incident potential field
coeff = 0
raise NotImplementedError
out = exp(1j * (self.omega * self.t))
degrees = len(coeff)
out *= sum([coeff[n] * (lpmv(2, n, cos(self.theta)))
for n in range(degrees)])
out *= -sin(self.theta)
out /= self.r
return out
def _compute_Phi_star(self) -> complex:
"""
Compute the velocity potential inside the sphere
according to (19) for travelling waves
and (30) for standing waves
"""
if self.r > self.R_0:
print("Warning: Phi_star is only defined inside the particle!")
# (19)
if self.wave_type == "traveling":
coeff = self.legendre_coeffs(lambda n:
......@@ -197,10 +382,66 @@ class ScatteringField(BaseYosioka, BaseScattering):
else:
# No incident potential field
coeff = 0
raise NotImplementedError
out = exp(1j * (self.omega * self.t))
out *= leg.first_cos_poly(self.theta, coeff)
return out
def _compute_d_Phi_star_dr(self) -> complex:
"""
Compute the velocity potential inside the sphere
according to (19) for travelling waves
and (30) for standing waves
"""
# (19)
if self.wave_type == "traveling":
coeff = self.legendre_coeffs(lambda n:
(2 * n + 1) * (-1j) ** n * self.B_n
* self.k_s *
sp.d_jn(self.n, self.k_s * self.r))
# (30)
elif self.wave_type == "standing":
coeff = self.legendre_coeffs(lambda n:
(2 * n + 1) * (-1j) ** n * self.D_n
* self.k_s *
sp.d_jn(self.n, self.k_s * self.r))
else:
# No incident potential field
coeff = 0
raise NotImplementedError
out = exp(1j * (self.omega * self.t))
out *= leg.first_cos_poly(self.theta, coeff)
return out
def _compute_d_Phi_star_dtheta(self) -> complex:
"""
Compute the velocity potential inside the sphere
according to (19) for travelling waves
and (30) for standing waves
"""
# (19)
if self.wave_type == "traveling":
coeff = self.legendre_coeffs(lambda n:
(2 * n + 1) * (-1j) ** n * self.B_n
* sp.jn(self.n, self.k_s * self.r))
# (30)
elif self.wave_type == "standing":
coeff = self.legendre_coeffs(lambda n:
(2 * n + 1) * (-1j) ** n * self.D_n
* sp.jn(self.n, self.k_s * self.r))
else:
# No incident potential field
coeff = 0
raise NotImplementedError
out = exp(1j * (self.omega * self.t))
degrees = len(coeff)
out *= sum([coeff[n] * (lpmv(2, n, cos(self.theta)))
for n in range(degrees)])
out *= -sin(self.theta)
out /= self.r
return out
# -------------------------------------------------------------------------
# User-faced functions
# -------------------------------------------------------------------------
......@@ -218,14 +459,8 @@ class ScatteringField(BaseYosioka, BaseScattering):
"""
self.t = t
self.theta = theta
epsilon = r / 1000.0
self.r = r + epsilon
out = self.Phi_1
self.r -= epsilon
# Phi1 is dependent on nu, so it is updated
out -= self.Phi_1
out /= epsilon
return out.real
self.r = r
return self.d_Phi_1_dr.real
def tangential_acoustic_fluid_velocity(self, r: float, theta: float,
t: float) -> float:
......@@ -239,17 +474,9 @@ class ScatteringField(BaseYosioka, BaseScattering):
:rtype: float
"""
self.t = t
self.r = r
self.theta = theta
epsilon = 1.0 / 1000.0
self.theta = theta + epsilon
out = self.Phi_1
self.theta -= epsilon
# Phi1 is dependent on theta, so it is updated
out -= self.Phi_1
out /= epsilon
out /= self.r
return out.real
self.r = r
return self.d_Phi_1_dtheta.real / self.r
def radial_particle_velocity(self, r: float, theta: float,
t: float) -> float:
......@@ -264,14 +491,8 @@ class ScatteringField(BaseYosioka, BaseScattering):
"""
self.t = t
self.theta = theta
epsilon = r / 1000.0
self.r = r
out = self.Phi_star
self.r -= epsilon
# Phi_star is dependent on nu, so it is updated
out -= self.Phi_star
out /= epsilon
return out.real
return self.d_Phi_star_dr.real
def tangential_particle_velocity(self, r: float, theta: float,
t: float) -> float:
......@@ -288,16 +509,9 @@ class ScatteringField(BaseYosioka, BaseScattering):
:rtype: float
"""
self.t = t
self.theta = theta
self.r = r
epsilon = 1.0 / 1000.0
self.theta = theta + epsilon
out = self.Phi_star
self.theta -= epsilon
# Phi_star is dependent on theta, so it is updated
out -= self.Phi_star
out /= epsilon
out /= self.r
return out.real
return self.d_Phi_star_dtheta.real / self.r
def d_Phi1_dt(self, r: float, theta: float, t: float) -> float:
self.r = r
......
......@@ -11,31 +11,50 @@ class TestYosiokaARF(BaseTest):
def setUp(self) -> None:
self.c_f = 1e3
# Geometry
self.R_0 = 50e-6
# Frequency
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_s = 1e6
self.rho_s = 1e3
self.f = 1e6 / 2 / pi
self.p0 = 1e5
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
self.R0 = 1e-6
n_max = 3
self.cls = Yosioka1955.ARF(frequency=self.f,
radius=self.R0, rho_s=self.rho_s,
radius=self.R_0, rho_s=self.rho_s,
rho_f=self.rho_f, c_f=self.c_f,
c_s=self.c_s, p_0=self.p0,
position=self.pos,
wave_type="traveling",
particle_type="analytical", n_max=n_max)
self.cls_num = Yosioka1955.ARF(frequency=self.f,
radius=self.R0, rho_s=self.rho_s,
radius=self.R_0, rho_s=self.rho_s,
rho_f=self.rho_f, c_f=self.c_f,
c_s=self.c_s, p_0=self.p0,
position=self.pos,
wave_type="traveling",
particle_type="numerical", n_max=n_max)
self.cls_special = Yosioka1955.ARF(frequency=self.f,
radius=self.R0, rho_s=self.rho_s,
radius=self.R_0, rho_s=self.rho_s,
rho_f=self.rho_f, c_f=self.c_f,
c_s=self.c_s, p_0=self.p0,
position=self.pos,
......@@ -43,25 +62,25 @@ class TestYosiokaARF(BaseTest):
particle_type="special", n_max=n_max
)
self.cls_bubble = Yosioka1955.ARF(frequency=self.f,
radius=self.R0, rho_s=self.rho_s,
radius=self.R_0, rho_s=self.rho_s,
rho_f=self.rho_f, c_f=self.c_f,
c_s=self.c_s, p_0=self.p0,
position=self.pos,
wave_type="traveling",
particle_type="bubble", n_max=n_max)
self.cls_sphere = Yosioka1955.ARF(frequency=self.f,
radius=self.R0, rho_s=self.rho_s,
radius=self.R_0, rho_s=self.rho_s,
rho_f=self.rho_f, c_f=self.c_f,
c_s=self.c_s, p_0=self.p0,
position=self.pos,
wave_type="traveling",
particle_type="sphere", n_max=n_max)
self.clsKing = King1934.ARF(self.f,
self.R0, self.rho_s,
self.R_0, self.rho_s,
self.rho_f, self.c_f,
self.p0, self.pos)
self.clsGorkov = Gorkov1962.ARF(self.f,
self.R0, self.rho_s, self.c_s,
self.R_0, self.rho_s, self.c_s,
self.rho_f, self.c_f,
'traveling', self.p0, self.pos)
......@@ -82,7 +101,7 @@ class TestYosiokaARF(BaseTest):
self.rho_s = self.randomly_change_number(self.rho_s)
self.p0 = self.randomly_change_number(self.p0)