Tau constraints on coupled boundaries (Ball and Shell)

[Thomas Michl]

Hello everyone,

for my bachelors id like to simulate an active droplet using dedalus. (The paper with the governing equations https://www.annualreviews.org/content/journals/10.1146/annurev-conmatphys-031115-011517). My idea to achieve this was to use spherical coordinates with a ball and shell basis, to simulate the inner and outer fields of the droplet. The inner and outer fields are coupled by matching shear-stress, pressure difference and matching flow fields on the boundary.

I have coded something, but i have trouble adding the tau terms, such that a square and solveable system is created. I have attached my (non-working) code.

I have tried some things with the tau terms, but non have created a solvable system. I am also not sure, what gauge conditions i need (i think i need atleast one pressure gauge, for the inner or the outer pressure).

Any ideas on how this can be achieved would be great.

dtype = np.float64
eta = 10
Nphi, Ntheta, Nr = 64,64,64

R_inner = 2
R_outer = 3000
scale = 1
scales = (scale,scale,scale)
dealias = 1


coords = d3.SphericalCoordinates('phi', 'theta', 'r')
dist = d3.Distributor(coords, dtype=dtype, mesh = None)
shell = d3.ShellBasis(coords, shape=(Nphi, Ntheta, Nr), radii=(R_inner, R_outer), dtype=dtype)
sphere = shell.outer_surface
sphere_outer = shell.inner_surface

ball = d3.BallBasis(coords, shape=(Nphi, Ntheta, Nr), radius=R_inner, dealias=dealias, dtype=dtype)
sphere_inner = ball.surface



# Fields_shell
u = dist.VectorField(coords, name='u', bases = shell)
p = dist.Field(name = 'p', bases = shell)
tau_1 = dist.VectorField(coords, name='tau_1', bases=sphere )
tau_p = dist.Field(name = 'tau_p')
tau_2 = dist.VectorField(coords, name='tau_2', bases=sphere)

# Fields_ball
u_2 = dist.VectorField(coords, name='u_2', bases = ball)
p_2 = dist.Field(name = 'p_2', bases = ball)
tau_1_2 = dist.VectorField(coords, name='tau_1_2', bases=sphere_inner )
tau_p_2 = dist.Field(name = 'tau_p_2')
tau_2_2 = dist.VectorField(coords, name='tau_2_2', bases=sphere_inner)

#unit vectors, coords Shell
phi, theta, r = dist.local_grids(shell)
er = dist.VectorField(coords, bases=shell)
ephi = dist.VectorField(coords, bases=shell)
etheta = dist.VectorField(coords, bases=shell)

#substitutions shell
rvec=dist.VectorField(coords, bases=shell.radial_basis)
rvec['g'][2] = r

lift_basis = shell.derivative_basis(1)
lift = lambda A, n: d3.Lift(A, lift_basis, n)
grad_u = d3.grad(u) + rvec*lift(tau_2, -1)

#unit vectors, coords Ball
phi_2, theta_2, r_2 = dist.local_grids(ball)
er_2 = dist.VectorField(coords, bases=ball)
ephi_2 = dist.VectorField(coords, bases=ball)
etheta_2 = dist.VectorField(coords, bases=ball)

#subs ball
rvec_2=dist.VectorField(coords, bases=ball.radial_basis)
rvec_2['g'][2] = r_2

lift_basis_2 = ball.derivative_basis(1)
lift_2 = lambda A, n: d3.Lift(A, lift_basis_2, n)
grad_u_2 = d3.grad(u_2) + rvec_2*lift_2(tau_2_2, -2)



#boundary flow
x_hat = dist.VectorField(coords, name='x_hat', bases=sphere)

x_hat['g'][2] = -np.cos(theta)  # r-component
x_hat['g'][1] = -np.sin(theta)  # theta-component
x_hat['g'][0] = 0                # phi-component


Shear_inner = d3.angular(d3.radial((eta*(d3.grad(u)+d3.grad(u).T))(r=R_inner), index =1))
Shear_outer= d3.angular(d3.radial((eta*(d3.grad(u_2)+d3.grad(u_2).T))(r=R_inner), index =1))

Gamma = 1

#Problem
problem = d3.LBVP([p,u,p_2,u_2,tau_1, tau_2, tau_p, tau_p_2,tau_1_2], namespace = locals())
#equations shell
problem.add_equation("-grad(p) + eta * div(grad_u)+ lift(tau_1, -1) = 0")
problem.add_equation("trace(grad_u) + tau_p= 0")
problem.add_equation("radial(u(r=R_inner)) = 0")
problem.add_equation("u(r=R_outer) = x_hat")
problem.add_equation("integ(p)=0")

#interface
problem.add_equation("Shear_outer - Shear_inner = 0")
problem.add_equation("-p(r=R_inner) + p_2(r=R_inner) = -2*Gamma/R_inner")
problem.add_equation("angular(u_2(r=R_inner)) - angular(u(r=R_inner)) =0")

#equation ball
problem.add_equation("-grad(p_2) + eta * lap(u_2) +  lift_2(tau_1_2, -1)= 0")
problem.add_equation("div(u_2) + tau_p_2 = 0")
problem.add_equation("radial(u_2(r=R_inner)) = 0")
problem.add_equation("integ(p_2)=0")




solver = problem.build_solver()
solver.solve()

Sincerely
Thomas