Hello, as a newcomer to FEniCS, I am attempting to solve a simple Laplace problem with a known analytical solution to verify the results obtained with FEniCS.
Consider a problem involving a uniform flow with a dipole. The uniform flow V is parallel to the x-axis, moving from left to right, and a dipole aligned along the negative x-axis is placed at the origin.
In this case, the potential function ϕ satisfies Δϕ=0.
An artificial far field is set where, at the outer circle radius, the velocity in the x-direction is 1 and in the y-direction is 0, thus ∂ϕ/∂x=1 and ∂ϕ/∂y=0.
At the same time, at the inner circle radius (determined by the dipole), a natural Neumann boundary condition is satisfied where grad(ϕ)⋅n=0.
I cannot think of any other boundary conditions.
In short, Δϕ=0
∂ϕ/∂x=1 and ∂ϕ/∂y=0 for R_out
grad(ϕ)⋅n=0 for R_in
My FEniCS code is:
from dolfin import *
import matplotlib.pyplot as plt
file = 'dipole_flow_uns'
mesh = Mesh(file + ".xml")
boundaries = MeshFunction("size_t", mesh, file + "_facet_region.xml")
# to prove the right Boundary, for R_out (2) in blue and R_in (3) in red
plt.figure()
plot(mesh, color='black')
boundary_vertices = boundaries.array()
mesh_coords = mesh.coordinates()
mesh_topo = mesh.topology()
colors = {2: 'blue', 3: 'red'}
for facet in facets(mesh):
if boundary_vertices[facet.index()] in colors:
vertices = facet.entities(0)
x = [mesh_coords[vertex][0] for vertex in vertices]
y = [mesh_coords[vertex][1] for vertex in vertices]
plt.plot(x, y, color=colors[boundary_vertices[facet.index()]], linewidth=2)
plt.show()
P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
W = FunctionSpace(mesh, P1)
u = TrialFunction(W)
v = TestFunction(W)
a = inner(grad(u), grad(v)) * dx
L = Constant(0) * v * dx # Right-hand side of Laplace equation
n = FacetNormal(mesh)
g_expression = Constant(1) # Neumann BC ∂phi/∂x = 1 for R_out
h_expression = Constant(0) # Neumann BC ∂phi/∂y = 0 for R_out
neumann_term = g_expression * v * dot(n, Constant((1, 0))) * ds(2) + h_expression * v * dot(n, Constant((0, 1))) * ds(2)
L += neumann_term
u_sol = Function(W)
solve(a == L, u_sol)
print("Solution range:", u_sol.vector().min(), u_sol.vector().max())
plot(u_sol, title="Solution u")
plt.show()
I have imported the files “dipole_flow_uns.xml” and “dipole_flow_uns_facet_region.xml” converted from .msh format. The download link is: https://drive.google.com/file/d/1AtehifKiD2Hh2_9HNldPgN2QRB_oynXO/view?usp=drive_link and
https://drive.google.com/file/d/1SmFaz18X3YxqcS0_sJX1CCCT_eVh4Kfl/view?usp=drive_link
I have set the marker for R_out as 2 and R_in as 3, and this has been confirmed through plotting, as shown in the figure below:
Subsequently, my FEniCS code is quite basic. Since the Neumann boundary condition for the inner circle is naturally satisfied, I have only applied the Neumann boundary condition for the outer circle, but I am not sure if it is correct.
Nevertheless, the result of the solution is completely unsuccessful, with a solution range of 0.0 to 0.0.
This means that the value throughout the computation domain is 0.0, which is absurd. Even though having only Neumann boundary conditions might involve differences due to a constant C, it shouldn’t be like this.
Given that this problem is not complex and there is an analytical solution to verify the results from FEniCS, how should I modify my FEniCS code?
Any help would be greatly appreciated.
Best regards.