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

ria01 · June 30, 2021, 3:06pm

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!

nate · June 30, 2021, 3:45pm

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

print(f"drag: {assemble(D)}")

ria01 · June 30, 2021, 3:52pm

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 ().

nate · June 30, 2021, 4:06pm

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.

ria01 · June 30, 2021, 4:21pm

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):
    return sym(nabla_grad(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)}")

ria01 · June 30, 2021, 5:49pm

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)

arash_goli · November 9, 2021, 7:43pm

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 :

name ‘fe’ is not defined, when running ds = fe.ds(subdomain_data=bc)
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?

arash_goli · November 9, 2021, 11:06pm

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()

Ben_Bird · December 6, 2021, 12:05pm

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

arash_goli · December 28, 2021, 11:46pm

I defined it. but i am still getting zero values.
Would you please share your entire code.

jcallaham · December 30, 2021, 2:16pm

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

Define the cylinder and other boundaries using SubDomain (rather than Expression as in the tutorial)
Define a MeshFunction and mark boundaries with integers (maybe helpful to name these, like CYLINDER=4)
Define BCs using the SubDomain objects
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.

arash_goli · December 30, 2021, 6:56pm

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

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
bc = MeshFunction(“size_t”, mesh, mesh.topology().dim() - 1, cylinder=4 ) and then delete this code
line: cylinder().mark(bc, 4)
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.
ds=Measure(‘ds’, domain=mesh, subdomain_data=bc (or 4 ?))

jcallaham · January 24, 2022, 12:46pm

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.

Topic		Replies	Views
How to obtain Drag and lift Coefficient for Flow Over a Cylinder - Navier Stokes, Benchmark 2D-2 (Re=100, periodic)	8	1659	November 20, 2021
Getting 0 drag coefficient and lift coefficient when working on the flow over a cylinder General	1	288	August 26, 2023
Calculating lift and drag from Navier stokes equation mathematics	1	321	September 25, 2023
The application of ds() and assemble() to calculate the drag and lift coefficients of the flow around a cylinder General	5	287	August 26, 2023
Problem with "Did you combine test or trial functions from different spaces?"	1	788	November 12, 2020

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

Related Topics