Nonlinear variational problem does not converge (Newton method)

Question

Hello,

The nonlinear variational problem that I am dealing with is the 1D model of a liquid jet stretched by gravity, developed by Eggers & Dupont (1994) and García & Castellanos (1994), and which reads,

$\frac{\partial h^2}{\partial t} + \frac{\partial (h^2 u)}{\partial z} = 0$,
$\frac{\partial u}{\partial t} + u \frac{\partial u}{\partial z} = 1 - \frac{\partial }{\partial z} \left( \nabla \cdot \mathbf{n}\right) + \frac{\Gamma}{h^2} \frac{\partial}{\partial z} \left(h^2 \frac{\partial u}{\partial z} \right)$,
where $h(z,t)$ is the radius of the jet, $u(z,t)$ its axial velocity and $\nabla \cdot \mathbf{n} = h^{-1} \, \left[ 1 + \left( \frac{\partial h}{\partial z} \right)^2 \right]^{-1/2} - \frac{\partial^2 h}{\partial z^2} \, \left[ 1 + \left( \frac{\partial h}{\partial z} \right)^2 \right]^{-3/2}$, is the full curvature of the liquid-air interface. The boundary conditions that need to be imposed are: $h(0,t) = \sqrt{Bo}$, and $u(0,t) = We^{1/2} Bo^{-1/4}$, where $We$ and $Bo$ are the Weber and the Bond number respectively.

In particular I am trying to solve the stationary state of the latter system of equations, which is a third order nonlinear ODE. From the first equation (continuity) we can obtain the relation between the radius and the axial velocity, $h^2 u = q = Bo^{3/4} We^{1/2}$, where $q$ is the dimensionless flow rate, being its value a function of the Weber and the Bond number. If we introduce this expression into the stationary axial momentum equation we obtain a third order nonlinear ODE,

$-\frac{2 q^2}{h^5} h' + (\nabla \cdot \mathbf{n})' - \frac{\Gamma}{h^2} \left( \frac{-2 q}{h} h' \right)' -1 = 0$, $h(0) = \sqrt{Bo}$ (the Bond number is fixed, e.g. $Bo =1$).

I have tried to implement this problem in FEniCS through a nonlinear variational problem, but I always obtain NaNs just in the first iteration. I attach the code,

from fenics import *

# Dimensionless parameters
Bo    = Constant(1.0)
We    = Constant(8E-3)
Gamma = Constant(16.67)
q     = (Bo**(3/4))*sqrt(We)

# Create mesh and define function space
mesh = IntervalMesh(60,0, 20)
V    = FunctionSpace(mesh, 'P', 2)

# Define boundary condition
inflow   = 'near(x[0],0)'
bc       = DirichletBC(V,sqrt(Bo),inflow)


# Define variational problem
h = Function(V)
v = TestFunction(V)

F = 2*(q**2)*(h**(-4))*grad(v)[0]*dx -10*(q**2)*h**(-5)*grad(h)[0]*v*dx-((h**(-1))*(1+grad(h)[0]**2)**(-1/2)-grad(grad(h)[0])[0]*(1+grad(h)[0]**2)**(-3/2))*grad(v)[0]*dx - Gamma*2*q*grad(h)[0]*(h**(-3))*grad(v)[0]*dx + Gamma*4*q*(grad(h)[0]**2)*(h**(-4))*v*dx + ((q**2)*h**(-4))*v*ds +((h**(-1))*(1+grad(h)[0]**2)**(-1/2)-grad(grad(h)[0])[0]*(1+grad(h)[0]**2)**(-3/2))*v*ds + Gamma*2*q*grad(h)[0]*(h**(-3))*v*ds - v*dx

J = derivative(F,h)

problem = NonlinearVariationalProblem(F, h, bc, J=J)
solver  = NonlinearVariationalSolver(problem)
prm = solver.parameters
prm["newton_solver"]["absolute_tolerance"] = 1E-6
prm["newton_solver"]["relative_tolerance"] = 1E-6
prm["newton_solver"]["maximum_iterations"] = 50
prm["newton_solver"]["relaxation_parameter"] = 0.1
solver.solve()

#solve(F == 0, h, bc, J=J)

# Plot solution and mesh
plot(h)

# Hold plot
interactive()

Thanks in advance,
Alejandro

Answer 1

The initial guess is very important for problems this horrifically nonlinear.

Note that you have terms like $ h^{-1} $. Ensure that you set a non-zero initial guess as follows:

h = Function(V)
h.vector()[:] = 1.0

or by interpolating an Expression, for example,

h = interpolate(Expression("x[0] + 1.0", element=V.ufl_element()), V)

I'm still unable to achieve convergence, but I'm optimistic that you know more about the system to be able to determine a more intelligently designed initial guess.

Comments

Thank you so much Nate! Now I have achived convergence. Starting from $h(z) = \sqrt{Bo}/(1+z)$ (balance between surface tension force with slender curvature approximation and gravity), the Newton algorithm converges.

I have another question, but related to the variational formulation. Since I have a Dirichlet boundary condition at $z = 0$, the contribution to the variational problem of this boundary disappears due to the test function has to be zero, and the contribution of the other boundary, $z = L$, remains. I have implemented this in FEniCS as follows (you can also see it in the code attached in the other post),

f(h)*f(v)*ds

Where f(h) and f(v) are functions of the jet radius and the test function respectively. ¿Is this correct?

Thanks in advance,
Alejandro