Electric current in anisotropic media

Question

I'm trying to model a bulk metal with anisotropic electrical conductivity such as graphite. I want to run a stationary current through the material and I am using the finite volume method to do this. My boundary conditions are fixed potentials at certain boundary positions of the simulated volume, just like electrical contacts. In a homogeneous isotropic medium, I simply use Poisson's equation (or Laplace to be more precise) to compute the potential and I use Ohm's law to derive the electric current from the potential.

Now on to an anisotropic medium: As far as I've understood, my conductivity becomes a tensor so that the current $j$ will be given as: $$j_i = \Sigma_{j} \sigma_{ij} E_j$$ This is where my confusion starts.

1. Wouldn't it make more sense to compute the current from the electric displacement field $\vec{D}$?
2. What happens to Poisson's equation? Do I also have to include an anisotropic permittivity in order to compute the electric field correctly? Or is this irrelevant for the electrical current? (It can't really be irrelevant since this would give me the same electric field for isotropic and anisotropic media)

UPDATE: Here is what I found out during further research and while running a few simulations on simple model systems:

Concerning 1.: Check the answer below.

Concerning 2.: The anisotropy needs to be included in Poisson's equation as well. Otherwise current is not conserved. In the case of graphite for example, the anisotropy coefficient for the conductivity can be applied to Poisson's equation as well. When the current is then derived from the gradient of the potential, I obtain the correct values.

Answer 1

Question 1:

One could define a linear relationship between $\vec{j}$ and $\vec{D}$ if one wanted to. This is because $\vec{D}$ will itself be related to $\vec{E}$ by a permittivity tensor: $$ D_i = \epsilon_{ij} E_j. $$
(I'll be using Einstein summation convention throughout this answer.) Assuming that the $\epsilon_{ij}$ is an invertible tensor, this implies that $$ E_i = (\epsilon^{-1})_{ij} D_j, $$ and so we can write $$ j_i = \sigma_{ij} (\epsilon^{-1})_{jk} D_k. $$ One could then define a new tensor $Q_{ik} = \sigma_{ij} (\epsilon^{-1})_{jk}$, which would give a linear relationship between $\vec{j}$ and $\vec{D}$. As to why this is not usually done, I suspect it's because the conductivity for isotropic media is usually written as $\vec{j} = \sigma \vec{E}$, and the definition you've written above is simply a generalization of this.

[As an aside, I'm not sure what you mean by "$\vec{D}$ is the physical quantity and not $\vec{E}$"; I'm much more used to a perspective that says that $\vec{E}$ is the more fundamental field of the two. See this answer for more discussion of this.]

Question 2:

(You already seem to have answered this, but I'll include it for completeness.)

According to conservation of charge, you should have $\vec{\nabla} \cdot \vec{j} = \partial_i j_i = 0$ in steady state. We also have $\vec{E} = -\vec{\nabla} \phi$ by definition. Combining the two, we see that the true condition on $\vec{\phi}$ in an anisotropic medium should be $$ \partial_i \left( \sigma_{ij} \partial_j \phi \right) = 0. $$ If the conductivity tensor is constant over the bulk, then we can pull it out of the derivative, obtaining $$ \sigma_{ij} \partial_i \partial_j \phi = 0. $$ Now, if it happens that the medium is isotropic, then we have $\sigma_{ij} = \sigma \delta_{ij}$, where $\delta_{ij}$ is the identity matrix. In this case, the above condition reduces to $$ \delta_{ij} \partial_i \partial_j \phi = \partial_i \partial_i \phi = \nabla^2 \phi = 0, $$ and $\phi$ obeys Laplace's equation. However, if the medium is anisotropic, one has to stick with the more complicated relationship involving $\sigma_{ij}$.

Note also that if one is trying to solve for $\phi$ in a medium with insulating boundary conditions, the form of these boundary conditions $\phi$ changes in the presence of an anisotropic medium. Specifically, the boundary condition is that $\hat{n} \cdot \vec{j} = 0$, where $\hat{n}$ is a normal vector to the surface; this guarantees that current is not flowing out of the medium. In terms of $\phi$, this becomes $$ n_i j_i = n_i \sigma_{ij} E_j = - n_i \sigma_{ij} \partial_j \phi = 0. $$ This means that $\partial_j \phi$ must be orthogonal to the vector $\sigma_{ij} n_j$, which is not necessarily parallel to $\hat{n}$ itself. It will be if the vector $\hat{n}$ is one of the principal axes of $\sigma_{ij}$; but one could easily imagine a situation where this is not the case.