Weak form: Implementing source term - Strain history

JoelTomasek · December 2, 2021, 9:56am

Dear all,

I would like to implement the following equation in a Navier Stokes/Convection-diffusion solver in FEniCS:

The goal of this equation is to get an insight in the “strain-history” of particles inside the velocity field.To that end, the scalar stress is defined as the Frobenius norm of the rate of strain tensor, which is added as a source term in the convection diffusion equation.

I tried to implement this, but I am not certain if my approach and specifically the variational formulation is correct. Could you tell me if at first glance I made any mistakes, or whether this is the correct approach? I added the code below:


V = VectorFunctionSpace(mesh, 'P', 2)
Q = FunctionSpace(mesh, 'P', 1)
S = FunctionSpace(mesh,'P',1)

# Define trial and test functions
u = TrialFunction(V)
v = TestFunction(V)
p = TrialFunction(Q)
q = TestFunction(Q)
s = TrialFunction(S)
z = TestFunction(S)

# Define functions for solutions at previous and current time steps
u_n = Function(V)
u_  = Function(V)
p_n = Function(Q)
p_  = Function(Q)
s_n = Function(S)
s_ = Function(S)

# Define symmetric gradient
def epsilon(u):
    return sym(nabla_grad(u))

def strainNorm(u):
    return Constant(sqrt(2.)*4e-3)*sqrt(inner(epsilon(u),epsilon(u)))

# Define stress tensor
def sigma(u, p):
    return 2*mu*epsilon(u) - p*Identity(len(u))
	
# Define expressions used in variational forms
U  = 0.5*(u_n + u)
n  = FacetNormal(mesh)
f  = Constant((0, 0))
fs  = Constant(0.)
k  = Constant(dt)
mu = Constant(mu)
rho = Constant(rho)

# Define variational problem for step 1
F1 = rho*dot((u - u_n) / k, v)*dx \
   + rho*dot(dot(u_n, nabla_grad(u_n)), v)*dx \
   + inner(sigma(U, p_n), epsilon(v))*dx \
   + dot(p_n*n, v)*ds - dot(mu*nabla_grad(U)*n, v)*ds \
   - dot(f, v)*dx
a1 = lhs(F1)
L1 = rhs(F1)


# Define variational problem for step 2
a2 = dot(nabla_grad(p), nabla_grad(q))*dx
L2 = dot(nabla_grad(p_n), nabla_grad(q))*dx - (1/k)*div(u_)*q*dx

# Define variational problem for step 3
a3 = dot(u, v)*dx
L3 = dot(u_, v)*dx - k*dot(nabla_grad(p_ - p_n), v)*dx

# Define variational problem for step 4

F4 = rho*((s - s_n) / k)*z*dx \
   + rho*dot(u_n, grad(s))*z*dx \
   - fs*z*dx
a4 = lhs(F4)
L4 = rhs(F4)

# Assemble matrices
A1 = assemble(a1)
A2 = assemble(a2)
A3 = assemble(a3)
# Apply boundary conditions to matrices
[bc.apply(A1) for bc in bcu]
[bc.apply(A2) for bc in bcp]

# Time-stepping
t = 0
for n in range(num_steps):

    # Update current time
    t += dt

    # Step 1: Tentative velocity step
    b1 = assemble(L1)
    [bc.apply(b1) for bc in bcu]
    solve(A1, u_.vector(), b1, 'bicgstab', 'hypre_amg')

    # Step 2: Pressure correction step
    b2 = assemble(L2)
    [bc.apply(b2) for bc in bcp]
    solve(A2, p_.vector(), b2, 'bicgstab', 'hypre_amg')

    # Step 3: Velocity correction step
    b3 = assemble(L3)
    solve(A3, u_.vector(), b3, 'cg', 'sor')
	
    # Step 4: Scalar solve step
    A4 = assemble(a4)
    [bc.apply(A4) for bc in bcs]
    b4 = assemble(L4)
    [bc.apply(b4) for bc in bcs]
    solve(A4,s_.vector(),b4)

    # Update previous solution
    u_n.assign(u_)
    p_n.assign(p_)
    s_n.assign(s_)

