Hello. I tried to calculate charge using volume integral of Laplacian of potential and got strange results.
In order to simplify the task, I used 2d case here (capacitor model). The results are similar to the 3d case. First, they are inaccurate. Second, they depend on the type of finite elements I use.
In the example below, capacitance is calculated in two ways: using energy and using charge. The first approach gives close results to the ideal capacitor (neglecting errors due to usage of the 2d model, etc.). The second approach gives several times smaller capacitance (when P2 is used). For other finite elements (P3, P4, etc.), the results are different.
Should I use some particular finite element type for this kind of problem? Or I should use some completely different approach to calculate the charge (if it’s possible)?
(Also, I don’t know why, but I got wrong results with vcoef=1, but that’s not important for me currently)
Output:
Capacitance (ideal) = 4.425e-11
Capacitance (energy) = 4.4103e-11
Capacitance (integr) = 8.6263e-12
Code:
int CA=3, CK=4, CE=5, CV=6;
real w2=1.0, h=0.01, d2=0.2;
border bottomA(t=-w2,w2) {x=t; y=d2; label=CA;};
border rightA(t=d2,d2+h) {x=w2; y=t; label=CA;};
border topA(t=w2,-w2) {x=t; y=d2+h; label=CA;};
border leftA(t=d2+h,d2) {x=-w2; y=t; label=CA;};
border bottomK(t=-w2,w2) {x=t; y=-d2-h; label=CK;};
border rightK(t=-d2-h,-d2) {x=w2; y=t; label=CK;};
border topK(t=w2,-w2) {x=t; y=-d2; label=CK;};
border leftK(t=-d2,-d2-h) {x=-w2; y=t; label=CK;};
real vcoef = 0.999;
border V1(t=-w2,w2) {x=t; y=vcoef*d2; label=CV;}
border V2(t=0,vcoef) {x=w2; y=vcoef*d2-2*d2*t; label=CV;}
border V3(t=w2,-w2) {x=t; y=-vcoef*d2; label=CV;}
border V4(t=0,vcoef) {x=-w2; y=-vcoef*d2+2*d2*t; label=CV;}
border enclosure(t=0,2*pi) {x=5*cos(t); y=5*sin(t); label=CE;}
mesh Th = buildmesh(enclosure(100)+
bottomA(-10)+topA(-10)+rightA(-2)+leftA(-2)+
bottomK(-10)+topK(-10)+rightK(-2)+leftK(-2)+
V1(-10)+V2(-20)+V3(-10)+V4(-20));
fespace Vh(Th,P2);
Vh u,v;
solve a(u,v) = int2d(Th) (dx(u)*dx(v)+dy(u)*dy(v)) + on(CA,u=1.0) + on(CK,u=0);
cout << "Capacitance (ideal) = " << (w2*2)*(8.85e-12)/(d2*2) << endl;
Vh envol = (x<w2)*(x>-w2)*(y<d2)*(y>-d2);
real energy = 0.5*(8.85e-12)* int2d(Th) (envol*(dx(u)^2+dy(u)^2));
cout << "Capacitance (energy) = " << 2.0*energy/(1.0^2) << endl;
Vh smvol = ( (y>(d2-h*10)) * (y<(d2+h*10)) * (x>(-w2*2)) * (x<(w2*2)) );
real integr = int2d(Th)( smvol*(dxx(u)+dyy(u)) );
cout << "Capacitance (integr) = " << integr*(8.85e-12) << endl;