How to plot the pressure distribution on the airfoil surface

1

Dear all, I am new to FENICS and want to get some help from you. I tried to use FENICS to model airflows around airfoil at Re=20. After computation, I want to get the pressure distribution on the airfoil surface, but have no idea how to do that, could you give me some clues

def build_space(airfoil, Rect, N_bulk, u_in):
    """Prepare data for DGF benchmark. Return function
    space, list of boundary conditions and surface measure
    on the airfoil."""
    Xmin, Xmax, Ymin, Ymax = Rect[0], Rect[1], Rect[2], Rect[3]
    bc_pt = [Point(x,y) for x,y in airfoil]
    airfoil = Polygon(bc_pt)
    rect = Rectangle(Point(Xmin, Ymin), Point(Xmax, Ymax))
    geometry = rect - airfoil
    mesh = generate_mesh(geometry, N_bulk)
    #plot(mesh, title='Mesh of the total domain')
    #plt.show()

    EPS = 1e-4
    # Construct facet markers
    class Airfoil(SubDomain):
        def inside(self, x, on_boundary):
            return on_boundary and x[0] > -EPS and  x[0] < 1.0+EPS and x[1] > -0.1+EPS and x[1] < 0.15-EPS
    class Inlet(SubDomain):
        def inside(self, x, on_boundary):
            return near(x[0], Xmin) and on_boundary
    class Outlet(SubDomain):
        def inside(self, x, on_boundary):
            return near(x[0], Xmax) and on_boundary
    class Wall(SubDomain):
        def inside(self, x, on_boundary):
            return  (near(x[1], Ymin) or near(x[1], Ymax) ) and on_boundary
    sub_domains = MeshFunction("size_t", mesh, mesh.topology().dim()-1)
    sub_domains.set_all(0) # label all facets as 0
    inlet = Inlet()
    inlet.mark(sub_domains, 1)
    outlet = Outlet()
    outlet.mark(sub_domains, 2)
    wall = Wall()
    wall.mark(sub_domains, 3)
    air_foil = Airfoil()
    air_foil.mark(sub_domains, 4)
        # Build function spaces (Taylor-Hood)
    P2 = VectorElement("P", mesh.ufl_cell(), 2)
    P1 = FiniteElement("P", mesh.ufl_cell(), 1)
    TH = MixedElement([P2, P1])
    # P2 = VectorElement("Lagrange", triangle, 2)
    # P1 = FiniteElement("Lagrange", triangle, 1)
    # TH = P2 * P1
    W = FunctionSpace(mesh, TH)

    # Prepare Dirichlet boundary conditions
    bc_walls = DirichletBC(W.sub(0), (0, 0), sub_domains, 3)
    bc_airfoil = DirichletBC(W.sub(0), (0, 0), sub_domains, 4)
    bc_in = DirichletBC(W.sub(0), u_in, sub_domains, 1)
    bc_out = DirichletBC(W.sub(1), 0, sub_domains, 1)
    bcs = [bc_airfoil, bc_walls, bc_in]


    # Prepare surface measure on airfoil
    ds_airfoil= Measure("ds", subdomain_data=sub_domains, subdomain_id=4)

    return W, bcs, ds_airfoil


def solve_navier_stokes(W, nu, bcs):
    """Solve steady Navier-Stokes and return the solution"""
    # Define variational forms
    v, q = TestFunctions(W)
    w = Function(W)
    u, p = split(w)
    F = nu*inner(grad(u), grad(v))*dx + dot(dot(grad(u), u), v)*dx - p*div(v)*dx - q*div(u)*dx
    # Solve the problem
    solve(F == 0, w, bcs, solver_parameters = {'newton_solver': {'linear_solver': 'mumps'}})
    # update parameters
    return w
2

Hi jhssyb!

Take a look at the end of the Stokes equations with an iterative solver demo in the FEniCS docs for a similar problem case. The last section demonstrates the use of the split() method that allows you to separate the velocity and pressure functions contained in w.

I would use a code like this for saving .pvd files for visualization with Paraview:

# setup and solve problem
u,p = w.split()
File("u.pvd") << u
File("p.pvd") << p

Hope that helps!

1 Like
3

Thanks for your suggestion, Grantn. I am sorry as I did not explain it clearly. My computation domain is a rectangle with an embedded airfoil. After computation, I want to plot the pressure along the 2D airfoil curve. Even if I can get the pressure distribution of the whole domain, I cannot extract pressure information along the airfoil boundary. Could you give me some clues?

4

Here is a minimal example of how to extract the the pressure from a boundary surface:

from dolfin import *

mesh = UnitSquareMesh(10, 10)


class Boundary(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 0) and on_boundary


el1 = VectorElement("Lagrange", mesh.ufl_cell(), 2)
el2 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)

W = FunctionSpace(mesh, MixedElement([el1, el2]))

w = Function(W)

# Assign some data to mixed space
p = Function(W.sub(1).collapse())
p.interpolate(Expression("2*x[1]", degree=1))
fa = FunctionAssigner(W.sub(1), p.function_space())
fa.assign(w.sub(1), p)

# Mark all vertices on boundary with 1
vf = MeshFunction("size_t", mesh, 0)
Boundary().mark(vf, 1)

# Get sub-component
w_1 = w.sub(1, deepcopy=True)
v_to_d = vertex_to_dof_map(w_1.function_space())
boundary_vertices = vf.where_equal(1)
vertex_coords = mesh.coordinates()[boundary_vertices]
for (coord, entity) in zip(vertex_coords, boundary_vertices):
    print(coord, w_1.vector().get_local()[v_to_d[entity]])

yielding:

[ 0.  0.] 0.0
[ 0.   0.1] 0.2
[ 0.   0.2] 0.4
[ 0.   0.3] 0.6
[ 0.   0.4] 0.8
[ 0.   0.5] 1.0
[ 0.   0.6] 1.2
[ 0.   0.7] 1.4
[ 0.   0.8] 1.6
[ 0.   0.9] 1.8
[ 0.  1.] 2.0

which you can append to an array an plot with matplotlib.

5

Many thanks for your help, Dokken. The minimal example is extremely useful. Just a small question. Since w_1 is a scalar, do we have to use w_1.vector() in the last line?

6

The vector is the vector containing the values at all degrees of freedom, not a vector in the (x,y,z) coordinate sense

7

Thanks for your explanation, great help