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)}")
1 Like

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):
    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)}")

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)
1 Like