# Capacitance in freefem++

I would like to simulate a capacitor in 2d with freefem++. This is the code I used:

int C0=0;
int C1=1;
int C2=2;

border C01(t=0,4.25){x=t; y=0; label=C0;}
border C02(t=0,2.25){x=4.25; y=t;label=C0;}
border C03(t=4.25,0){x=t; y=2.25; label=C0;}
border C04(t=2.25,0){x=0; y=t; label=C0;}

border C11(t=0.25,0){x=1.75+t; y=0.75; label=C1;}
border C12(t=0,0.75){x=1.75; y=0.75+t; label=C1; }
border C13(t=0,0.25){x=1.75+t; y=1.50;label=C1; }
border C14(t=0,0.75){x=2; y=1.5-t;  label=C1;}

border C21(t=0.25,0){x=2.25+t; y=0.75; label=C2;}
border C22(t=0,0.75){x=2.25; y=0.75+t; label=C2; }
border C23(t=0,0.25){x=2.25+t; y=1.50;label=C2; }
border C24(t=0,0.75){x=2.5; y=1.5-t;  label=C2;}

Th=buildmesh(C01(50)+C02(50)+C03(50)+C04(50)+C11(+20)+C12(+40)+C13(+20)+C14(+40)+C21(+20)+C22(+40)+C23(+20)+C24(+40));

fespace Vh(Th,P2);
Vh uh,vh;

real voltc1=1;
real voltc2=-1;

//1 e -1 volt dirichlet conditions
problem Poisson(uh,vh)=int2d(Th)(dx(uh)*dx(vh)+dy(uh)*dy(vh))+on(C0,uh=0)+on(C1,uh=voltc1)+on(C2,uh=voltc2);

Poisson;

plot(uh,wait=true,fill=true,value=1);

real energy=int1d(Th) ((dx(uh))^2+(dy(uh))^2);

cout<<"Energy (double)"<<endl<<"E =  "<<energy<<endl;

cout<<"Capacitance"<<endl<<"C =  "<<energy/(voltc1-voltc2)^2<<endl;


Runnng this code it gives me a value for the capacitance of 36.319 F. Is it this value reliable? Using the simplest formula $C=\epsilon \frac{A}{d}$ with $A$ area and $d$ the distance between the two plates the result is completely different. What is wrong in all of this?

• Are you sure about the way you compute the energy? You use int1d, but Th is a 2D domain. – Dr_Sam Oct 18 '13 at 7:58
• Yes this is already a mistake. Thanks a lot. Changing to int2d the new value for the capacitance is $C=5.156 F$. I think the right value should be around 3 because as $A$ I took a value of 0.75 while $d=0.25$ and $\epsilon=1.0$. – Rob Oct 18 '13 at 9:27
• For what it's worth, the computed capacitance will always be higher than your L/d formula, because the latter ignores any additional energy stored in the "fringing fields" outside the region between the plates. As the answer below suggests, your formula is really an approximation for large plates separated by small thickness, so you might have better agreement when your discrete model is shaped more like that. – rchilton1980 Oct 18 '13 at 12:52
• Also, a picture is worth a thousand words here - it's not immediately clear to me what your capacitor's dimensions are from your source code. – rchilton1980 Oct 18 '13 at 12:55
• Yes You were right...is the thickness the problem...I tried reducing it and now I have a value really close to 3. Thanks. – Rob Oct 18 '13 at 14:54

One of the problems (except for the 1d integration in the code) is that you use a formula that is not valid for your test case.

Indeed, the formula that you use is made for two parallel plates

• infinitly thin
• with a finite area

so, they should be like the two super thin slices of bread of a square burger.

What you represent with your 2D case are two parallel plates

