From 08d088fc15d3c160cf2bce64c4a86a8de1c4a043 Mon Sep 17 00:00:00 2001 From: John Date: Wed, 1 Apr 2020 14:55:13 -0500 Subject: [PATCH] added hydrodynamic interface code --- stoked/__init__.py | 2 +- stoked/hydrodynamics.py | 35 +++++++++++++++++++++++++++++++++++ stoked/solver.py | 4 +++- 3 files changed, 39 insertions(+), 2 deletions(-) diff --git a/stoked/__init__.py b/stoked/__init__.py index 2e5487b..d188ff0 100644 --- a/stoked/__init__.py +++ b/stoked/__init__.py @@ -14,7 +14,7 @@ from .electrodynamics import point_dipole_electrodynamics from .van_der_waals import van_der_waals_sphere, van_der_waals_sphere_interface from .gravity import gravity, gravity_sphere, gravity_ellipsoid -from .hydrodynamics import grand_mobility_matrix +from .hydrodynamics import grand_mobility_matrix, interface from .constraints import constrain_position, constrain_rotation from .forces import pairwise_central_force, pairwise_force from .common import lennard_jones diff --git a/stoked/hydrodynamics.py b/stoked/hydrodynamics.py index 658160f..1aad723 100644 --- a/stoked/hydrodynamics.py +++ b/stoked/hydrodynamics.py @@ -1,5 +1,14 @@ import numpy as np +class interface: + """A no-slip interface""" + def __init__(self, z=0): + """ + Arguments: + z z-position of the interface + """ + self.z = z + def levi_civita(): """return the levi-civita symbol""" @@ -8,6 +17,32 @@ def levi_civita(): eijk[0, 2, 1] = eijk[2, 1, 0] = eijk[1, 0, 2] = -1 return eijk +def particle_wall_self_mobility(position, interface, viscosity, radius): + """ + Construct the particle wall self-mobility matrix for a single particle + + Arguments: + position[3] position of particle + interface interface object + viscosity dynamic viscosity ยต of surrounding fluid + radius particle radius + """ + M = np.zeros([2, 2, 3, 3], dtype=float) + h = position[2] - interface.z + + gamma_T = 6*np.pi*viscosity*radius + gamma_R = 6*np.pi*viscosity*radius**3 + + a = 1/(16*gamma_T)*(9/h - 2/h**3 + 1/h**5) + b = 1/(8*gamma_T)*(9/h - 4/h**3 + 1/h**5) + M[0,0] = np.diag([a,a,b]) + + a = 15/(64*gamma_R)*(1/h**3) + b = 3/(32*gamma_R)*(1/h**3) + M[1,1] = np.diag([a,a,b]) + + return M + def grand_mobility_matrix(position, drag_T, drag_R, viscosity): """ Construct the grand mobility matrix for a given cluster diff --git a/stoked/solver.py b/stoked/solver.py index 72e9667..6a38b50 100644 --- a/stoked/solver.py +++ b/stoked/solver.py @@ -63,7 +63,7 @@ class stokesian_dynamics: """ def __init__(self, *, temperature, dt, position, drag, orientation=None, force=None, torque=None, interactions=None, constraint=None, - hydrodynamic_coupling=True): + interface=None, hydrodynamic_coupling=True): """ Arguments: temperature system temperature @@ -75,6 +75,7 @@ def __init__(self, *, temperature, dt, position, drag, orientation=None, torque(t, r, q) external torque function given time t, position r[N,D], orientation q[N] and returns torque T[N,D] (can be a list of functions) interactions particle interactions (can be a list) constraint particle constraints (can be a list) + interface no-slip boundary interface (default: no interface) hydrodynamic_coupling if True, include hydrodynamic coupling interactions (default: True) """ self.temperature = temperature @@ -84,6 +85,7 @@ def __init__(self, *, temperature, dt, position, drag, orientation=None, self.drag = drag self.Nparticles = len(self.position) self.ndim = self.position.shape[1] + self.interface = interface self.hydrodynamic_coupling = hydrodynamic_coupling if orientation is not None: