Discussion Forum

Hi all,

I have been struggling with the following problem for a long time, so I'd appreciate any help.

The problem consists of (a) determine the free surface of a liquid in a cylindrical container as you rotate it (think about a small glass of water that you rotate) and (b) studying its stability.

I've managed to solve (a) but I don't know which way to go for (b).

For (a), I use a stationary PDE modes in a circle (2D) to solve for the interface of the free surface, z = H0(x,y).

The formulation for (b) is to insert to the Stokes flow equations a perturbation around the steady state found in (a). This results in an eigenvalue problem for the 3 components of the perturbed velocity, (u(x,y,z),v(x,y,z),w(x,y,z)), perturbed pressure p(x,y,z) and perturbed free interface hh(x,y) in a cylinder C whose top face is given by z = H0(x,y).
The problem looks like this:
px = uxx + uyy + uzz
py = vxx + vyy + vzz
pz = wxx + wyy + wzz
0 = ux + vy + wz

The perturbed interface satisfies: lambda*hh(x,y) = w(x,y,H0) - u(x,y,H0)*H0x - v(x,y,H0)*H0y, hh = 0 at the boundary (circle).

And the boundary conditions for the 3D problem for u, v, w, p are
1) u = v = w = 0 at the container boundaries (all cylinder boundaries expect the top, which is open)
2) -p + N(H0,u,v,w) = kappa(H0,hh), at z = H0
3) T1(H0,u,v,w) = 0, at z = H0
4) T2(H0,u,v,w) = 0, at z = H0

1) represents the no slip condition at the solid boundaries, and 2), 3) and 4) are the normal and shear stress balances accross the free interface, respectively.

To solve problem b), I first create in Matlab the cylindrical geometry C using the solution of a) and then import it into Comsol.
I have attempted to solve it in two ways:

B1) Treating hh as an unknown: I use a multiphysics approach coupling a 3D geometry with a 2D geometry, and use extrusion couple variables. The eigenvalue only appears in the equation for hh. Model eigstab2D3D.mph (attached) implements this (to start with, in here I was playing with the simplest case in which H0 = 0, flat interface).
I manage to obtain a solution but I think it's wrong since in the source term in the equation for hh does not depend on hh, so the solver might be ignoring it. (page 258 from Modeling Guide says "If the general or weak solution forms are used, the source terms are not ignored if they depend on the solution components.") And I'm not sure the solver interprets an extrusion variable as a solution component...

B2) Since hh only appears in the boundary condition 2), I write hh(x, y) = (w(x,y,H0) - u(x,y,H0)*H0x - v(x,y,H0)*H0) / lambda, insert this in 2) and multiply the whole expression by lambda (for the case lambda = 0). This results in an eigenvalue problem for u, v, w, and p only in the 3D geometry, for which lambda only appears in one of the boundary conditions, but not in the equations. Thus, the mass matrix is zero, and comsol doesn't seem to deal with this kind of problems. This is implemented in model geomcyl.mph (this one has a H0 nontrivial).

I would be very very grateful for any hint or advise that you could offer.

Best wishes,
Maria