kamensky · December 2, 2021, 3:41pm

If F4 is meant to correspond to the advection equation in the post, you would need to replace fs with \eta instead of Constant(0). Treating the velocity explicitly, as in the advection term of F4, this would be

fs = strainNorm(u_n)

JoelTomasek · December 3, 2021, 9:53am

Thank you for your reply and checking the implementation. Indeed I meant to write fs = strainNorm(u_n), but I forgot to make the edit.

JoelTomasek · December 3, 2021, 10:28am

Upon running this implementation in a benchmark case, I notice that I get much coarser results for the scalar solution than the results from this paper Potential fluid mechanic pathways of platelet activation | SpringerLink
With the implementation above, is it true that I only calculate strain for one point in a cell?

If I want to get a higher resolution for the strain computation and the scalar solution, is it sufficient to change S to FunctionSpace(mesh,‘P’, 2)?

kamensky · December 3, 2021, 3:47pm

I get much coarser results for the scalar solution than the results from this paper

It looks like the linked paper is using Lagrangian particles instead of an Eulerian discretization of advection. This is more complicated to implement, but introduces less diffusion. (See, e.g., LEoPart for a FEniCS-based implementation of Lagrangian particles.)

If I want to get a higher resolution for the strain computation and the scalar solution, is it sufficient to change S to FunctionSpace(mesh,‘P’, 2)?

I don’t think this is likely to improve results without first working on other aspects of the discretization. The first-order time discretization may be introducing significant error and, for pure advection, it’s best to include some form of stabilization in the spatial discretization.

JoelTomasek · December 4, 2021, 12:42pm

Thank you for your suggestions. I have decided to take a look at the implementation of LEoPart. For my full application, I am using Oasis to solve incompressible NS in optimized fashion (GitHub - mikaem/Oasis). Is it possible to implement LEoPart in this formulation, or does it require that I solve NS with a discontinuous Galerkin approach?

In the meantime I will attempt to introduce stabilization in the spatial discretization and improve time discretization for the Eularian approach.
Below I attached my new version of the code, with the addition of SUPG stabilization terms.
I never implemented SUPG before, I am currently performing some tests to see how this implementation performs. Would you say that the code below provides a correct implementation?

# Define function spaces
V = VectorFunctionSpace(mesh, 'P', 2)
Q = FunctionSpace(mesh, 'P', 1)
S = FunctionSpace(mesh,'P',1)

# Define inlet flow profile
plug_velocity=0.02e-6/(math.pi*0.001375**2)
inflow_profile = ('2*plug_velocity*velocity_scaling*((0.001375-x[1])*(0.001375+x[1])) / pow(0.001375, 2)', '0')
vely=Expression(inflow_profile, plug_velocity=plug_velocity,velocity_scaling=velocity_scaling, degree=2)

# Define boundary conditions
bcu_inflow = DirichletBC(V, vely, inflow)
bcu_walls = DirichletBC(V, Constant((0, 0)), walls)
bcp_outflow = DirichletBC(Q, Constant(0), outflow)
bcs_inflow = DirichletBC(S, Constant(0.),inflow)
bcu = [bcu_inflow, bcu_walls]
bcp = [bcp_outflow]
bcs = [bcs_inflow]

# Define trial and test functions
u = TrialFunction(V)
v = TestFunction(V)
p = TrialFunction(Q)
q = TestFunction(Q)
s = TrialFunction(S)
z = TestFunction(S)

# Define functions for solutions at previous and current time steps
u_n = Function(V)
u_  = Function(V)
p_n = Function(Q)
p_  = Function(Q)
s_n = Function(S)
s_ = Function(S)

# Define symmetric gradient
def epsilon(u):
    return sym(nabla_grad(u))
