# How to Get Drag and Lift Coefficients for Flow Around a Cylinder

So far for my problem I have been successful in getting the Umax and the inflow\outflow flux. I am using Solving PDEs in Python - <br> The FEniCS Tutorial Volume I as a framework to solve the problem. My last step is to obtain the drag and lift values, and as of now printing them in this way seen below gives the following result.

Input:

``````force = dot(sigma(U, p_), n)
D = (force[0]/0.002)*ds(4)
L = (force[1]/0.002)*ds(4)
print('drag:', D)
print('lift:', L)
``````

Output:

``````u max: 2.1079706751320226
p max: 2.99549776658916
inflow flux: -0.4100347906645285
outflow flux: 0.40923874357201123
drag: { ((({ A | A_{i_{130}, i_{131}} = -1 * ({ A | A_{i_{128}, i_{129}} = I[i_{128}, i_{129}] * f_39 })[i_{130}, i_{131}] }) + ({ A | A_{i_{126}, i_{127}} = (sym(nabla_grad({ A | A_{i_8} = 0.5 * (v_1 + f_30)[i_8] })))[i_{126}, i_{127}] * 2 * f_45 })) . (n))[0] / 0.002 } * ds(<Mesh #7>[4], {})
lift: { ((({ A | A_{i_{130}, i_{131}} = -1 * ({ A | A_{i_{128}, i_{129}} = I[i_{128}, i_{129}] * f_39 })[i_{130}, i_{131}] }) + ({ A | A_{i_{126}, i_{127}} = (sym(nabla_grad({ A | A_{i_8} = 0.5 * (v_1 + f_30)[i_8] })))[i_{126}, i_{127}] * 2 * f_45 })) . (n))[1] / 0.002 } * ds(<Mesh #7>[4], {})
``````

So, my question is how do I convert this Form() object to a scalar value? I appreciate any help!

Those forms are integrals, which you can compute using DOLFIN’s `assemble()` function. Something like

``````print(f"drag: {assemble(D)}")
``````
2 Likes

Thank you for your assistance, however I am getting an error using this function. It reads the following …

``````Traceback (most recent call last):
File "main.py", line 154, in <module>
print(f"drag: {assemble(D)}")
File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/assembling.py", line 202, in assemble
dolfin_form = _create_dolfin_form(form, form_compiler_parameters)
File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/assembling.py", line 60, in _create_dolfin_form
return Form(form,
File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/form.py", line 43, in __init__
ufc_form = ffc_jit(form, form_compiler_parameters=form_compiler_parameters,
File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/jit/jit.py", line 50, in mpi_jit
return local_jit(*args, **kwargs)
File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/jit/jit.py", line 100, in ffc_jit
return ffc.jit(ufl_form, parameters=p)
File "/usr/lib/python3/dist-packages/ffc/jitcompiler.py", line 217, in jit
module = jit_build(ufl_object, module_name, parameters)
File "/usr/lib/python3/dist-packages/ffc/jitcompiler.py", line 130, in jit_build
module, signature = dijitso.jit(jitable=ufl_object,
File "/usr/lib/python3/dist-packages/dijitso/jit.py", line 165, in jit
header, source, dependencies = generate(jitable, name, signature, params["generator"])
File "/usr/lib/python3/dist-packages/ffc/jitcompiler.py", line 65, in jit_generate
code_h, code_c, dependent_ufl_objects = compile_object(ufl_object,
File "/usr/lib/python3/dist-packages/ffc/compiler.py", line 142, in compile_form
return compile_ufl_objects(forms, "form", object_names,
File "/usr/lib/python3/dist-packages/ffc/compiler.py", line 185, in compile_ufl_objects
analysis = analyze_ufl_objects(ufl_objects, kind, parameters)
File "/usr/lib/python3/dist-packages/ffc/analysis.py", line 89, in analyze_ufl_objects
form_datas = tuple(_analyze_form(form, parameters)
File "/usr/lib/python3/dist-packages/ffc/analysis.py", line 89, in <genexpr>
form_datas = tuple(_analyze_form(form, parameters)
File "/usr/lib/python3/dist-packages/ffc/analysis.py", line 169, in _analyze_form
form_data = compute_form_data(form,
File "/usr/lib/python3/dist-packages/ufl/algorithms/compute_form_data.py", line 407, in compute_form_data
check_form_arity(preprocessed_form, self.original_form.arguments(), complex_mode)  # Currently testing how fast this is
File "/usr/lib/python3/dist-packages/ufl/algorithms/check_arities.py", line 177, in check_form_arity
check_integrand_arity(itg.integrand(), arguments, complex_mode)
File "/usr/lib/python3/dist-packages/ufl/algorithms/check_arities.py", line 159, in check_integrand_arity
arg_tuples = map_expr_dag(rules, expr, compress=False)
File "/usr/lib/python3/dist-packages/ufl/corealg/map_dag.py", line 26, in map_expr_dag
result, = map_expr_dags(function, [expression], compress=compress)
File "/usr/lib/python3/dist-packages/ufl/corealg/map_dag.py", line 75, in map_expr_dags
r = handlers[v._ufl_typecode_](v, *[vcache[u] for u in v.ufl_operands])
File "/usr/lib/python3/dist-packages/ufl/algorithms/check_arities.py", line 48, in sum
raise ArityMismatch("Adding expressions with non-matching form arguments {0} vs {1}.".format(_afmt(a), _afmt(b)))
ufl.algorithms.check_arities.ArityMismatch: Adding expressions with non-matching form arguments ('v_1',) vs ().``````

Make sure you’re computing the lift and drag of the solution function. Not a trial or test function. It’s difficult to debug without a MWE.

1 Like

I am not using the test or trial to compute the values. I have a reference to the code I am based my solution off of in the original post, but I am unable to pinpoint the issue if I can get assistance that would be wonderful. I am posting my MWE now below, thank you!

``````# Function Space
V = VectorFunctionSpace(mesh, 'CG', 2)
Q = FunctionSpace(mesh, 'CG', 1)

# Boundary Conditions
class inflow(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 0) and on_boundary
class outflow(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 2.2) and on_boundary
class walls(SubDomain):
def inside(self, x, on_boundary):
return (near(x[1], 0) or near(x[1], 0.41)) and on_boundary
class cylinder(SubDomain):
def inside(self, x, on_boundary):
return x[0]>=0.15 and x[0]<=0.25 and x[1]>=0.15 and x[1]<=0.25 and on_boundary

inflow_profile = ('4.0*1.5*x[1]*(0.41 - x[1]) / pow(0.41, 2)', '0')

bc = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
bc.set_all(0)
inflow().mark(bc, 1)
outflow().mark(bc, 2)
walls().mark(bc, 3)
cylinder().mark(bc, 4)

bcu_inflow = DirichletBC(V, Expression(inflow_profile, degree=2), inflow())
bcp_outflow = DirichletBC(Q, Constant(0), outflow())
bcu_walls = DirichletBC(V, Constant((0, 0)), walls())
bcu_cylinder = DirichletBC(V, Constant((0, 0)), cylinder())

bcu = [bcu_inflow, bcu_walls, bcu_cylinder]
bcp = [bcp_outflow]
ds = fe.ds(subdomain_data=bc)

# Trial/Test Functions
u = TrialFunction(V)
v = TestFunction(V)
p = TrialFunction(Q)
q = TestFunction(Q)

# Time-Step Function Solutions
u_n = Function(V)
u_  = Function(V)
p_n = Function(Q)
p_  = Function(Q)

# Functions
def epsilon(u):
def sigma(u, p):
return 2*mu*epsilon(u) - p*Identity(len(u))

... (variational problem)

t = 0
for j in range(num_steps):
# Update current time
t += dt

# Step 1: Tentative velocity step
b1 = assemble(L1)
[bc.apply(b1) for bc in bcu]
solve(A1, u_.vector(), b1, 'bicgstab', 'hypre_amg')

# Step 2: Pressure correction step
b2 = assemble(L2)
[bc.apply(b2) for bc in bcp]
solve(A2, p_.vector(), b2, 'bicgstab', 'hypre_amg')

# Step 3: Velocity correction step
b3 = assemble(L3)
solve(A3, u_.vector(), b3, 'cg', 'sor')

# Update previous solution
u_n.assign(u_)
p_n.assign(p_)

force = dot(sigma(U, p_n), n)
D = (force[0]/0.002)*ds(4)
L = (force[1]/0.002)*ds(4)
print(f"drag: {assemble(D)}")
print(f"lift: {assemble(L)}")``````

I figured out the issue, it was due to not defining the force correctly, I changed it to be the code below and it worked. Thank you for the assistance!

``````I = Identity(u_.geometric_dimension())
force = dot(-p_*I + 2.0*nu*sym(grad(u_)), n)``````
3 Likes

Hello Ria

I am using Solving PDEs in Python -
The FEniCS Tutorial Volume I
as a framework to solve the problem as well.
I am trying to obtain the drag and lift values but as you have done, but I end up with the following errors :

1. name ‘fe’ is not defined, when running ds = fe.ds(subdomain_data=bc)
2. Dot product requires non-scalar arguments, got arguments with ranks 2 and 0 , when running
force = dot(-p_I + 2.0nu*sym(grad(u_)), n)
Would you please let me know how did you fix it?

I figured them out as follows :
ds = ds(subdomain_data=bc)
And I changed n by j in my for loop.

But why are you using ds(4)? I get zero values for drag and lift when ds(4), However, different values when ds()

Hi there,

Solution to problem 2):

You need to define

n = FacetNormal(mesh)

This should solve your problem and also give you non-zero values for drag and lift when you use ds(4).
Keep j in your for loop so there are no conflicting variable names.

Ben

1 Like

I defined it. but i am still getting zero values.

The `ds(4)` is indicating that the integral should be over the subdomain labeled “4”. The integer labels are arbitrary, but the cylinder is marked as “4” in `cylinder().mark(bc, 4)`. But if you just use `ds()` I think it will integrate over all boundaries.

So the boundaries are handled slightly different from the tutorial, but I think the full steps (which are all in the code from ria01) should be

1. Define the cylinder and other boundaries using `SubDomain` (rather than `Expression` as in the tutorial)
2. Define a `MeshFunction` and mark boundaries with integers (maybe helpful to name these, like `CYLINDER=4`)
3. Define BCs using the `SubDomain` objects
4. Define the measure `ds` using your `MeshFunction` as the `subdomain_data`

Then it’s just the form and assembly, which it sounds like you have figured out already.

1 Like

Would you please check the following steps if I did correctly :

1. class cylinder(SubDomain):
def inside(self, x, on_boundary):
return x[0]>=0.15 and x[0]<=0.25 and x[1]>=0.15 and x[1]<=0.25 and on_boundary

2. bc = MeshFunction(“size_t”, mesh, mesh.topology().dim() - 1, cylinder=4 ) and then delete this code
line: cylinder().mark(bc, 4)

3. As the same ria01 code : # Boundary Conditions
class inflow(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 0) and on_boundary
and so the others.

4. ds=Measure(‘ds’, domain=mesh, subdomain_data=bc (or 4 ?))

Looks okay at a glance, but it’d be hard to say if it’s correct or not without running the full code and comparing against something. Usually the best way is to validate against known benchmarks. For cylinder flow this paper has a nice summary towards the end for various Reynolds numbers (tables 2 and 3).

Obviously if you do try to match results like this be careful with domain size, boundary conditions, non-dimensionalization, and so on. And check with different mesh resolution, time steps, and all of that standard CFD stuff.