Impose electric current magnitude / Floating potential

variational formulation
1

Dear all,

I’m working on a static electric conduction problem, governed by a simple Poisson equation:

-\nabla \cdot [\sigma(x,y)\nabla u(x,y)] = 0

with \sigma the electric conductivity and u the electric potential.

I understand that if I impose Dirichlet Boundary Conditions (BC) on a boundary, the electric potential u will be fixed on this boundary, while a Neumann BC will impose \frac{\delta u}{\delta n} on the boundary (with n the normal vector to the boundary).

For a known conductivity distribution \sigma(x,y), the Neumann boundary conditions could somehow be linked to imposing the normal component of the current density at the boundary. With the general current density distribution being defined as: J(x,y) = -\sigma (x,y) \nabla u.

However, I would like to impose the current magnitude I_0 and not the current density on a boundary \Gamma_1. Do you have some recommendations on how to achieve that? I saw in some commercial FEA software that they impose:

U\equiv const \\ \int_{\Gamma_1} -n\cdot J \; dS = I_0

Do you have some recommendations/tips on how to achieve that in FEniCS?

Thank you for your time,
Best
Adrien

2

In weak form, the problem is find u \in V such that

\int_{\Omega} \sigma \nabla u \cdot \nabla v \, d\Omega - \int_{\Gamma} \sigma \frac{\partial u}{\partial n} v \, d\Gamma= 0 \quad \forall v \in V \text{ .}

Using the definition you gave for the current density,

- \int_{\Gamma} \sigma \frac{\partial u}{\partial n} \, d\Gamma = \int_{\Gamma} \mathbf{J} \cdot \mathbf{n} \, d\Gamma = - I_{0} \text{ .}

Therefore, it may be easiest to work back from I_{0} (I am assuming you have that since you mentioned that you wanted to impose it) and calculate the \frac{\partial u}{\partial n} for

\int_{\Gamma} \sigma \frac{\partial u}{\partial n} \, d\Gamma = I_{0} \text{ .}

This can then be inserted into the weak form and written up in FEniCS using UFL, as usual: see, for example, the demo Poisson equation with Neumann boundary conditions.

1 Like
3

Thank you for your answer.

The link you shared didn’t work, but I think I’ve found a similar one: Poisson equation with pure Neumann boundary conditions
Is there a more detailed development of the construction of the weak form for this problem?
I don’t understand what are the mathematical steps described by the text below:

However, the variational problem expressed below is well-posed as the Lagrange multiplier introduced to satisfy the condition \int_\Omega u \, dx=0 effectively redefines the right-hand side such that it safisfies the necessary condition \int_\Omega f \, dx = -\int_{\partial\Omega} g \, ds

Nevertheless, I don’t really get what you’re suggesting. How can I add this condition into the weak form of the problem? Is there a suggested workflow to enforce custom conditions on the boundaries?

Furthermore, there is another condition to enforce. Due to the “floating potential” constraint, we should also have:

U \equiv const \Rightarrow \nabla u \times n = 0 \quad {\text{on } \Gamma_1}.
4

Not sure why that link stopped working. This demo shows the application of a Neumann condition for a Poisson equation.

My point was that you implied you knew what I_0 was; therefore, you can use that I_0 to find \frac{\partial u}{\partial n} on the boundary since

\int_{\Gamma} \sigma \frac{\partial u}{\partial n} \, d\Gamma = I_{0} \text{ .}

Now you know \frac{\partial u}{\partial n} on the boundary, you can insert that into the weak form:

\int_{\Omega} \sigma \nabla u \cdot \nabla v \, d\Omega - \int_{\Gamma} \sigma \frac{\partial u}{\partial n} v \, d\Gamma= 0 \quad \forall v \in V \text{,}

which is how you implement a Neumann boundary condition.

In regards to having \nabla u \times \mathbf{n} = 0 on \Gamma, I am not too sure I understand how that would fully define a unique solution. In general, it is true that a Poisson equation with pure Neumann conditions only defines a solution up to a constant, which explains the need for an additional constraint in the example Poisson with pure Neumann conditions.

5

I think the problem only really makes sense if \Gamma_1 is a strict subset of \partial\Omega. Otherwise, only I_0 = 0 would be compatible with a zero source term, and the solution would just be u \equiv \text{const} in the entire domain. I hacked together an ad hoc Nitsche-like formulation for the problem on a unit square, with a nonzero Dirichlet BC on the left, insulating (J=0) Neumann BCs on the top and bottom, and this special BC on the right, in terms of an average current density,

J_\text{avg} := -I_0 / \vert\Gamma_1\vert\text{ ,}

