Hello,
I’m working on a free surface problem using the ALE method, and so far I’ve been able to update the surface curvature and normal vector by redefining my variational problem at each time step. I apologize as couldn’t figure out how to develop a MWE for this application, so I kindly ask to refer here for the complete code. The code is still being developed, so it is still computationally inneficient (and a little fuzzy). I’m using the fenics tutorial for the navier-stokes equation using the IPCS approach, found here. I have the following stress boundary condition to represent the free passive surface:
Currently I’m redefining the variational problem at each time step using the python function below to apply the boundary condition:
# Interface stress
def IST(sigma, kappa, n):
return sigma*kappa*n
The variational formulation becomes
F1 = rho*dot((u - u_n) / k, v)*dx \
+ rho*dot(dot(u_n, nabla_grad(u_n)), v)*dx \
+ inner(TT(U, p_n, mu, I), DD(v))*dx \
- dot(mu*nabla_grad(U)*n, v)*ds \
- dot(f, v)*dx \
+ dot(p_n*n, v)*ds
b_int = inner(IST(sigma, kappa, n), v)*ds(surface_marker)
F1 = F1 - b_int
By using this method I generate physically coherent results, such as the velocity field and pressure field of a perturbed interface as described below. The results present a pressure difference due to the surface curvature, expressing the capillarity effects that are expected for this transport problem. The capillary effects will thus induce a fluid displacement in both on the surface and the bulk of the domain.
Velocity field:
Pressure field:
Now i’m trying to implement a UserExpression to be updated at each time step, so that I don’t need to redefine the variational problem for every step of the temporal loop. I’m following this post, in special this solution to implement a UserExpression class. I am also studying this post on how to update time-dependent boundary conditions, as it is also related to my question. I came up to the following class for my problem:
class InterfaceTension(UserExpression):
def __init__(self, sigma, kappa, n, **kwargs):
super().__init__(kwargs)
self.sigma = sigma
self.kappa = kappa
self.n = n
def eval(self, values, x):
values = self.sigma * self.kappa * self.n
def update(self, t, u):
self.kappa = kappa
self.n = n
def value_shape(self):
return (2,)
Then I instantiate the class by
IST = InterfaceTension(sigma, kappa, n)
The variational form then becomes
F1 = rho*dot((u - u_n) / k, v)*dx \
+ rho*dot(dot(u_n, nabla_grad(u_n)), v)*dx \
+ inner(TT(U, p_n, mu, I), DD(v))*dx \
- dot(mu*nabla_grad(U)*n, v)*ds \
- dot(f, v)*dx \
+ dot(p_n*n, v)*ds
b_int = inner(IST, v)*ds(surface_marker)
F1 = F1 - b_int
After solving the problem for the first time-step, I need to update the new curvature and normal vector associated with the surface (inside the time loop):
# inside the temporal loop
kappa, n = curvature(mesh, ds, surface_marker) # see github to access this function
IST.kappa = kappa
IST.n = n
This approach returns a non-physical result, different from the results obtained by using the first approach (redefining the variational problem at each time step):
Velocity field from the second approach:
Pressure field from the second approach:
Does anyone know what can be happening to create this difference between the results obtained using the two approaches? Honestly I have no idea why this is happening. I am very grateful for any insights.
Thanks in advance.