[Help] Compute lift of a NACA0012 using dolfin (new to FEniCS)

1

Hey there, I’ve recently taken these edx courses by KTH on Finite Element Modelling with FEniCS, and for the final project we need to develop some code ourselves. My partner and I were willing to compute the lift and drag of a NACA0012 with different angles of attack.

However, we checked if the lift was 0 at 0º of angle of attack (0012 is a symmetric airfoil, no lift at 0º), and it is not. It gives different results depending on the method used to solve (implicit Euler or midpoint rule) and the mesh resolution, but never zero, and always significantly larger than the drag.

I would be very grateful if you could have a look to my code and tell me if you notice any errors or if you just know why it does not work:

import math
import numpy as np

#Define angle of attack (degrees)
degrees = 0.;

#Define domain and mesh
XMIN = -0.5; 
XMAX = 2.; 
YMIN = 4.; 
YMAX = 6.; 
G = [XMIN, XMAX, YMIN, YMAX]; 
eps = 1e-5; 

mresolution = 100;

step = 0.01
points = []

for x in np.arange(0,1+step,step):
    x = x + eps
    y =  5 - 0.12/0.2 * (0.2969*sqrt(x) - 0.126*x - 0.3516*pow(x,2) + 0.2843*pow(x,3) - 0.1015*pow(x,4))
    point = Point(x,y)
    points.append(point)

for x in np.arange(1,0-step,-step):
    x = x + eps
    y =  5 + 0.12/0.2 * (0.2969*sqrt(x) - 0.126*x - 0.3516*pow(x,2) + 0.2843*pow(x,3) - 0.1015*pow(x,4))
    point = Point(x,y)
    points.append(point)


aoa = -degrees*math.pi/180
ox = 0.5; oy = 5.
rotated_points = []

for point in points:
    px = point.x()
    py = point.y()

    qx = ox + math.cos(aoa) * (px - ox) - math.sin(aoa) * (py - oy)
    qy = oy + math.sin(aoa) * (px - ox) + math.cos(aoa) * (py - oy)

    point = Point(qx,qy)
    rotated_points.append(point)


airfoil = Polygon(rotated_points)

mesh = generate_mesh(Rectangle(Point(G[0], G[2]), Point(G[1], G[3])) - airfoil, mresolution)



# Define function spaces and functions
VE = VectorElement("CG", mesh.ufl_cell(), 1); 
QE = FiniteElement("CG", mesh.ufl_cell(), 1); 
h = CellDiameter(mesh); # FEM functions
WE = MixedElement([VE,QE])
W = FunctionSpace(mesh, WE); 
V = FunctionSpace(mesh, VE); 
Q = FunctionSpace(mesh, QE) # Make a mixed space
(v, q) = TestFunctions(W); 
w = Function(W); 
(u, p) = (as_vector((w[0], w[1])), w[2]);
u0 = Function(V)

# Inflow velocity
#uin = Expression(("-96 + 40*x[1] - 4*x[1]*x[1]", "0.1"), YMAX=YMAX, element = V.ufl_element())
uin = Expression(("4.", "0"), element = V.ufl_element())

# Mark regions for boundary conditions 
om = Expression("x[0] > XMAX - eps ? 1. : 0.", XMAX=XMAX, eps=eps, element = Q.ufl_element()) 
im = Expression("x[0] < XMIN + eps ? 1. : 0.", XMIN=XMIN, eps=eps, element = Q.ufl_element())
nm = Expression("x[0] > XMIN + eps && x[0] < XMAX - eps ? 1. : 0.", XMIN=XMIN, XMAX=XMAX, eps=eps, element = Q.ufl_element())

# Initialization of variables needed for plotting
stepcounter = 0; pl = None; ax = None

# Timestep, viscosity and stabilization parameters
k = 0.1; nu = 1e-6; d = .2*h**(3./2.) 
# Time interval and penalty parameter
t, T = 0., 10.; gamma = 10*1./h 

# 0.5 - Midpoint rule,f 1.0 - Implicit Euler, 0.0 - Explicit Euler
theta = 1. 

# Time-stepping loop
while t < T: 
  um = theta*u + (1.0-theta)*u0 
  # Weak residual of stabilized FEM for Navier-Stokes eq.
  r = ((inner((u - u0)/k + grad(p) + grad(um)*um, v) + nu*inner(grad(um), grad(v)) + div(um)*q)*dx +
      gamma*(om*p*q + im*inner(u - uin, v) + nm*inner(u, v))*ds + # Weak boundary conditions
      d*(inner(grad(p) + grad(um)*um, grad(q) + grad(um)*v) + inner(div(um), div(v)))*dx) # Stabilization
  solve(r==0, w)  # Solve the Navier-Stokes PDE (one timestep)
  
  # Plot all quantities (see implementation above)
  pl, ax = plot_compact(u, t, stepcounter, Q, pl, ax) 
  
  # Shift to next timestep
  t += k; stepcounter += 1; u0 = project(u, V); 