where \vert\Gamma_1\vert is the length of the right edge. The formulation is: Find (u,u_0)\in S\times\mathbb{R} such that, \forall(v,v_0)\in V\times\mathbb{R},

\int_\Omega\sigma\nabla u\cdot\nabla v\,d\Omega\\ - \int_{\Gamma_1}v\sigma\nabla u\cdot\mathbf{n}\,d\Gamma + \int_{\Gamma_1}(u-u_0)\sigma\nabla v\cdot\mathbf{n}\,d\Gamma \\ + \int_{\Gamma_1}(J_\text{avg} + \sigma\nabla u\cdot\mathbf{n})v_0\,d\Gamma = 0\text{ ,}

where u_0 is the free constant that u must equal on \Gamma_1, S is the trial space satisfying the Dirichlet BC on the left, and V is the corresponding test space. The solution in the following FEniCS implementation appears to converge toward satisfaction of both u = u_0 on \Gamma_1 and \int_{\Gamma_1}\mathbf{J}\cdot\mathbf{n}\,d\Gamma = -I_0, as the mesh is refined:

from dolfin import *
N = 64
mesh = UnitSquareMesh(N,N)
cell = mesh.ufl_cell()
uE = FiniteElement("CG",cell,1)
u0E = FiniteElement("Real",cell,0)
W = FunctionSpace(mesh,MixedElement([uE,u0E]))
x = SpatialCoordinate(mesh)

# u is the potential, while u0 is a global scalar unknown, which is the
# constant potential over the right face of the domain.
u,u0 = TrialFunctions(W)
v,v0 = TestFunctions(W)

# Apply some arbitrary Dirichlet BC to the potential on the left face:
bc = DirichletBC(W.sub(0),
                 Expression("sin(3*pi*x[1])",degree=2),
                 "x[0]<DOLFIN_EPS")

# Characteristic function for right face:
right_char = conditional(gt(x[0],Constant(1.0-DOLFIN_EPS)),1.0,Constant(0))

# Formulation, using a Nitsche-like enforcement of u=u0 on the right
# face of the domain, and weak enforcement of an average current.
sigma = Constant(1)
def J(u):
    return -sigma*grad(u)
n = FacetNormal(mesh)
I_0 = 1.0
meas_Gamma1 = assemble(right_char*ds)
J_avg = Constant(-I_0/meas_Gamma1)
res = inner(-J(u),grad(v))*dx \
    + right_char*(dot(J(u),n)*v - dot(J(v),n)*(u-u0)
                  + (J_avg - dot(J(u),n))*v0)*ds
w = Function(W)
solve(lhs(res)==rhs(res),w,[bc,])
u_out,u0_out = w.split()

# Check that u=u0 on the right boundary:
e = u_out-u0_out
print(sqrt(assemble(right_char*e*e*ds)))

# Check that the total current through the right boundary
# has average value J_avg:
print(abs(assemble(right_char*(dot(J(u_out),n) - J_avg)*ds)))

# Visualize output:
u_out.rename("u","u")
File("u.pvd") << u_out
2 Likes
6

Thank you very much @nav and @kamensky for your answers.
It is true that I didn’t specify all the boundary conditions (BC) of my problem, sorry for that.

\partial\Omega = \Gamma_D \,\cap \, \Gamma_{0} \,\cap \, \Gamma_{1}

with \Gamma_D the boundaries onto which Dirichlet BC are applied, \Gamma_{0} the boundaries onto which null Neumann BC are imposed, and \Gamma_{1} the boundary onto which I desire to enforce a floating potential u\equiv const, as well as a constant electric current I_0 = \int_{\Gamma_1} \sigma \frac{\partial u}{\partial n} dS.

The code appears to do exactly that, thanks a lot!

7

Hai,

I wanted to understand the implementation of the floating potential condition in fenics. As per the above code, I notice the value of current is chosen to be 1. But what in case if we don’t know the value of current?.
Can you explain the implementation of the floating potential in more detail, which will be really helpful.

Regards,
Gunasheela

8

I think the easiest way to understand the formulation is to start with the simpler problem of assuming that u_0 is a known constant on \Gamma_1. Then the formulation reduces to: Find u\in S such that, for all v\in V,

\int_\Omega\sigma\nabla u \cdot\nabla v\,d\Omega - \int_{\Gamma_1}v\sigma\nabla u\cdot\mathbf{n}\,d\Gamma + \int_{\Gamma_1}(u - u_0)\sigma\nabla v\cdot\mathbf{n}\,d\Gamma = 0\text{ .}

