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

172

views

1

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}$

where

In my example I used project() and split() on

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

I tried the following according to the above 'general' equation (1) but did not succeed

I get the following error message:

Thanks in advance,

HoWil

Simplified model (exported from jupyter-notebook)

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}}=`ε`_{0}`E`_{i}`E`_{j}−12 (`ε`_{0}`E`^{2})`δ`_{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.
```

Thanks in advance,

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
F = (inner(a0*grad(u), grad(v))*dx(0)
- 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]:
# Plot solution and gradient
plt.figure()
p = plot(u, title="u")
plt.colorbar(p)
plt.figure()
E = grad(u)
p_grad = plot(E, title="E-field - Projected grad(u)")
# # 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

### 1 Answer

2

If you replace the first

`dot(E,E)`

by `outer(E,E)`

?
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`=`a`_{1}`b`_{1}+`a`_{2}`b`_{2}+`a`_{3}`b`_{3} 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)?

Please double-check your answers with the code I posed above.

Thank you in advance,

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.

Thank you in advance,

BR Howil

EDIT:

Finally found it myself... just insert the following code in my example directly below

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.

Thank you in advance,

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

Please login to add an answer/comment or follow this question.

Thanks four your answer, but it gives a different error message: