# Calculating lift and drag from Navier stokes equation

Hello everyone,

I am new to using Fenics and also FEM/ fluid dynamics in general and I am struggling to complete code that will compute lift and drag around a cyliner inside a fluid flow.

So far this is my code, but I dont know how to start the computation of the lift and drag. Also I am not entireely sure if my code is working correctly, so tips on that would also be very welcome.

Thank you in advance for the help!

``````# -*- coding: utf-8 -*-
"""2D_lift_drag_NSE_clean.ipynb

Automatically generated by Colaboratory.

Original file is located at
"""

# Installation of FEniCS in Colab via
try:
import dolfin
except ImportError:
!wget "https://fem-on-colab.github.io/releases/fenics-install.sh" -O "/tmp/fenics-install.sh" && bash "/tmp/fenics-install.sh"
import dolfin

import csv
import pandas as pd
import numpy as np
from fenics import *
from mshr import *
import h5py

parameters["std_out_all_processes"] = False;

#DEFINING THE DOMAIN

L_domain  = 23.0 #lenght of the domain
H_domain  = 12.0 #height of the domain
X_disc  = 6.0 #position of the centre of the disc along the x-axis
Y_disc = H_domain/2
r_disc  = 2.0 #radius of the disc
resolution_disc = 200 #4*16 #number of fragrments

rectangle_domain   = Polygon([Point(0,0),Point(L_domain,0),Point(L_domain,H_domain),Point(0,H_domain)])
disc   = Circle(Point(X_disc,Y_disc),r_disc,resolution_disc)
domain = rectangle_domain - disc
mesh   = generate_mesh(domain,32) #integer value presents mesh refinment

mesh

# Define function spaces (P2-P1)
velocity_function_space = VectorFunctionSpace(mesh, "Lagrange", 2)  #velocity (2d vector)
pressure_function_space = FunctionSpace(mesh, "Lagrange", 1) #pressure (scalar)

u_trial = TrialFunction(velocity_function_space)
p_trial = TrialFunction(pressure_function_space)
v_test = TestFunction(velocity_function_space)
q_test = TestFunction(pressure_function_space)

# Define boundary markers
boundary_markers = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundary_markers.set_all(0)

class TopBottomBoundary(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and (near(x[1], 0, DOLFIN_EPS) or near(x[1], H_domain, DOLFIN_EPS))

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

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

class DiscBoundary(SubDomain):
def inside(self, x, on_boundary):
disc_x_left = X_disc - r_disc
disc_x_right = X_disc + r_disc
disc_y_bottom = Y_disc - r_disc
disc_y_top = Y_disc + r_disc
return on_boundary and (disc_x_left < x[0] < disc_x_right and disc_y_bottom < x[1] < disc_y_top)

top_bottom_boundary = TopBottomBoundary()
left_boundary = LeftBoundary()
outflow_boundary = OutflowBoundary()
disc_boundary = DiscBoundary()

top_bottom_boundary.mark(boundary_markers, 1)
left_boundary.mark(boundary_markers, 2)
outflow_boundary.mark(boundary_markers, 3)
disc_boundary.mark(boundary_markers, 4)

#DEFINE BOUNDARY CONDITIONS

#velocity boundary conditions
u_inflow = Constant((5.0, 0.0))
u_top_bottom = Constant((5.0, 0.0))
u_zero = Constant((0.0, 0.0))

top_bottom_boundary_conditions = DirichletBC(velocity_function_space, u_top_bottom, top_bottom_boundary)
left_boundary_conditions = DirichletBC(velocity_function_space, u_inflow, left_boundary)
disc_boundary_conditions = DirichletBC(velocity_function_space, u_zero, disc_boundary)
velocity_boundary_conditions  = [top_bottom_boundary_conditions, left_boundary_conditions, disc_boundary_conditions]

#pressure boundary conditions
p_outflow = Constant((0.0))

bco = DirichletBC(pressure_function_space, p_outflow, outflow_boundary)

#SIMULATION PARAMETERS
kinematic_viscosity = Constant(0.001) #Re =  (flow speed * charachteristic linear dimension) / kinematic viscosity
density = Constant(1000)
dynamic_viscosity = Constant(kinematic_viscosity * density)

total_simulation_time = 0.01
n_steps = 100
dt = total_simulation_time/n_steps
print("dt is: " +str(dt))

dt = Constant(dt)

#DEFINE SOLUTION FIELDS

u_prev = Function(velocity_function_space) #velocity in previous step
u_tent = Function(velocity_function_space) #velocity in tentative step
u_next = Function(velocity_function_space)  #velocity in next step
p_next = Function(pressure_function_space) #pressure in next time step

#WEAK FORM OF THE MOMENTUM EQUATION

#dx means we integrate over entire domain - inbuilt function
momentum_weak_form_residuum = (
(1/dt) * inner(u_trial - u_prev, v_test)*dx
+
+
)

momentum_weak_form_lhs = lhs(momentum_weak_form_residuum)
momentum_weak_form_rhs = rhs(momentum_weak_form_residuum)

#WEAK FROM OF THE PRESSURE POISSON PROBLEM

pressure_poisson_weak_form_rhs = -(1/dt) * div(u_tent)*q_test*dx

#WEAK FORM OF THE VELOCITY UPDATE EQUATION

velocity_update_weak_form_lhs = inner(u_trial, v_test)*dx
velocity_update_weak_form_rhs = (
inner(u_tent, v_test)*dx
-
)

# ASSEMBLE THE MATRICES
momentum_assembled_system_matrix = assemble(momentum_weak_form_lhs)
pressure_poisson_assembled_system_matrix = assemble(pressure_poisson_weak_form_lhs)
velocity_update_assembled_system_matrix = assemble(velocity_update_weak_form_lhs)

#RUN SIMULATION
# Use amg preconditioner if available
prec = "amg" if has_krylov_solver_preconditioner("amg") else "default"

for i in range(n_steps):

if i%10==0:
print ("N = " + str(i) + ", t=" +str(round(i*dt, 3)))

# COMPUTE TENTATIVE VELOCITY STEP
begin("Computing tentative velocity")

momentum_assembled_rhs = assemble(momentum_weak_form_rhs)

[bound_cond.apply(momentum_assembled_system_matrix, momentum_assembled_rhs) for bound_cond in velocity_boundary_conditions]

solve(momentum_assembled_system_matrix,
u_tent.vector(),
momentum_assembled_rhs,
"gmres",
"ilu"
)
#end()

# COMPUTE FOR PRESSURE
begin("Computing pressure correction")

pressure_poisson_assembled_rhs = assemble(pressure_poisson_weak_form_rhs)

bco.apply(pressure_poisson_assembled_system_matrix, pressure_poisson_assembled_rhs)

solve(pressure_poisson_assembled_system_matrix,
p_next.vector(),
pressure_poisson_assembled_rhs,
"cg",
prec
)
#end()

# COMPUTE VELOCITY CORRECTION
begin("Computing velocity correction")

velocity_update_assembled_rhs = assemble(velocity_update_weak_form_rhs)

[bound_cond.apply(velocity_update_assembled_system_matrix, velocity_update_assembled_rhs) for bound_cond in velocity_boundary_conditions]

solve(velocity_update_assembled_system_matrix,
u_next.vector(),
velocity_update_assembled_rhs,
"gmres",
"ilu"
)

# Plot solution
plot(p_next, title="Pressure", rescale=True)
plot(u_next, title="Velocity", rescale=True)

u_prev.assign(u_next)

``````

You can find a discussion on the computation of drag & lift here: How to Get Drag and Lift Coefficients for Flow Around a Cylinder