Errors solving a 1D problem (fiber spinning)

Question

Hello,
I am completely new to FEniCS (I have just test the tutorial examples), and now I am trying to solve a fluid mechanics 1D problem, in particular the 1D model developed in the field of fiber spinning in the mid 60s, which reads,

$\frac{\partial h}{\partial t} + \frac{\partial(u \, h)}{\partial x} = 0, \, \, 0\leq x \leq 1,$
$\frac{\partial}{\partial x} \left( h \frac{\partial u}{\partial x} \right) = 0, \, \, 0\leq x \leq 1,$
$h(0,t) = u(0,t) =1, \, \, u(1,t) = D_r,$

where $D_r$ is the draw ratio and it is a constant, and the initial conditions are the stationary state of the 1D equations,

$u(x,0) = D_r^x, \, \, h(x,0) = D_r^{-x}$

My main issue with this kind of 1D problems is that I always get mismatch errors. In particular here I get this type of error when solving a_1 == L_1. I guess this is due to "dtdiv(htest)u_nhdx", but I do not know any other form to define the derivative of $h(x,t)$ respect to $x$.

Besides, I do not know if it is appropriate to define two function spaces, one for the velocity $u(x,t)$ and another one for the fiber radius $h(x,t)$.

I attach the code,

from __future__ import print_function
from fenics import *

T = 2.0            # final time
num_steps = 1000     # number of time steps
dt = T / num_steps # time step size

# Dimensionless variables
q     = 1

# Create mesh and define function space
mesh = IntervalMesh(30,0,5)
V    = FunctionSpace(mesh, 'P', 2)
Q    = FunctionSpace(mesh, 'P', 2)

# Define boundary condition
inflow   = 'near(x[0],0)'
outflow  = 'near(x[0],5)'

bch_inflow   = DirichletBC(V, Constant(1), inflow)
bcu_inflow   = DirichletBC(Q, Constant(1), inflow)
bcu_outflow  = DirichletBC(Q, Constant(23), outflow)

bcu = [bcu_inflow, bcu_outflow]
bch = [bch_inflow]

# Define the initial value
h_0 = Expression('pow(Dr,-x[0])', degree=2,Dr = 23)
h_n = interpolate(h_0,V)
u_0 = Expression('pow(Dr,x[0])', degree=2,Dr = 23)
u_n = interpolate(u_0,Q)

# Define variational problem
h     = TrialFunction(V)  
htest = TestFunction(V)
u     = TrialFunction(Q)
utest = TestFunction(Q)
a1 = h*htest*dx - dt*div(htest)*u_n*h*dx
L1 = h_n*htest*dx

a2 =-h_n*inner(grad(u),grad(utest))*dx

# Time-stepping
h = Function(V)
u = Function(Q)
t = 0

for n in range(num_steps):
    # Update current time
    t += dt
    # Compute solution
    solve(a1 == L1, h, bch)
    solve(a2 == 0, u, bcu)
    #Plot solution

    plot(u)
    print('t = %.2f' % (t))
    # Update previous solution
    u_n.assign(u)
    h_n.asign(h)
# Hold plot

interactive()

Thanks in advance,
Alejandro

Answer 1

I don't really know what to expect from a solution of your problem, so I am unable to validate my answer.

Regarding your issues with div(htest) I refer you to here.

Furthermore, ufl formulations of the type a2 == 0 will assume that a2 is a semiliear residual arising from the variational formulation of a nonlinear problem.

You also had a typo in .assign.

Consider the following:

from __future__ import print_function
from fenics import *

T = 2.0            # final time
num_steps = 1000     # number of time steps
dt = T / num_steps # time step size

# Dimensionless variables
q     = 1

# Create mesh and define function space
mesh = IntervalMesh(30,0,5)
V    = FunctionSpace(mesh, 'P', 2)
Q    = FunctionSpace(mesh, 'P', 2)

# Define boundary condition
inflow   = 'near(x[0],0)'
outflow  = 'near(x[0],5)'

bch_inflow   = DirichletBC(V, Constant(1), inflow)
bcu_inflow   = DirichletBC(Q, Constant(1), inflow)
bcu_outflow  = DirichletBC(Q, Constant(23), outflow)

