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
https://colab.research.google.com/drive/1m4XEwJ-wz1Fmd0SztggoOemI0MQaTYia
"""
# 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
from google.colab import drive
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
+
inner(grad(u_prev) * u_prev, v_test)*dx
+
kinematic_viscosity * inner(grad(u_trial), grad(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_lhs = inner(grad(p_trial), grad(q_test))*dx
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
-
dt * inner(grad(p_next), 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)