Electrical conductivity

Next: Elastic moduli Up: Finite element theory Previous: General aspects and

Electrical conductivity

Variables in programs

u_r = voltage at the r'th node in a pixel
σ_pq = conductivity tensor (differs in general from element to element)
= (E_x, E_y, E_z) = external electric field applied to image
e = (e_x, e_y, e_z) = local field at a point (x,y,z) inside a pixel
D_rs = stiffness matrix in a pixel
N_r(x,y,z) = shape array for cubic pixel

The purpose of this kind of finite element method is to express the total energy of the system in terms of the voltages only at the nodes, thus effectively discretizing the continuum. There are many ways in which to do this. We choose here to use a linear interpolation method, in what is often called a tri-linear scheme. Essentially we write out the energy of one pixel, which involves an integral over that pixel, and then combine the pixels of the whole system into a global energy functional. The fields inside the pixel are expressed as functions of (x,y,z) via a linear interpolation scheme in terms of the nodal voltages. The integral is carried out, and we end up with an expression for the energy of the pixel that is quadratic in the nodal voltages. The global energy is then also a quadratic functional of the nodal voltages. This expression is minimized with respect to the nodal voltages, using a conjugate gradient scheme. Using the set of nodal voltages found, the interpolation scheme can then be re-used to find the average current, total energy, etc. that can be used to define the effective conductivity and other quantities of interest (see Section 5.1).

Focus on a single pixel for now. Define V(x,y,z) to be the voltage within the given pixel, where 0 < x,y,z <1. The origin is taken where the pixel label is 1 (see Table 1 and Fig. 1). V(x,y,z) is determined by linear interpolation of the nodal voltages, such that

V(x,y,z) = N_ru_r

(2)

and N_r = N_r(x,y,z) . The actual functional forms of N_r in 3-D are:

N₁ = (1 − x)(1 − y)(1 − z)

N₂ = x(1 − y)(1 − z)

N₃ = xy(1 − z)

N₄ = (1 − x) y(1 − z)

N₅ = (1 − x)(1 − y) z

N₆ = x(1 − y) z

N₇ = xyz

N₈ = (1 − x)yz

Now the local field in the pixel, , is given in terms of derivatives of this local voltage, via

In terms of the nodal voltages, this expression for e_p (the p'th component of the local field

) then becomes

where − ∂N_r / ∂ x_p ≡ n_pr is a 3 x 8 matrix linking the 8−vector of nodal voltages to the 3−vector of local fields, and u_r has no explicit x dependence. Table 2 gives the components of this matrix.

Table 2: Formulas for components of [∂N_r / ∂x_i] matrix

Now, the total energy dissipated in the pixel is

In terms of the nodal voltages, this equation becomes

The indicated integrations can be easily carried out exactly using Simpson's rule, since they are at most quadratic in x, y, or z, and Simpson's rule is exact for quadratic expressions. These expressions could also be evaluated analytically, and then put into the program. When these programs were developed, it was decided to use Simpson's rule for simplicity, and to avoid doing all the derivations necessary to get the "stiffness" matrices in analytical form. The subroutine FEMAT in the finite element programs carries out this integration. The energy per pixel is then

where the matrix D_rs, an 8 x 8 matrix in 3-D, is defined by comparing eq. (7) to eq. (6). This matrix is known as a "stiffness" matrix, a term originating in finite element treatments of linear elasticity problems. The above term is added up for all pixels to give the global energy, each pixel using its own value of the conductivity tensor. This sum is then the total energy dissipated of the system under the applied field, which must be minimized with respect to the nodal voltages to find the solution to this discretized version of the original problem.

In order to be able to link up local energy with the global energy, it is important to see how the local numbering scheme for the pixel nodes links up with the global numbering scheme. The digital images are in the form pix(m), where pix is a 2-byte integer vector (defined INTEGER*2 in FORTRAN 77), and m is a one dimensional label. The digital image in 3-D is defined by (i,j,k), where i = 1, nx, j = 1, ny, and k = 1, nz, with the origin of the image being the (1,1,1) node. The one dimensional labelling scheme is defined by:

m = nx · ny · (k − 1) + nx · (j − 1) + i (8)

The label m refers to the m'th node, and the m'th pixel. The m'th pixel has the m'th node at its corner that is labelled "1" in the 1-8 local finite element numbering scheme. In 3-D , every node is part of 8 pixels. All the nodes in these pixels will need to be known with respect to the node, since they are connected to it in the global energy via the local stiffness matrices, and so their labels are encoded in the array ib. If the m'th node corresponds to (i,j,k), then the 27 nodes that are part of these 8 pixels are given by the m labels corresponding to adding or subtracting 0 or 1 from the i,j, k labels, so that ib = ib(m,27) has 27 entries for each value of m. Table 3 gives this neighbor labelling scheme. For example, ib(m,15) would be the m label of the node located at (−1,0,−1) with respect to the (i,j,k) position of the m'th node. Table 3 also gives the neighbor labelling scheme for 2-D, where Δk = 0, since there is no third dimension. Note that in the interior of an image, the m labels of the neighbors would always be simply known anyway, because of their fixed relationship in terms of i, j, and k, but with periodic boundary conditions (discussed below) a "neighbor" could be at the other side of the image. The relationship between the local (1-8) labelling, and the global (1-27) neighbor labelling used in array ib, is given in Table 4.

