### Correct use of variables to calculate Maxwell stress-tensor

172
views
1
3 months ago by
Hello,
I have prepared a simplified 3D model for computing the electrostatic field in a air-cube (basically the poisson example, please see below). This works really nicely!
Now I need to compute the force on one of the boundaries using the Maxwell-stress-tensor which can be calculated as
$\sigma_{\left\{ij\right\}}=\varepsilon_0E_iE_j-\frac{1}{2}\left(\varepsilon_0E^2\right)\delta_{ij}$σ{ij}=ε0EiEj12 (ε0E2)δij (1)
where E is grad(u) as shown in my example.
In my example I used project() and split() on E to have access to the x-, y- and z- components of E to compute the x-component of the force on the face using
Fx = epsilon_0 * assemble( ((E_x*E_x) * n_x + (E_x*E_y) * n_y + (E_x*E_z) * n_z - 1/2.* dot(E_p, E_p) * n_x ) *ds(2))
where nx, ny and nz are the components of the normal vector of the boundary.

Now the question is: how do I correctly use the E and n without projecting and splitting them?

I tried the following according to the above 'general' equation (1) but did not succeed
# 3 x 3 identiy matrix
I = Identity(3)
# Kronecker delta
delta_ij = I[i,j]

epsilon_0 = 8.854e-12

sigma = epsilon_0*( dot(E, E) ) + ( -1/2*epsilon_0* inner(E, E)* delta_ij )

I get the following error message:
Can't add expressions with different free indices.
Traceback (most recent call last):
File "poisson_3D_capacitor_minimal-v1p3.py", line 173, in <module>
sigma = epsilon_0*( dot(E, E) ) + ( -1/2*epsilon_0 * inner(E, E) * delta_ij)
File "/usr/lib/python2.7/dist-packages/ufl/exproperators.py", line 212, in _add
return Sum(self, o)
File "/usr/lib/python2.7/dist-packages/ufl/algebra.py", line 55, in __new__
error("Can't add expressions with different free indices.")
File "/usr/lib/python2.7/dist-packages/ufl/log.py", line 172, in error
raise self._exception_type(self._format_raw(*message))
ufl.log.UFLException: Can't add expressions with different free indices.​

HoWil

Simplified model (exported from jupyter-notebook)
from __future__ import print_function
# In[1]:

#get_ipython().magic(u'matplotlib notebook')
import matplotlib.pyplot as plt

# In[2]:

from dolfin import *

# # Setting up and solving the model

# In[3]:

length = 10.0  # x-dir in m, distance between plates
depth = 10.0  # y-dir in m, depth of capacitor plate
height = 10.0  # z-dir in m, heigth of capacitor plate

V1 = -3. # V, potential on first plate
V2 = 10. # V, potential on second plate

# Create classes for defining parts of the boundaries and the interior
# of the domain
class Left(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 0.0)

class Right(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], length)

# Initialize sub-domain instances
left = Left()
right = Right()

# In[4]:

# Define mesh in a uni-cube
mesh = BoxMesh(Point(0.0, 0.0, 0.0), Point(length, depth, height), 12, 12, 12)

# Initialize mesh function for interior domains
domains = CellFunction("size_t", mesh)
domains.set_all(0)

# Initialize mesh function for boundary domains
boundaries = FacetFunction("size_t", mesh)
boundaries.set_all(0)
left.mark(boundaries, 1)
right.mark(boundaries, 2)

# In[7]:

# Define input data
a0 = Constant(1.0)  # epsilon_r
g_L = 0
g_R = 0
f = 0

# Define function space and basis functions
V = FunctionSpace(mesh, "CG", 1)
u = TrialFunction(V)
v = TestFunction(V)

# Define Dirichlet boundary conditions at left and right boundaries
bcs = [DirichletBC(V, V1, boundaries, 1),
DirichletBC(V, V2, boundaries, 2)]

# In[9]:

# Define new measures associated with the interior domains and
# exterior boundaries
dx = Measure('dx', domain=mesh, subdomain_data=domains)
ds = Measure('ds', domain=mesh, subdomain_data=boundaries)

# In[10]:

# Define variational form
- g_L*v*ds(1) - g_R*v*ds(2)
- f*v*dx(0) )

# Separate left and right hand sides of equation
a, L = lhs(F), rhs(F)

# Solve problem
u = Function(V)
solve(a == L, u, bcs)

# # Plotting the solution

# In[23]:

plt.figure()
p = plot(u, title="u")
plt.colorbar(p)

plt.figure()

# # Comupting the force acting on a plate

# In[24]:

degree = V.ufl_element().degree()
W = VectorFunctionSpace(mesh, 'P', degree)

# In[25]:

E_p = project(E, W)

# ### Analytics

# The general form of the maxwell stress tensor is $$\sigma_{i j} = \epsilon_0 E_i E_j + \frac{1} # {{\mu _0 }}B_i B_j - \frac{1}{2}\bigl( {\epsilon_0 E^2 + \tfrac{1} # {{\mu _0 }}B^2 } \bigr)\delta _{ij}$$
#

# Without B it simplifies to $$\sigma_{i j} = \epsilon_0 E_i E_j - \frac{1}{2} \bigl( {\epsilon_0 E^2 } \bigr) \delta _{ij}$$

# In[26]:

E_x, E_y, E_z = E_p.split(deepcopy=True)  # extract components

n = FacetNormal(mesh)
n_x = variable(n)[0]
n_y = variable(n)[1]
n_z = variable(n)[2]

epsilon_0 = 8.854e-12

# ## Force acting on right boundary in x-direction

# In[27]:

Fx = epsilon_0 * assemble( ((E_x*E_x) * n_x + (E_x*E_y) * n_y + (E_x*E_z) * n_z - 1/2.* dot(E_p, E_p) * n_x ) *ds(2))
Fx # in Newton

# In[28]:

A = depth *height  # m2, area of capacitor plates
E_ = (V2-V1)/length  # V/m, electostatic field in x-dir

# In[29]:

F =  epsilon_0/2. * E_**2 *A  # N, resulting force
print("F_x = "+str(F)+" N")
​
Community: FEniCS Project

2
3 months ago by
If you replace the first dot(E,E) by outer(E,E)?
Hello Martin,
Can't add expressions with different shapes.

UFLExceptionTraceback (most recent call last)
<ipython-input-17-d1af42398493> in <module>()
----> 1 sigma = epsilon_0*( outer(E, E) ) + ( -1/2*epsilon_0* inner(E, E)* delta_ij )

210     if not isinstance(o, _valid_types):
211         return NotImplemented
--> 212     return Sum(self, o)
213
214

/usr/lib/python2.7/dist-packages/ufl/algebra.pyc in __new__(cls, a, b)
51         fid = a.ufl_index_dimensions
52         if b.ufl_shape != sh:
---> 53             error("Can't add expressions with different shapes.")
54         if b.ufl_free_indices != fi:
55             error("Can't add expressions with different free indices.")

/usr/lib/python2.7/dist-packages/ufl/log.pyc in error(self, *message)
170         "Write error message and raise an exception."
171         self._log.error(*message)
--> 172         raise self._exception_type(self._format_raw(*message))
173
174     def begin(self, *message):

UFLException: Can't add expressions with different shapes.

​
written 3 months ago by HoWil
If you use I directly instead of delta_ij?
written 3 months ago by Martin Genet

Dear Martin,
I do not really understand why should I use the outer product here.
Since the scalar product is  $\vec{a}\text{ }\cdot\vec{b}=a_1b_1+a_2b_2+a_3b_3$a ·b=a1b1+a2b2+a3b3   and I need (E_x*E_x) * n_x + (E_x*E_y) * n_y + (E_x*E_z)

Regaridng "If you use I directly instead of delta_ij?" This helps solving the equation with the use of outer product but I am not sure if the solution is correct!
How can I check the x-component of sigma on d(2)?

BR,
HoWil

written 3 months ago by HoWil
outer(a,b) = c with c_ij = a_i b_j. sigma = epsilon_0 * outer(E,E) - (epsilon_0/2) * inner(E,E) * I is the dolfin version of Equation (1).
written 3 months ago by Martin Genet
Thank you Martin,
I will test it again.
Do you have any advice how I can extract the x-component of Sigma. This would help me to compare it with Fx in my original post.
BR Howil

EDIT:
Finally found it myself... just insert the following code in my example directly below solve(a == L, u, bcs)...
E = grad(u)

epsilon_0 = 8.854e-12

I = Identity(v.geometric_dimension())

T = TensorFunctionSpace(mesh,'DG',2)

sigmaTensor = project(epsilon_0 * outer(E,E) - (epsilon_0/2) * inner(E,E) * I, T)

res_Fx = assemble(sigmaTensor[0,0]*ds(2))
res_Fx​

Thanks again!

written 3 months ago by HoWil