Axisymmetric problems: Nearly incompressible linear viscoelastic solid sphere with surface tension

General
1

Deal all,

I am currently trying to implement a nearly incompressible linear viscoelastic solid sphere with surface tension.
The viscoelastic part is running ok but introducing surface tension is creating problems, as the stress is higher at the top node and leads to wrong or diverging results.

The code is quite long, but the solver uses:

############################### VISCOELASTIC RESIDUAL ################################

    dx = ufl.Measure("dx",domain=domain,metadata={"quadrature_degree": deg_quad, "quadrature_scheme": "default"})
    Residual_viscoelastic = ufl.dot(eps_Mandel_axi(u_,r), sig) *2*np.pi*r* dx
    viscoelastic_tangent_form = ufl.dot(eps_Mandel_axi(v,r), ufl.dot(Ct, eps_Mandel_axi(u_,r))) *2*np.pi*r*dx

for viscoelasticity (I am using this tutorial for viscoelasticity).
Surface tension and later the solver write:

################################## SURFACE TENSION RESIDUAL #######################

    # Surface tension coefficient 
    sigma = fem.Constant(domain, surface_tension_value)

    # Deformation gradient
    F = ufl.Identity(2) + ufl.grad(u)

    # Reference normal to the boundary
    n = ufl.FacetNormal(domain)

    t_ref = ufl.as_vector([-n[1], n[0]])   # unit tangent, since n is unit

    # actual meridian stretch ratio (current length / reference length)
    Ft = ufl.dot(F, t_ref)
    J_surf = ufl.sqrt(ufl.dot(Ft, Ft))     

    # Current configuration: both radius and arc length change
    r_current = r + u[0]  # Current radius    

    Surface_tension_energy = sigma * 2 * np.pi  * r_current * J_surf * ds(1)     

    # Residual and tangent from surface tension
    Residual_surface_tension = ufl.derivative(Surface_tension_energy, u, u_)
    surface_tension_tangent_form = ufl.derivative(Residual_surface_tension, u, v)

The solver reads:

###################### PROBLEM AND RESOLUTION ############################
    
    tangent_problem = CustomLinearProblem(
        tangent_form,
        -Residual,
        u=du,
        bcs=bcs,
        petsc_options={
            "ksp_type": "preonly",
            "pc_type": "lu",
            "pc_factor_mat_solver_type": "mumps",
            "ksp_monitor": None,          # print linear residuals per iteration
            #"ksp_view": None,             # print KSP solver info
            #"ksp_converged_reason": None  # why linear solver converged or failed
        },
        petsc_options_prefix="tangent_"
    )
    
    newton = CustomNewtonSolver(tangent_problem)
    
    newton.callback = constitutive_update

with boundary conditions classical for axisymmetry.
A specific image summarizing the problem would be this one, at the beginning of the simulation:

image
image802×565 27 KB

Would you have any idea what to look for please ?

Note: the \theta\theta Mandel strain contribution is very slightly regularized to avoid divergence at r=0.