Hi everybody,
I’m trying to implement a simulation for the acoustic-elastic interaction in 2D. Thanks to this community and some articles on the subject, I was able to pose mathematically a simplified version of the problem I want to solve. Unfortunately, I’m still encountering some problems in order to run it on FEniCS. For a start, let us consider a solid domain \Omega_s inside a bounded fluid domain \Omega_a. The domains, boundaries and surface-normal vector are represented in this figure.
The problem is: given p on \partial B_{left}, find the pressure p and the displacement \mathbf{u} such that
\cases{ \Delta p + \kappa^2 p = 0 \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \text{in } \Omega_a\\
\nabla \mathbf{\sigma(u)} + \omega^2 \mathbf{u} = 0 \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \text{ in } \Omega_s\\
\partial _{\mathbf{n}_1} p = \rho_a \omega^2 \mathbf{n}_1 \cdot \mathbf{u} \text{ and} -p \mathbf{n}_1 = \mathbf{\sigma(u)} \cdot \mathbf{n}_1 \quad \space \space \text{on } \Gamma_s\\
p = 0 \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \space \space \text{on } \partial B \neq \partial B_{left}}
where \kappa is the wave number, \omega the angular frequency, \rho_a the fluid density, \mu and \lambda the Lamè parameters for which \mathbf{\sigma(u)} = 2 \mu \mathbf{\epsilon(u)} + \lambda tr(\mathbf{\epsilon(u)})I and \mathbf{\epsilon(u)} = \frac{1}{2} (\nabla \mathbf{u} + \nabla \mathbf{u} ^T)
The complete weak form should be
\int_{\Omega_a}(\nabla p \cdot \nabla \overline{q} - \kappa^2 p \overline{q})d \mathbf{x} + \int_{\Omega_s}(\mathbf{\sigma(u)} : \nabla \overline{\mathbf{v}} - \omega^2 \mathbf{u} \cdot \overline{\mathbf{v}})d \mathbf{x}
+ \int_{\Gamma_s} \rho_a \omega^2 (\mathbf{n}_1 \cdot \mathbf{u}) \overline{q} ds + \int_{\Gamma_s}(p \mathbf{n}_1) \cdot \overline{\mathbf{v}}ds = 0
Now in FEniCS:
1. Import libraries and define geometrical and material properties
from fenics import *
import matplotlib.pyplot as plt
x_1, y_1 = [-2, -1]
x_2, y_2 = [2, 1]
nx = 200
ny = 200
E, nu = Constant(1e5), Constant(0.)
mu = E/2./(1+nu)
lmbda = E*nu/(1+nu)/(1-2*nu)
rho = Constant(1.2)
omega = 1
k = 15
2. Create the mesh and the appropriate function spaces
mesh = RectangleMesh(Point(x_1, y_1), Point(x_2, y_2), nx, ny, "crossed")
V = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
Vv = VectorElement("Lagrange", mesh.ufl_cell(), 2)
W = FunctionSpace(mesh, V*Vv)
TRF = TrialFunction(W)
TTF = TestFunction(W)
(p, u) = split(TRF)
(q, v) = split(TTF)
3. Define a square subdomain for the solid obstacle and mark the domains and the boundaries
class Obstacle(SubDomain):
def inside(self, x, on_boundary):
return (between(x[1], (-0.75, 0.75)) and between(x[0], (0.5, 0.6)))
obstacle = Obstacle()
domains = MeshFunction("size_t", mesh, mesh.topology().dim(), 0)
domains.set_all(0)
obstacle.mark(domains, 1)
boundaries = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)
boundaries.set_all(0)
obstacle.mark(boundaries, 1)
4. Define boundary conditions: p |_{\partial B_{left}} = cos(\kappa x_1)
tol = 1E-14
def LEFT(x, on_boundary):
return on_boundary and abs(x[0] - x_1) < tol
bound_x = Expression(("cos(k*x[0])"), degree=1, k=k)
bcs = [DirichletBC(W.sub(0), Constant(0.0), "on_boundary"),
DirichletBC(W.sub(0), bound_x, LEFT)]
5. Define the variational problem, solve and plot
dX = Measure('dx', domain=mesh, subdomain_data=domains)
dS = Measure('ds', domain=mesh, subdomain_data=boundaries)
sigma = 2.0*mu*sym(grad(u)) + lmbda*tr(sym(grad(u)))*Identity(u.geometric_dimension())
n = FacetNormal(mesh)
#Fluid domain
g = Constant(0.0)
a_f = (inner(grad(p),grad(q)) - k**2 * p*q)*dX(0)
L_f = inner(g,q)*dS(0)
#Solid domain
T = Constant((0.0, 0.0))
a_s = (inner(sigma, grad(v)) - rho*omega**2*inner(u,v))*dX(1)
L_s = inner(T,v)*dS(1)
#Interface fluid-solid
a_i = (rho*omega**2 * inner(n('+'), u('+'))*q('+') + p*inner(n('+'), v('+')))*dS(1)
L_i = Constant(0.0)*dS(1)
#Weak form
a = a_f + a_s + a_i
L = L_f + L_s + L_i
#Solve
s = Function(W)
solve(a == L, s, bcs)
#Plot
plot(s)
plt.show()
This code produces the following error on solve(a == L, s, bcs)
ArityMismatch: Adding expressions with non-matching form arguments ('v_0',) vs ().
I believe the problem is in the #Interface fluid-solid definition, however I’d need someone more experienced to help me sort it out.
Many thanks for your time!