# Drag force calculation 
psimarker = Expression("x[0] > 0-eps && x[0] < 1+eps && x[1] > -0.12/0.2*(0.2969*sqrt(x[0])-0.126*x[0]-0.3516*pow(x[0],2)+0.2843*pow(x[0],3)-0.1015*pow(x[0],4))+ 5-eps && x[1] < 0.12/0.2*(0.2969*sqrt(x[0])-0.126*x[0]-0.3516*pow(x[0],2)+0.2843*pow(x[0],3)-0.1015*pow(x[0],4))+5 + eps ? 1. : 0.", XMAX = XMAX, XMIN = XMIN, ox = ox, oy = oy, aoa = aoa, eps = eps, degree = 1) 
n = FacetNormal(mesh)
I = Identity(2)
def epsilon(z):
  return 0.5*(grad(z) + grad(z).T)
sigma = p*I - 2*nu*epsilon(u)
theta = Constant((1.0, 0.0))


M1 = psimarker*p*n[0]*ds # Drag (only pressure)
M2 = psimarker*p*n[1]*ds # Lift (only pressure)
#M3 = inner(g, u)*dx # Mean of the velocity in a region
M4 = psimarker*dot(dot(sigma, n), theta)*ds # Drag (full stress)
M5 = u[0]*dx # Mean of the x-velocity in the whole domain 

drag_force_pressure = assemble(M1)
print("Drag force (pressure) = ", drag_force_pressure)

drag_force_stress = assemble(M4)
print("Drag force (full stress) = ", drag_force_stress)

lift_force = assemble(M2)
print("Lift force = ",lift_force)

Have a nice day and thank you in advance!

Note: the imports and two functions are defined in another cell (I’m using Google Colab):

from google.colab import files
try:
    from dolfin import *; from mshr import *
except ImportError as e:
    !apt-get install -y -qq software-properties-common python-software-properties module-init-tools
    !add-apt-repository -y ppa:fenics-packages/fenics
    !apt-get update -qq
    !sed -e 's:artful:bionic:' /etc/apt/sources.list.d/fenics-packages-ubuntu-fenics-artful.list > temp
    !mv temp /etc/apt/sources.list.d/fenics-packages-ubuntu-fenics-artful.list
    !sed -e 's:artful:bionic:' /etc/apt/sources.list > temp
    !mv temp /etc/apt/sources.list
    !apt-get update -qq
    !apt install -y --no-install-recommends fenics    
    from dolfin import *; from mshr import *
    
import matplotlib.pyplot as plt
from IPython.display import clear_output, display, update_display
import time
import dolfin.common.plotting as fenicsplot 

def plot_compact(u, t, stepcounter, QQ, pl, ax, soltit = "Velocity"): # Compact plot utility function
  if stepcounter == 0:
    pl, ax = plt.subplots(); display(pl); clear_output(); # Plotting setup
  if stepcounter % 5 == 0:
    #u.set_allow_extrapolation(True)
    uEuclidnorm = project(sqrt(inner(u, u)), QQ)
    ax.cla()
    fig = plt.gcf()
    fig.set_size_inches(16, 4)
    plt.subplot(1, 2, 1)
    pp = plot(uEuclidnorm, cmap="coolwarm")
    plt.title("%s at t=%f" % (soltit, t)) # Plot norm of solution
    if t == 0:
      plt.axis(G)
      plt.colorbar(pp, shrink=0.5)
    plt.subplot(1, 2, 2);
    if t == 0:
      plot(QQ.mesh())
      plt.title("Mesh") # Plot mesh
    plt.tight_layout(); dpl = display(pl, display_id="test");
  
  return (pl, ax)

def plot_compact_static(u, QQ, label = "Velocity"): # Compact plot utility function
  pl, ax = plt.subplots();
  uEuclidnorm = project(sqrt(inner(u, u)), QQ); ax.cla(); fig = plt.gcf(); fig.set_size_inches(16, 4)
  plt.subplot(1, 2, 1); pp = plot(uEuclidnorm, cmap="coolwarm"); plt.title("%s" % (label)) # Plot norm of velocity
  plt.axis(G); plt.colorbar(pp, shrink=0.5); 
  plt.subplot(1, 2, 2);
  plot(QQ.mesh()); plt.title("Mesh") # Plot mesh
  plt.tight_layout(); dpl = display(pl, display_id="test");
  
  return (pl, ax)
2

First step is to look at the solution of the problem graphically (for instance using Paraview), does it look like you are solving the correct problem?

One suggestion is to move the residual out of the while loop and assign u0 to u using the FunctionAssigner, see for instance this test for intended usage. This would speed up your computations quite drastically. I would also suggest using the constant wrappers for variables such as theta, nu and k.

I would also use u,p= split(w) instead of

I would suggest that you start by verifying your variational formulation by using the method of manufactured solutions, for instance solving the Taylor-Green vortex problem, verifying that you obtain the expected spatial convergence rates.

My next step would be to solve the DFG benchmark problem, which can be used to compare drag and lift coefficients.

Solving the Taylor-Green problem will help you verify the implementation of the navier stokes variational problem, while solving the DFG benchmark will help you verify your lift computation.

1 Like
3

Okay thank you, I’ll have a look to all of these!