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