This is just the non-symmetric variant of Nistche’s method for enforcing the Dirichlet condition u = u_0 on \Gamma_1. The first and second terms come from integrating \int_\Omega-\nabla\cdot(\sigma\nabla u)v\,d\Omega by parts without assuming that v goes to zero on \Gamma_1. The third term is clearly consistent with the Dirichlet condition, since it will be zero when u=u_0, and its contribution to the formulation’s left-hand-side bilinear form would cancel that of the second term when plugging in the same function into the test and trial arguments, to test for stability.

If we further assume that u_0\in \mathbb{R} is unknown, we need another equation to have a non-singular system. That is where the third line of my formulation comes in. If the formulation is satisfied \forall (v,v_0)\in V\times\mathbb{R}, then it is satisfied for the specific choice of (v,v_0) = (0,1). Plugging in that choice, the formulation reduces to

\int_{\Gamma_1}(J_\text{avg} + \sigma\nabla u\cdot\mathbf{n})\,d\Gamma = 0\quad\Rightarrow\quad -I_0 + \int_{\Gamma_1}\sigma\nabla u\cdot\mathbf{n}\,d\Gamma = 0\quad\Rightarrow\quad\int_{\Gamma_1}\mathbf{J}\cdot\mathbf{n}\,d\Gamma = -I_0\text{ ,}

where \mathbf{J} = -\sigma\nabla u.

If I_0 were also unknown, you would need yet another scalar equation to close the system. (Otherwise, there would be infinitely-many solutions, since you could use the above formulation to find u and u_0 for every possible choice of I_0.)

1 Like
9

Hi,

Thanks for the explanation. I tried implementing the floating potential BC on the multi-physical piezo model.
The governing equations of the piezo model are shown in the image below:

PZ
PZ1105×734 79.8 KB

Currently am currently trying to conduct the stationary analysis of cube structure with piezoelectric material. I have the following boundary conditions:

  1. fixed displacement - bottom boundary
  2. Zero Voltage/Grounded - bottom boundary
    and I also want to implement the third BC, which is
  3. Floating potential or the surface should have equal potential distribution on the top boundary as defined in the figure above.
    I have currently implemented the piezoelectric effect with the first and second BCS. Am struggling to implement the floating potential impact on the top boundary marked. I have attached the test code below. I have tried to implement the floating potential effect as you suggested for the previous case. But I notice that the results are not matching with the test case. Looking forward to your suggestions.
from dolfin import *
import matplotlib.pyplot as plt
import numpy as np

#units: m, kg, s, V, K

################################## Mesh part ##################################

mesh = UnitCubeMesh(1,1,1)
Ue = VectorElement("CG", mesh.ufl_cell(), 1) # displacement vector element
Ve = FiniteElement("CG", mesh.ufl_cell(), 1) # voltage finite element
V0e = FiniteElement("Real", mesh.ufl_cell(),0) 
W = FunctionSpace(mesh, MixedElement([Ue, Ve, V0e]))
print(f"Number of dofs global: {W.dim()}")

# Parameters of Anisotropic Piezomaterial (PIC-255)
c_11 = (1.23e11)
c_12 = (7.67e10)
c_13 = (7.023e10)
c_33 = (9.711e10)
c_44 = (2.226e10)

# Relative permitivity matrix parameters
eps_0 = (8.854e-12) # relative permitivity of free space
eps_11 = (1649 * eps_0)
eps_33 = (1750 * eps_0)

# Coupling matrix parameters
e_13 = (-7.15)
e_33 = (13.75)
e_15 = (11.91)

# Order of matrices similar to ANSYS notation: x, y, z, xy, yz,xz
# IEEE notation: x, y, z, yz, xz, xy
# Elasticity tensor c^E
#IEEE Notation
C = as_tensor([[c_11, c_12, c_13, 0., 0., 0.], 
               [c_12, c_11, c_13, 0., 0., 0.],
               [c_13, c_13, c_33, 0., 0., 0.],
               [0., 0., 0., c_44, 0., 0],
               [0., 0., 0., 0., c_44, 0.],
               [0., 0., 0., 0., 0., (c_11-c_12)/2]])

# piezoelectric coupling tensor e
et_piezo = as_tensor([[0.0, 0.0, 0.0, 0.0, e_15, 0.0],
                     [0.0, 0.0, 0.0, e_15, 0.0, 0.0],
                     [e_13, e_13, e_33, 0.0, 0.0, 0.0]])

# transpose form of piezo tensor
e_piezo = as_tensor([[0.0, 0.0, e_13],
                    [0.0, 0.0, e_13],
                    [0.0, 0.0, e_33],
                    [0.0, e_15, 0.0],
                    [e_15, 0.0, 0.0],
                    [0.0, 0.0, 0.0]])