def strainNorm(u):
    return Constant(sqrt(2.)*4e-3)*sqrt(inner(epsilon(u),epsilon(u)))
# Define stress tensor
def sigma(u, p):
    return 2*mu*epsilon(u) - p*Identity(len(u))

# Define expressions used in variational forms
U  = 0.5*(u_n + u)
s_mid = 0.5*(s_n + s)
n  = FacetNormal(mesh)
f  = Constant((0, 0))
fs  = strainNorm(u_n)
k  = Constant(dt)
mu = Constant(mu)
rho = Constant(rho)
h = CellDiameter(mesh)

# Define variational problem for step 1
F1 = rho*dot((u - u_n) / k, v)*dx \
   + rho*dot(dot(u_n, nabla_grad(u_n)), v)*dx \
   + inner(sigma(U, p_n), epsilon(v))*dx \
   + dot(p_n*n, v)*ds - dot(mu*nabla_grad(U)*n, v)*ds \
   - dot(f, v)*dx
a1 = lhs(F1)
L1 = rhs(F1)

# Define variational problem for step 2
a2 = dot(nabla_grad(p), nabla_grad(q))*dx
L2 = dot(nabla_grad(p_n), nabla_grad(q))*dx - (1/k)*div(u_)*q*dx

# Define variational problem for step 3
a3 = dot(u, v)*dx
L3 = dot(u_, v)*dx - k*dot(nabla_grad(p_ - p_n), v)*dx

# Define variational problem for step 4
F4 = ((s - s_n) / k)*z*dx \
   + dot(u_, grad(s_mid))*z*dx \
   - fs*z*dx
   
#Add SUPG stabilization terms
vnorm = sqrt(dot(u_, u_))
# Residual
r = (s - s_n)/k + (dot(u_, grad(s_mid)) - fs)
F4 += (h/(2.0*vnorm))*dot(u_, grad(z))*r*dx

a4 = lhs(F4)
L4 = rhs(F4)

# Assemble matrices
A1 = assemble(a1)
A2 = assemble(a2)
A3 = assemble(a3)

# Apply boundary conditions to matrices
[bc.apply(A1) for bc in bcu]
[bc.apply(A2) for bc in bcp]

s=s_n

# Time-stepping
t = 0
for n in range(num_steps):

    # Update current time
    t += dt
    if n%int(0.05*1/dt)==0:
        print(t)
    flowsec=0.04e-6
    baseflow = 0.06e-6
    vely.plug_velocity = (baseflow+flowsec*math.sin(4*math.pi*t-0.5*math.pi))/(math.pi*0.001375**2)

    # Step 1: Tentative velocity step
    b1 = assemble(L1)
    [bc.apply(b1) for bc in bcu]
    solve(A1, u_.vector(), b1, 'bicgstab', 'hypre_amg')

    # Step 2: Pressure correction step
    b2 = assemble(L2)
    [bc.apply(b2) for bc in bcp]
    solve(A2, p_.vector(), b2, 'bicgstab', 'hypre_amg')

    # Step 3: Velocity correction step
    b3 = assemble(L3)
    solve(A3, u_.vector(), b3, 'cg', 'sor')
	
    # Step 4: Scalar solve step
    A4 = assemble(a4)
    [bc.apply(A4) for bc in bcs]
    b4 = assemble(L4)
    [bc.apply(b4) for bc in bcs]
    solve(A4,s_.vector(),b4)

    # Update previous solution
    u_n.assign(u_)
    p_n.assign(p_)
    s_n.assign(s_)

Topic		Replies	Views
Weak formulation for square of gradient variational formulation	2	1348	August 15, 2019
Challenging variational formulation variational formulation	1	364	April 24, 2020
Nonlinear terms with many derivatives variational formulation	4	854	March 16, 2021
Help with a weak form	4	494	December 2, 2021
Variation formulation for a second order PDE with first order derivative term variational formulation	1	397	March 9, 2021

Weak form: Implementing source term - Strain history

Related Topics