• which are thick (0.25 thick)
• with an inifite area, because you neglect the third dimension (so it's like it would be infinite in that direction).

So the formula just does not apply.

To check your 2D implementation, you should better create a test for the coaxial cable (formula found here on wiki) for which it is much easier to restrict to a 2D test case.

• I think his formula (or a slight variant) is still applicable in 2D, to determine the capacitance per-unit-length, no? – rchilton1980 Oct 18 '13 at 12:46
• There might be a formula yes, but in my opinion, not the one he uses: the area in the formula is the area that one plate shows to the parallel one, not the area of a cut. – Dr_Sam Oct 18 '13 at 13:03
• I think there should be an analogous formula yes. For the coaxial case the 2d geometry would be rapresented by two concentric circles. The problem is that the value of $l$ would be infinity also in this case; so how can I apply the formula in the wiki for the capacitance? – Rob Oct 18 '13 at 13:15
• Yes, you need $l$, sorry for that, but I think that you should be able to extract the capacitance by considering that the solution is the same for all the cuts of the cable. So, finally, it's like computing the capacitance per unit length as @rchilton1980 proposed for the plates (but it is maybe easier to realize how to create the geometry). – Dr_Sam Oct 18 '13 at 13:20

If I use the new geometry for the coaxial case defined like this:

border C01(t=0,5){x=t; y=0; label=C0;}
border C02(t=0,5){x=5; y=t;label=C0;}
border C03(t=5,0){x=t; y=5; label=C0;}
border C04(t=5,0){x=0; y=t; label=C0;}
border C11(t=0,2*pi){x=2.5+sin(t);y=2.5+cos(t);label=C1;};
border C21(t=0,2*pi){x=2.5+1.1*sin(t);y=2.5+1.1*cos(t);label=C1;};
border C31(t=0,2*pi){x=2.5+0.3*sin(t);y=2.5+0.3*cos(t);label=C2;};


I obtained with the same code a value for the capacitance of $C=6.96877 F$. I still don't see the relation with the formula in the wiki...

• The coaxial cable is a little different - the outer conductor isolates the problem from the exterior. (You should only have the black region, the red region is not needed). The poission problem you should solve is over the black region, with a dirichlet value of 1 on the inner conductor and a dirichlet value of 0 on the outer (yielding a potential difference of 1 volt). Solve for the unknown potentials in the black region, take the gradient of that field to determine the E-field everywhere. Integrate over E to find the stored energy. With energy and voltage, you can determine C. – rchilton1980 Oct 18 '13 at 17:20
• The formula on wikipedia is for the total capacitance of a length l of coaxial cable. Your 2D code does not solve that - it finds capacitance per unit length. (So drop the factor of l from the wikipedia formula to make things match). Your capacitance per unit length should be 2*pi*epsilon/ln(R2/R1). And your code should compute this almost exactly, because the coaxial case does not have any parasitic fringing effects. – rchilton1980 Oct 18 '13 at 17:25
• This information should also be edited into your original post. – Geoff Oxberry Oct 18 '13 at 17:56
• Yes still following your instruction I obtained a perfect match! Thanks a lot. Can I ask You if You know some links where I can find more electromagnetic test cases to work with? I would also like to implement things in 3d and also with other kind of elements like Nedelec and/or Raviart Thomas. – Rob Oct 18 '13 at 19:49
• You can piece together the 3D sphere case immediately - it's just a 3d shell instead of a 2D shell. As far as next steps, it depends upon what you want to model. Nedelec elements are used mainly for magnetostatics (electric motors/machines/solenoids) and the vector helmoltz equation(waveguide discontinuities/antennas). RT is more used for surface integral equations (scattering/antennas). There's also the time domain/frequency domain schism when you start talking about non-static problems, which influences what you should study on the linear algebra side (iterative/gmres? direct/multifrontal?) – rchilton1980 Oct 18 '13 at 20:38

Ok for the plate capacitor I solved following the suggestion of @rchilton1980 using a new geometry like this one. The thickness is really small compared to the side of the plate so the approximation works and I obtained a capacitance of around 3.

• This information should be edited into your original post. – Geoff Oxberry Oct 18 '13 at 17:56