Hello, everyone, I am working in fenics to solve a complex fsi problem and encounter some problems. Here is a simple example :
I need to solve a fluid equation (N-S equation) and a structure equation (line elastic equation) .So i build two mesh . First i get the fluid stress s1 on the mesh_f ,next writing s1 on the elastic equation 's right hand .There is my code
# FSI simple example
from fenics import *
from mshr import *
import numpy as np
import matplotlib.pyplot as plt
# physical constants
rho_f = 1.1
mu = 1e-3
# mesh
domain1 = Rectangle(Point(0.0, 0.0), Point(6.0, 1.0))
#domain2 = Circle(Point(0.5, 0.5), 0.15)
domain3 = Rectangle(Point(2.95, 0.0), Point(3.05, 0.5))
mesh_f = generate_mesh(domain1 - domain3, 30)
mesh = generate_mesh(domain1, 40)
plot(mesh_f)
plt.show()
n_f = FacetNormal(mesh_f)
u_elem = VectorElement("P", mesh_f.ufl_cell(), 2)
p_elem = FiniteElement('P', mesh_f.ufl_cell(), 1)
up_elem = MixedElement([u_elem, p_elem])
V_f = FunctionSpace(mesh_f, up_elem)
u, p = TrialFunctions(V_f)
v, q = TestFunctions(V_f)
up_old = Function(V_f)
u_old, p_old = split(up_old)
up_fem = Function(V_f)
def walls(x, on_boundary):
return on_boundary and near(x[1], 0) or near(x[1] , 1.0)
def in_flow(x, on_boundary):
return on_boundary and near(x[0], 0)
def out_flow(x, on_boundary):
return on_boundary and near(x[0] , 6.0)
def rectangle(x, on_boundary):
return on_boundary and near(x[0], 2.95) or near(x[0], 3.05) or near(x[1], 0.5)
#u_in = Expression(('100*x[1]*(1.0 - x[1])*pow(cos(pi*10*t), 2)/pow(0.5, 2)', '0.0'), degree = 2, t= 0.0)
u_in = Expression(('6.0*x[1]*(1.0 - x[1])/pow(0.5, 2)', '0.0'), degree = 2)
bcu_in = DirichletBC(V_f.sub(0), u_in, in_flow)
bcu_walls = DirichletBC(V_f.sub(0), Constant((0.0, 0.0)), walls)
bcu_cylinder = DirichletBC(V_f.sub(0), Constant((0.0, 0.0)), rectangle)
bcp_out = DirichletBC(V_f.sub(1), Constant(0.0), out_flow)
bcu = [bcu_in, bcu_walls, bcu_cylinder, bcp_out]
def epsilon(u):
return sym(nabla_grad(u))
def sigma_f(u, p):
return 2*mu*epsilon(u) - p*Identity(2)
rho_f = Constant(rho_f)
mu = Constant(mu)
f = Constant((0.0, 0.0))
F_ns = rho_f*dot(dot(u_old, nabla_grad(u)), v)*dx + \
rho_f*dot(dot(u, nabla_grad(u_old)), v)*dx + \
inner(epsilon(u), epsilon(v))*dx -\
p*div(v)*dx+\
div(u)*q*dx - \
dot(f, v)*dx -\
rho_f*dot(dot(u_old, nabla_grad(u_old)), v)*dx
a_f, L_f = lhs(F_ns), rhs(F_ns)
# stress function space
VT = TensorFunctionSpace(mesh_f, 'P', 2)
solve(a_f == L_f, up_fem, bcu)
u_fem, p_fem = up_fem.split()
stress = sigma_f(u_fem, p_fem)
s1 = project(stress, VT)
print("go, on!", s1(2.90, 0.5))
plot(u_fem)
plt.show()
# fluid stress s1
# write s1 in the elastic equation's right hand
mesh_s = generate_mesh(domain3, 10)
plot(mesh_s)
plt.show()
n_s = FacetNormal(mesh_s)
n_f = FacetNormal(mesh_f)
# physical constant
# physical constants
mu_s = 1e3
lambda1 = 1.25
# solve a stable elastic beam problem
# get the displacement data
# functionspace
V_s = VectorFunctionSpace(mesh_s, 'P', 1)
eta_s = TrialFunction(V_s)
phi_s = TestFunction(V_s)
eta_s_fem = Function(V_s)
# boundary condition Dirichlet on the bottom
class Bottom_s(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[1], 0.0)
bottom_s = Bottom_s()
bc_s = DirichletBC(V_s, Constant((0.0, 0.0)), bottom_s)
# weak form
def sigma_s(eta):
"""
elastic stress tensor
"""
return 2*mu_s*epsilon(eta) + lambda1*tr(epsilon(eta))*Identity(2)
a_s = inner(sigma_s(eta_s), epsilon(phi_s))*dx
#a = 2*mu*inner(epsilon(eta_s), epsilon(phi_s))*dx + \
# lambda1*innre(tr(epsilon(eta_s))*Identity(2), epsilon(phi_s))*dx
s1.set_allow_extrapolation(True)
L_s = dot(dot(s1, n_f), phi_s)*ds
# solve problem
solve(a_s == L_s, eta_s_fem, bc_s)
plot(eta_s_fem, title='Displacement', mode='displacement')
plt.show()
But i get an error that
→ 150 L_s = dot(dot(s1, n_f), phi_s)*ds
151 # solve problem
152 solve(a_s == L_s, eta_s_fem, bc_s)
/usr/local/lib/python3.6/dist-packages/ufl/measure.py in rmul(self, integrand)
435 domain = self.ufl_domain()
436 if domain is None:
→ 437 domains = extract_domains(integrand)
438 if len(domains) == 1:
439 domain, = domains
/usr/local/lib/python3.6/dist-packages/ufl/domain.py in extract_domains(expr)
353 for t in traverse_unique_terminals(expr):
354 domainlist.extend(t.ufl_domains())
→ 355 return sorted(join_domains(domainlist))
356
357
TypeError: ‘<’ not supported between instances of ‘Mesh’ and ‘Mesh’
Thank you very much for your help!