How to add coupled body force per volume unit in Navier-Stokes eqts?


I am trying to include the force por volume unit in a Navier-Stokes problem caused by a concentration of ions in the fluid in an electric field.

First I am calculating the electric field as a vector field:

from fenics import *
import numpy as np

# Create mesh
channel = Rectangle(Point(0.0, 0.0), Point(5.0, 2.0))
mesh = generate_mesh(channel, 64)

# Define function spaces
scalar_order = 1
nodal_space = FunctionSpace(mesh, 'Lagrange', scalar_order)

# Define boundary conditions
left_side  = 'near(x[0], 0)'
right_side = 'near(x[0], 5)'

bcleft  = DirichletBC(nodal_space, Constant(100), left_side)
bcright  = DirichletBC(nodal_space, Constant(-100), right_side)

# Define trial and test function
L_i= TestFunction(nodal_space)
L_j= TrialFunction(nodal_space)

# Define constants, functions and expression used in variational problem. Normal value: 4.00701379*10^6 C/m^3 
rho_carga = Constant(0.0)
perm = Constant(8.8541*pow(10,-12))
phi = Function(nodal_space)
c = phi.vector()
A_ij = inner(grad(L_i), grad(L_j))*dx
b_ij = (rho_carga/perm)*L_j*dx

# Assemble potential matrices
A = assemble(A_ij)
b = assemble(b_ij)

#Apply boundary conditions
bcleft.apply(A, b)
bcright.apply(A, b)

# Solve potential problem

# Electric field
E = - grad(phi)

Then, using the navier-stokes formulation from the FEniCS Tutorial I would like to inlcude the coupling term f =q*Np*E with E being the electric vector field, q the charge of each ion (electron charge) and Np the concentration obtained with the diffusion-advection equation, also form the tutorial.

I haven’t found any examples where the term of external forces pero volume unit on the fluid isn’t Constant(0.0, 0.0)

Is there a way to do this?

P.S. I am new to FEniCS and this discourse site, sorry in advance for all my mistakes when posting.