bcu = [bcu_inflow, bcu_outflow]
bch = [bch_inflow]

# Define the initial value
h_0 = Expression('pow(Dr,-x[0])', degree=2,Dr = 23)
h_n = interpolate(h_0,V)
u_0 = Expression('pow(Dr,x[0])', degree=2,Dr = 23)
u_n = interpolate(u_0,Q)

# Define variational problem
h     = TrialFunction(V)  
htest = TestFunction(V)
u     = TrialFunction(Q)
utest = TestFunction(Q)

a1 = h*htest*dx - dt*grad(htest)[0]*u_n*h*dx
L1 = h_n*htest*dx

a2 =-h_n*inner(grad(u),grad(utest))*dx
L2 = Constant(0.0)*utest*dx

# Time-stepping
h = Function(V)
u = Function(Q)
t = 0

for n in range(num_steps):
    # Update current time
    t += dt
    # Compute solution
    solve(a1 == L1, h, bch)
    solve(a2 == L2, u, bcu)
    #Plot solution

    plot(u)
    print('t = %.2f' % (t))
    # Update previous solution
    u_n.assign(u)
    h_n.assign(h)
# Hold plot

interactive()

Comments

Thank you so much, now the code is running without any error messages! However the program does not converge. According to the mathematical model presented above, the variational problem reads,

$\int_0^1 \frac{\partial \, h}{\partial t} h_{\text{test}} \, \text{d}x - \int_0^1 u \, h \frac{\partial h_{\text{test}}}{\partial x} \, \text{d}x + u\,h \, h_{\text{test}}|_{1} = 0$

$-\int_0^1 h \, \frac{\partial u}{\partial x} \frac{\partial u_{\text{test}}}{\partial x} \, \text{d} x = 0.$

If the initial condition is $h(x,0) = D_r^{-x}$, $u(x,0) D_r^x$, the result that I expect from this equations is that the solution has to be absolutely stable, due to this initial condition is the analytical equilibrium state of the model. If this condition is slightly perturbed, with a value of $D_r$ above $22$ (do not know the exact value), the solution has to oscillate eventually reaching a limit cycle. However, as I said before the program is not numerically stable.

I attach you the last version of the code,

from __future__ import print_function
from fenics import *

T = 2.0            # final time
num_steps = 10000     # number of time steps
dt = T / num_steps # time step size

# Create mesh and define function space
mesh = IntervalMesh(57,0,1)
V    = FunctionSpace(mesh, 'P', 2)
Q    = FunctionSpace(mesh, 'P', 2)

# Define boundary condition
inflow   = 'near(x[0],0)'
outflow  = 'near(x[0],1)'

bch_inflow   = DirichletBC(V, Constant(1), inflow)
bcu_inflow   = DirichletBC(Q, Constant(1), inflow)
bcu_outflow  = DirichletBC(Q, Constant(23), outflow)

bcu = [bcu_inflow, bcu_outflow]
bch = [bch_inflow]

# Define the initial value
h_0 = Expression('pow(Dr,-x[0])', degree=2,Dr = 23)
h_n = interpolate(h_0,V)
u_0 = Expression('pow(Dr,x[0])', degree=2,Dr = 23)
u_n = interpolate(u_0,Q)

# Define variational problem
h     = TrialFunction(V)  
htest = TestFunction(V)
u     = TrialFunction(Q)
utest = TestFunction(Q)

a1 = h*htest*dx- dt*grad(htest)[0]*u_n*h*dx + dt*u_n*h*htest*ds(1)
L1 = h_n*htest*dx

a2 =-grad(u)[0]*grad(utest)[0]*h_n*dx
L2 = Constant(0.0)*utest*dx

# Time-stepping
h = Function(V)
u = Function(Q)
t = 0

for n in range(num_steps):
    # Update current time
    t += dt
    # Compute solution
    solve(a1 == L1, h, bch)
    solve(a2 == L2, u, bcu)
    #Plot solution

    plot(h)
    print('t = %.2f' % (t))
    # Update previous solution
    u_n.assign(u)
    h_n.assign(h)
# Hold plot

interactive()

Thanks in advance,
Alejandro