SRC="img86.gif"

Table 3: Neighbor labelling in 2-D and 3-D

Table 4: Relation between finite element labelling and neighbor labeling

Periodic boundary conditions mean that if a neighbor is outside the digital image, it is periodically continued on the opposite side of the system. For example, consider the node at i = 10, j = 12, k = 20, in an nx = ny = nz = 20 digital image. The m label of this node (pixel) would be 7830, via eq. (8). Neighbor number 18 would be located at i = 11, j = 13, and k = 21. However, k = 21 > nz, so the k label is changed to k − nz = 21 − 20 = 1. The m label of i = 11, j = 13, k = 1 is m = 251, so that ib(7830,18) = 251. The process is similar for the other faces (i = nx, j = ny), the edges (i = nx and j = ny, i = nx and k = nz, j = ny and k = nz) and the corner point (i = nx, j = ny, k = nz).

Before going into the solution technique, one must first consider the boundary conditions of the problem. Equation (7) is positive definite, and so has its minimum value when all the voltages are zero, an uninteresting result. We use periodic boundary conditions to apply an electric field, so that the energy is minimized with respect to an applied field, a more useful case.

Consider a pixel as before, with i = nx, j < ny, k < nz, and label m. When we go to evaluate the energy expression above, we find no voltages at the 2,3,6, and 7 nodes, because these would have i labels of (nx + 1), which is undefined. We use the voltages from across the system, numbers 1,4,5, and 8 from the pixel with the same j and k label but with i = 1, and label M = m − nx + 1. However, there is a jump of nx in length from where these voltages are defined. If the applied field (E_x,E_y, E_z) is present, then on average there is a drop in voltage of −E_x nx between these nodes. Therefore, the voltages to be used in the energy expression for the m'th pixel are:

u₁ = u₁(m),
u₂ = u₁(M) − E_x nx,
u₃ = u₄(M) − E_x nx,
u₄ = u₄(m),
u₅ = u₅(m),
u₆ = u₅(M) − E_x nx,
u₇ = u₈(M) − E_x nx,
and u₈ = u₈ (m).

We can write this, in general for a pixel at a boundary, as u_r = U_r + δ_r, where U_r is an 8−vector of the voltages that are given by the ib(m,n) labels, and δ_r is an 8−vector that corrects them to what they should be in the pixel in question. The neighbor array ib is designed to pick up the appropriate nodal voltages from across the image when evaluating the energy of pixels on faces, edges, and corner.

Inserting this choice of voltages into the expression for the energy involving the stiffness matrix, eq. (7), we get for the given pixel

which gives a term quadratic in the nodal voltages, a term linear in the nodal voltages, and a term constant with respect to the nodal voltages. In the particular case discussed, i = nx, j < ny, k < nz, the components of the 8−vector δ_r are: δ₁ = 0, δ₂ = −E_x nx, δ₃ = − E_x nx, δ₄ = 0, δ ₅ = 0, δ₆ = − E_x nx, δ₇ = −E_x nx, δ₈ = 0.

We can rewrite eq. (9) as

Adding these up over every pixel is done in the subroutine FEMAT, giving a global array b that gives the term in the energy that is linear in the voltages, and a global constant C. It is important to remember that the only contributions to b and C come from pixels having nodes at the unit cell boundaries and having a non-zero stiffness matrix. That makes it easy to make non-periodic boundary conditions by surrounding the system of interest by a layer of insulating material, one pixel thick, which effectively gets rid of the periodic boundary conditions and makes b and C equal to zero (see Sec. 6.2). Table 5 shows the values of the δ_r variable for the various faces and edges of the digital image.

Table 5: Values of Δ_r variable for system

Once the energy equation is set up, all that remains is to find the set of voltages that minimize the electrical energy dissipated. This is done by a conjugate gradient method, similar to that described in Ref. [5]. The subroutine DEMBX contains this routine, and is the same for both electric and elastic problems, aside from some minor labelling differences of indices, due to the fact that voltage is a scalar while elastic displacement is a 3-D vector. (see Section 4.2). No preconditioning of the matrix to be solved is done. It is not clear, since the full Hessian matrix is never stored, whether it would even be possible to carry out pre-conditioning in this case.

In the subroutine ENERGY, the gradient of the energy is computed, as this must become very small in order to solve the problem, which is exactly solved when the gradient is zero. Remember that the gradient is a vector containing all the partial derivatives of the energy with respect to all the nodal voltages. Using the above expression, the gradient of the energy will then be

where now A_mn and b_m should be thought of as global quantities, so that all the terms connected to u_m are included. The global matrix A_mn is of course built up from the individual D_rs matrices of the eight pixels that touch the node labelled m. The matrix A_mn is in principle large, but sparse. The way the finite element programs save memory is by computing A_mnu_n using only the small single-pixel stiffness matrices and the appropriate labelling scheme. This algorithm can be seen in subroutines ENERGY and DEMBX of the finite element programs. The correct D_rs term is used in place of the A_mn term that would be needed. In essence, the matrix A is re-built every time it is needed, without being stored.

Next: Elastic moduli Up: Finite element theory Previous: General aspects and