Re: [Libmesh-users] additional term in laplace equation

[Mark B. <mb...@au...>]

Hi Michael,

I am using libmesh to solve for the electric potential of a current source 
within a 3d cubic domain. The calculations are used for a geoelectrical 
forward solver. On the top side of the domain, I apply a Neuman boundary 
condition (du/dn=0) and on the other sides I use Dirichlet boundary 
conditions (u=0).
The second term 
div( (sigma-sigma0) grad Un)
is added to remove the singularity arising due to the point current source. 
When I discretize the second term with 

 Fe(i) += - (sigma-sigma0) * JxW[qp] * (gradVn[qp] * dphi[i][qp])
(of course, as John noticed, I forgot to mention the surface integral term, 
but in my case it is zero anyway cause sigma=sigma0 at the outer boundary. )

I get wrong results: The resulting potential U is the potential that does not 
contain the singular potential Un=1/(2*pi) *1/r anymore.
For the case of a homogenous cube with a high conducting cube centered inside, 
the introduction of the discretization above results in a strange dipole 
field within the high conducting cube, which physically doesnt make any 
sense. 
Thanks for the help so far,

Mark



Am Montag 10 Januar 2005 19:28 schrieb Michael Schindler:
> Hi Mark,
>
> On 10.01.05, Mark Blome wrote:
> > Hi everybody,
> >
> > I am trying to descretize a laplace equation which contains a second
> > term. In my opinion the problem should be very easy, but I become more
> > and more confused the more I think about it. The laplace equation reads:
> > div(sigma grad U) + div( (sigma-sigma0) grad Un) =0
> > where Un is a potential for which the equation is known ( Un=j0/(2 * PI *
> > sigma0)  * 1/r), hence the term has to go to the right hand side.
> >
> > For the term div(sigma grad U)=0 I simply get matrix terms:
> > 	    Ke(i,j) += - sigma * JxW[qp]*(dphi[i][qp]*dphi[j][qp]);
> > But what to add for div( (sigma-sigma0) grad Un)?
>
> The FEM method is to find a weak solution with a test-function phi
>
>   int div(sigma grad U) phi  =  int div( (sigma-sigma0) grad Un) phi
>
> Depending on sigma you can use the Galerkin-ansatz where phi is the
> same function you approximate U with. If sigma is constant, omitting
> the surface terms.
>
>   - sigma int grad(U)*grad(phi) = - (sigma-sigma0) int grad(Un)*grad(phi)
>
> This is what you tried first.
>
> > I tried:
> >         Fe(i) += - (sigma-sigma0) * JxW[qp] * (gradVn[qp] * dphi[i][qp])
> > ; where gradVn[qp] simply is grad(Un) analytically calculated at the
> > quadrature points. But it doesnt work like this.
>
> What actually does not work? What boundary conditions do you use?
>
> Michael.