How to solve uniformed flow + dipole

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

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)
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 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.

Hello again!

Would love your thoughts on this. Thanks for taking a look!

I think I know the reason; it’s because all my equations are governed by Neumann boundary conditions, thus possessing infinitely many solutions. Therefore, it’s necessary to introduce additional conditions to fix the constant C.

However, I don’t want to impose too many constraints. Is it possible in FEniCS to specify fewer additional boundary conditions for PDEs with full Neumann boundary conditions? For example, could we only constrain the value at a single point within the domain, since the only task is to determine the constant C?

Best regards.

You can add a DirichletBC at a single point/dof. See for instance: Error in method = pointwise after definition of facets - #2 by Leo