# Permittivitats tensor  epsilon^S
eps_s = as_tensor([[eps_11, 0., 0.],
                    [0., eps_11, 0.],
                    [0., 0., eps_33]])

# rewrite the tensor into a vector using Voigt notation
def strain3voigt(ten):
    # FEniCS does not know anything about Voigt notation, 
    # so one need to access the components directly as eps[0, 0] etc.   
    return as_vector([ten[0,0],ten[1,1],ten[2,2],2*ten[1,2],2*ten[0,2],2*ten[0,1]])

# rewrite a vector into a tensor using Voigt notation
def voigt3stress(vec):
    return as_tensor([[vec[0], vec[5], vec[4]], 
                      [vec[5], vec[1], vec[3]],
                      [vec[4], vec[3], vec[2]]])

# first define the function for the strain tensor in fenics notation
def B(u):
    return 0.5*(nabla_grad(u) + nabla_grad(u).T)

#define the electric field vector                 
def E_field_vector(Ei): 
    return as_vector([Ei[0], Ei[1], Ei[2]])
                  
#define the electric field in relation to potential
def E_field(v):
    return -grad(v)

#definition of the different components of my Stress tensor equation (sigma=C^E:Bu-e_piezo.E)

def sigma(u, v):
   return voigt3stress(dot(C, strain3voigt(B(u))) - e_piezo * E_field_vector(E_field(v)))

# the electric displacements vector
def disp_D(Di):
    return as_vector([Di[0], Di[1], Di[2]])

# definition of the different components of my electric displacements equation (D=e^s:Bu+eps^S.E)
# elect_disp_D for coupled linear elasicity and electrostatics
def D_field(u, v):
    D =   et_piezo * strain3voigt(B(u)) + eps_s * E_field(v)
    return disp_D(D) 

#Defining the "Trial" functions
#w = Function(W)
soln = Function(W)
(u, v, v0) = split(soln)

#Defining the test  functions
(wu, wv, wv0) = TestFunctions(W)

# Mark boundary subdomians
# Create classes for defining parts of the boundaries and the interior of the domain
class Bottom(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[2], 1)

class Top(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[2], 0.0)

# Initialize sub-domain instances
bottom = Bottom()
top = Top()

# Mark subdomains 
boundaries = MeshFunction("size_t", mesh, mesh.topology().dim()-1)
boundaries.set_all(0)
bottom.mark(boundaries, 2) # Displacement and Voltage BC
top.mark(boundaries, 3)  # Floating potential
ds = Measure('ds',subdomain_data=boundaries)

# Dirichlet boundary condition for fixed support and ground voltage
bcu = DirichletBC(W.sub(0), Constant((0,0,0)), bottom)
bcv1 = DirichletBC(W.sub(1), Constant(0), bottom)
bcs = [bcu, bcv1]

# Stationary analysis

n = FacetNormal(mesh)
Q0 = 1e-6 #C/m^2
meas_Gamma1 = Constant(1) # area of the edge or surface of current flow
D_avg = Constant(-Q0/meas_Gamma1)

Kuuv = inner(sigma(u,v), B(wu))*dx
Kvvu = inner(D_field(u,v),grad(wv))*dx  - (dot(D_field(u,v),n)*wv + dot(D_field(wu,wv),n)*(v-v0))*ds(3) + (D_avg - dot (D_field(u,v),n))*wv0*ds(3)

F = Kuuv + Kvvu 

# Set up the PDE
solve(F == 0, soln, bcs)
U, V, V0 = soln.split(deepcopy=True)
print(max(abs(U.vector().get_local())))
print(f"Voltage: {V.vector().get_local()}")

Results obtained from Fenics:
Displacement: 1.960682369212692e-08 m
Voltage: [ 0. 49.66534269 0. 0. 0. 46.62686436
38.87742427 46.62686436]

Correct results from commercial software:
Displacement: 2.0059095750665296E-8 m
Voltage: [ 0. 48.28872311077018 0. 0. 0. 48.28872311074099
48.28872311077021 48.28872311077022]

1 Like
10

Hi,

As explained in the previous reply, I tried implementing the floating potential condition in a multi-physical piezoelectric model. But the results don’t match. Can you check the implementation part? which will be really helpful.

11

Dear Sir,

The above implementation works fine within the FEnicS, but what in case we want to implement the same in DolfinX, as there is no option to use real elements space?

Regards,
GS

12

Dear @Gunasheela

I saw your code results are not the same as with commercial software. And your last reply said

The above implementation works fine within the FEnicS

May I ask where did you modified in the previous code?

Best regards