In the FeniCSx tutorial for linear elasticity, they have calculated the Von Mises stress as follows: (link)
s = sigma(uh) -1./3*ufl.tr(sigma(uh))*ufl.Identity(uh.geometric_dimension())
von_Mises = ufl.sqrt(3./2*ufl.inner(s, s))
V_von_mises = fem.FunctionSpace(domain, ("DG", 0))
stress_expr = fem.Expression(von_Mises, V_von_mises.element.interpolation_points())
stresses = fem.Function(V_von_mises)
stresses.interpolate(stress_expr)
I’m trying to calculate the fluid shear stress by following the same steps above.
That means:
In place of
s = sigma(uh) -1./3*ufl.tr(sigma(uh))*ufl.Identity(uh.geometric_dimension())
von_Mises = ufl.sqrt(3./2*ufl.inner(s, s))
I have (The calculation you need for the shear stresses):
def epsilon(u):
return sym(nabla_grad(u))
# Define stress tensor
def sigma(u, p):
return 2*mu*epsilon(u) - p*Identity(u.geometric_dimension())
n = FacetNormal(mesh)
#Stress tensor. Comes from the theory
T = -sigma(u_n, p_n)*n
#Normal stress (A scalar)
Tn = inner(T, n)
#Shear stress (A vector)
Tt = T - Tn*n
And then, for the function space given in the Vonn Meses stress:
V_von_mises = fem.FunctionSpace(domain, ("DG", 0))
I have the:
p_cg1 = FiniteElement('CG', mesh.ufl_cell(), 1)
scalar = FunctionSpace(mesh, p_cg1) #For the normal stress
v_cg2 = VectorElement('CG', mesh.ufl_cell(), 2)
vector = FunctionSpace(mesh, v_cg2) #For the shear stress
normal_stress = Function(scalar)
shear_stress = Function(vector)
Finally, in place of the calculation and the projection in Von Meses stress:
stress_expr = fem.Expression(von_Mises, V_von_mises.element.interpolation_points())
stresses = fem.Function(V_von_mises)
stresses.interpolate(stress_expr)
I will have (For the normal stress):
normal_stress_expr = Expression(Tn, scalar.element.interpolation_points())
normal_stress.interpolate(normal_stress_expr)
And I should have a similar thing for the shear stresses too.
However, even for this normal stress calculation, I’m getting the following error:
Looks like, there should be a different way to define the facet normals (Given by n)
I would really appreciate your suggestions.
I’m giving my entire code below if you need to reproduce everything. It is a bit lengthy, but it is just the exact copy (Except for the stress calculation) of the tutorial code in Poiseuille flow
import numpy as np
import gmsh
import sys
from dolfinx import plot
# Initialize gmsh:
from mpi4py import MPI
#from dolfinx import mesh
import pyvista
from dolfinx.io.gmshio import model_to_mesh
import matplotlib.pyplot as plt
from dolfinx import geometry
# from dolfinx import fem,mesh
import ufl
from petsc4py import PETSc
from prettytable import PrettyTable
import random
from dolfinx.mesh import create_unit_square
from ufl import (FacetNormal, FiniteElement, FacetArea,Identity,TestFunction, TrialFunction, VectorElement,
div, dot, ds, dx, inner, lhs, nabla_grad, rhs, sym)
from dolfinx.fem import Expression,Constant,Function, FunctionSpace, assemble_scalar, dirichletbc, form, locate_dofs_geometrical,locate_dofs_topological
from dolfinx.fem.petsc import assemble_matrix, assemble_vector, apply_lifting, create_vector, set_bc
from dolfinx.io import XDMFFile
from dolfinx.plot import create_vtk_mesh
mesh = create_unit_square(MPI.COMM_WORLD, 10, 10)
t = 0
T = 10
num_steps = 500
dt = T/num_steps
v_cg2 = VectorElement("CG", mesh.ufl_cell(), 2)
s_cg1 = FiniteElement("CG", mesh.ufl_cell(), 1)
V = FunctionSpace(mesh, v_cg2)
Q = FunctionSpace(mesh, s_cg1)
u = TrialFunction(V)
v = TestFunction(V)
p = TrialFunction(Q)
q = TestFunction(Q)
def walls(x):
return np.logical_or(np.isclose(x[1],0), np.isclose(x[1],1))
wall_dofs = locate_dofs_geometrical(V, walls)
u_noslip = np.array((0,) * mesh.geometry.dim, dtype=PETSc.ScalarType)
bc_noslip = dirichletbc(u_noslip, wall_dofs, V)
def inflow(x):
return np.isclose(x[0], 0)
inflow_dofs = locate_dofs_geometrical(Q, inflow)
bc_inflow = dirichletbc(PETSc.ScalarType(8), inflow_dofs, Q)
def outflow(x):
return np.isclose(x[0], 1)
outflow_dofs = locate_dofs_geometrical(Q, outflow)
bc_outflow = dirichletbc(PETSc.ScalarType(0), outflow_dofs, Q)
bcu = [bc_noslip]
bcp = [bc_inflow, bc_outflow]
u_n = Function(V)
u_n.name = "u_n"
U = 0.5 * (u_n + u)
n = FacetNormal(mesh)
f = Constant(mesh, PETSc.ScalarType((0,0)))
k = Constant(mesh, PETSc.ScalarType(dt))
mu = Constant(mesh, PETSc.ScalarType(1))
rho = Constant(mesh, PETSc.ScalarType(1))
# Define strain-rate tensor
def epsilon(u):
return sym(nabla_grad(u))
# Define stress tensor
def sigma(u, p):
return 2*mu*epsilon(u) - p*Identity(u.geometric_dimension())
# Define the variational problem for the first step
p_n = Function(Q)
p_n.name = "p_n"
F1 = rho*dot((u - u_n) / k, v)*dx
F1 += rho*dot(dot(u_n, nabla_grad(u_n)), v)*dx
F1 += inner(sigma(U, p_n), epsilon(v))*dx
F1 += dot(p_n*n, v)*ds - dot(mu*nabla_grad(U)*n, v)*ds
F1 -= dot(f, v)*dx
a1 = form(lhs(F1))
L1 = form(rhs(F1))
A1 = assemble_matrix(a1, bcs=bcu)
A1.assemble()
b1 = create_vector(L1)
# Define variational problem for step 2
u_ = Function(V)
a2 = form(dot(nabla_grad(p), nabla_grad(q))*dx)
L2 = form(dot(nabla_grad(p_n), nabla_grad(q))*dx - (1/k)*div(u_)*q*dx)
A2 = assemble_matrix(a2, bcs=bcp)
A2.assemble()
b2 = create_vector(L2)
# Define variational problem for step 3
p_ = Function(Q)
a3 = form(dot(u, v)*dx)
L3 = form(dot(u_, v)*dx - k*dot(nabla_grad(p_ - p_n), v)*dx)
A3 = assemble_matrix(a3)
A3.assemble()
b3 = create_vector(L3)
# Solver for step 1
solver1 = PETSc.KSP().create(mesh.comm)
solver1.setOperators(A1)
solver1.setType(PETSc.KSP.Type.BCGS)
pc1 = solver1.getPC()
pc1.setType(PETSc.PC.Type.HYPRE)
pc1.setHYPREType("boomeramg")
# Solver for step 2
solver2 = PETSc.KSP().create(mesh.comm)
solver2.setOperators(A2)
solver2.setType(PETSc.KSP.Type.BCGS)
pc2 = solver2.getPC()
pc2.setType(PETSc.PC.Type.HYPRE)
pc2.setHYPREType("boomeramg")
# Solver for step 3
solver3 = PETSc.KSP().create(mesh.comm)
solver3.setOperators(A3)
solver3.setType(PETSc.KSP.Type.CG)
pc3 = solver3.getPC()
pc3.setType(PETSc.PC.Type.SOR)
xdmf = XDMFFile(mesh.comm, "poiseuille.xdmf", "w")
xdmf.write_mesh(mesh)
xdmf.write_function(u_n, t)
xdmf.write_function(p_n, t)
for i in range(num_steps):
# Update current time step
t += dt
# Step 1: Tentative veolcity step
with b1.localForm() as loc_1:
loc_1.set(0)
assemble_vector(b1, L1)
apply_lifting(b1, [a1], [bcu])
b1.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
set_bc(b1, bcu)
solver1.solve(b1, u_.vector)
u_.x.scatter_forward()
# Step 2: Pressure corrrection step
with b2.localForm() as loc_2:
loc_2.set(0)
assemble_vector(b2, L2)
apply_lifting(b2, [a2], [bcp])
b2.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
set_bc(b2, bcp)
solver2.solve(b2, p_.vector)
p_.x.scatter_forward()
# Step 3: Velocity correction step
with b3.localForm() as loc_3:
loc_3.set(0)
assemble_vector(b3, L3)
b3.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
solver3.solve(b3, u_.vector)
u_.x.scatter_forward()
# Update variable with solution form this time step
u_n.x.array[:] = u_.x.array[:]
p_n.x.array[:] = p_.x.array[:]
# Write solutions to file
xdmf.write_function(u_n, t)
xdmf.write_function(p_n, t)
# Close xmdf file
xdmf.close()
#Stress calculation:
n = FacetNormal(mesh)
T = -sigma(u_n, p_n)*n
Tn = inner(T, n)
Tt = T - Tn*n
p_cg1 = FiniteElement('CG', mesh.ufl_cell(), 1)
scalar = FunctionSpace(mesh, p_cg1)
v_cg2 = VectorElement('CG', mesh.ufl_cell(), 2)
vector = FunctionSpace(mesh, v_cg2)
normal_stress = Function(scalar)
shear_stress = Function(vector)
normal_stress_expr = Expression(Tn, scalar.element.interpolation_points())
normal_stress.interpolate(normal_stress_expr)
# shear_stress_expr = Expression(Tt, vector.element.interpolation_points())
# shear_stress.interpolate(shear_stress_expr)
