# Weak form: Implementing source term - Strain history

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)

def epsilon(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 \
+ 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

# 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 \
- 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_)



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)


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.

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)?

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.

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)

def epsilon(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 \
+ 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

# 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 \
- fs*z*dx

vnorm = sqrt(dot(u_, u_))
# Residual
r = (s - s_n)/k + (dot(u_, grad(s_mid)) - fs)

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_)