Mutual inductance of two disk coils (MATLAB)

[tommydyhr]

Hi,

I'm been spending way too much time this weekend, wrestling with the equations provided in the paper [The mutual inductance of two thin coaxial disk coils in air, Babic S., DOI: 10.1109/TMAG.2004.824810], but I just can't seem to verify the results obtained in the examples in the journal. Unfortunately, due to copyright reasons, I can't paste the equations here, but for those of you who have access to scientific journals (most college students, I suppose) can access it from the DOI.

To solve equation 3 in the paper, I've made the following script:

Code: [Select]

clear
clc

R1 = 0.1;
R2 = 0.2;
R3 = 0.3;
R4 = 0.4;
dist = 0.1;
N1 = 100;
N2 = 100;
mu_0 = 4.*pi.*10.^(-7);

p = [R1 R1 R2 R2];
l = [R3 R4 R4 R3];

sum_result = 0;
for n = 1:4
    J1F = @(beta) ((beta>0)&(beta<(pi/2))).*(1./sin(beta)).*(atan((dist.^2.*cos(beta)+p(n).*l(n).*(sin(beta)).^2)./(dist.*sin(beta).* (sqrt(l(n).^2+p(n).^2+dist.^2-2.*l(n).*p(n).*cos(beta))) ))-atan((dist.^2.*cos(beta)-p(n).*l(n).*(sin(beta)).^2)./(dist.*sin(beta).*(sqrt(l(n).^2+p(n).^2+dist.^2+2.*l(n).*p(n).*cos(beta)))))) + (beta==0)*((sqrt((l(n)+p(n))^2+dist^2)-(sqrt((l(n)-p(n)).^2+dist.^2)))./dist) + (beta==pi/2).*2.*atan((l(n).*p(n))./(dist.*sqrt(l(n).^2+p(n).^2+dist.^2)));
    J2F = @(beta) asinh((l(n)+p(n).*cos(2.*beta))./sqrt(p(n).^2.*(sin(2.*beta)).^2+dist.^2));
    J3F = @(beta) asinh((p(n)+l(n).*cos(2.*beta))./sqrt(l(n).^2.*(sin(2.*beta)).^2+dist.^2));
    J1 = integral(J1F,0,pi/2);
    J2 = integral(J2F,0,pi/2);
    J3 = integral(J3F,0,pi/2);
    k = sqrt((4*p(n)*l(n))/((l(n)+p(n))^2+dist^2));
    m = 2*p(n)/(sqrt(p(n)^2+dist^2)+p(n));
    pn = 2*l(n)/(sqrt(l(n)^2+dist^2)+l(n));
    theta1 = asin(abs(dist)/(sqrt(p(n)^2+dist^2)+p(n)));
    theta2 = asin(sqrt((1-m^2)/(1-k^2)));
    theta3 = asin(abs(dist)/(sqrt(l(n)^2+dist^2)+l(n)));
    theta4 = asin(sqrt((1-pn^2)/(1-k^2)));
    T = ellipticE(k)*(-2*l(n)*p(n)*sqrt(l(n)*p(n))/k) + k*sqrt(l(n)*p(n))*(dist^2-l(n)^2-p(n)^2+p(n)*sqrt(p(n)^2+dist^2)+l(n)*sqrt(l(n)^2+dist^2))*ellipticK(k)-pi*abs(dist)*0.5*p(n)*sqrt(p(n)^2+dist^2)*(1-Heuman_Lambda(theta1,k)-sign(sqrt(p(n)^2+dist^2)-l(n))*(1-Heuman_Lambda(theta2,k)))-pi*abs(dist)*0.5*l(n)*sqrt(l(n)^2+dist^2)*(1-Heuman_Lambda(theta3,k)-sign(sqrt(l(n)^2+dist^2)-p(n))*(1-Heuman_Lambda(theta4,k)))-0.5*dist^3*J1+p(n)^3*J2+l(n)^3*J3;
    sum_result = sum_result + ((-1)^(n-1))*T;
end

M = ((mu_0*N1*N2)/(3*(R2-R1)*(R4-R3)))*sum_resultWhere I've defined the Heuman Lambda function in a separate file:

Code: [Select]

function [ HL ] = Heuman_Lambda( t,m )

mdash = (1-m);

[K,E] = ellipke(m);
[K2,E2] = ellipke(mdash);
incF = ellipticF(t,mdash);
incE = ellipticE(t,mdash);

Z = incE - (E2*incF)/K2; % Jacobi zeta function

HL = incF / K2 + (2/pi) * K * Z;
Whichever of the examples I try, I always seem to get a complex result, as opposed to the 1.226 mH from the paper.

If anyone can spot any glaringly obvious mistakes from the code, I'd be forever in your debt!

I suppose it's worth noting that with the values provided (ie. R1 = 10cm, R2 = 20cm, R3 = 30cm, R4 = 40cm, dist = 10cm), theta1 and theta2 returns complex values due to arcsin(x) not being defined in R for x > 1; eg.



I also tried solving the examples with the simple summation equation (4), which was realized with the following:

Code: [Select]

clear
clc
R1 = 0.0762; R2 = 0.1594; R3 = 0.0762; R4 = 0.1594;
RI = (R2+R1)/2; RII = (R3+R4)/2;
hI = R2-R1; hII = R4-R3;
dist = 0.0468;
N1 = 516; N2 = 516;
mu_0 = 4*pi*10^(-7);
N = 5; n = 5;
sum = 0;
for h = -N:N
  for l = -n:n
      Rh = RI + (hI*h)/(2*N+1);
      Rl = RII + (hII*l)/(2*n+1);
      k = sqrt(4*Rh*Rl/((Rh+Rl)^2+dist^2));       
        sum = sum + (2*mu_0*sqrt(Rh*Rl))/k * ((1-0.5*k^2) * ellipticK(k) - ellipticE(k));
    end
end
M = ((N1*N2)/((2*N+1)*(2*n+1)))*sum This should yield an M (mutual inductance) of 0.0366 H , but no matter how I arrange the equation, I always get a result of 0.0528 H.

Any help or suggestions as to some other way to solve this problem would be greatly appreciated!

[German_EE]

Damn! My migraine is back

[gregariz]

4 coils?

I don't have time to dig in but I would attempt to validate the 2 coil case with maxwell's circular filament equations first. (his 2nd Volume IIRC)

[tommydyhr]

Sorry about that, German_EE. I would offer to take some of that pain, but I'm afraid my head doesn't have capacity for any more pain at this point! :<

gregariz, I'm just using two coils. Initially I'm just examining the case of two coaxial disk coils separated by some distance. If I eventually do succeed, I'll move on to study the noncoaxial case. My second attempt (the code at the bottom of my original post) is actually realizing Maxwell's filament method, but alas, that didn't work either.

[T3sl4co1l]

I don't think equations can be copywritten; anyway, small excerpts are generally acceptable under fair use.

And, that's MATLAB / Octave code, I think?  Looks familiar, though I don't recognize the @(...) expression.

Tim

[Bored@Work]

And, that's MATLAB / Octave code, I think?  Looks familiar, though I don't recognize the @(...) expression.

It is an anonymous function. The fact that the OP defines three anonymous functions, dynamically in a loop (making it 12 in total), tells me that code is too clever for its own good.

Further, the use of constructs like

Code: [Select]

(beta==0) or

Code: [Select]

(beta==pi/2) is literally begging for trouble with floating point numbers. There is also constant changing between scalar and vector operations. I have a hunch that a few of the operators are of the wrong type.

And with this collection of trigonometric functions I would also suspect more numeric instabilities and rounding errors all over the place.

So, how to debug that? One expression after the other, inside out. It would be tremendous helpful to have the intermediate results from one working example to compare the intermediate results from that code with. How to get one working example? Manual calculation ... Good luck with that.

[KJDS]

I've sacked people for coding like that.

Having both R1 and RI as variables?

[IanB]

Unfortunately, due to copyright reasons, I can't paste the equations here

Why do you think that? Equations can't be copyrighted, they are just mathematical expressions. It will be very hard for you to get any help unless you post the equations you are trying to solve along with the numerical parameters of the problem.

[tommydyhr]

I should've prefaced the OP by saying I'm by no means experienced with Matlab at all. I've worked with the symbolic toolbox (MyPad) for a while, but that's about the extent of it.

Thanks for all the responses; much appreciated! I'm not that savvy with copyright-laws, so I'd rather be safe than sorry. The equations are derived in the paper (and that's what the actual paper is about, really), so I wasn't sure as to if I could post them verbatim.

KJDS, I know that R1/RI and R2/RII could be rather confusing, but I named them like that intentionally - not to be confusing, but to directly match the equations and variables in the paper.

Bored@Work, what would happen if I were to, say, perform a vector operation on a scalar (ie. a.^b as opposed to a^b, a and b being scalars)?

For the example I initially tried to examine, the following were given: R1 = 10cm, R2 = 20cm, R3 = 30cm, R4 = 40cm, zQ = 10cm, N1 = 100, N2 = 100. This should yield a mutual inductance of M = 1.226 mH.

The complete equation is defined as:

[tommydyhr]

I suspect there may be something about elliptic integrals I don't quite understand.

Given the following equation:



where K(f) and E(f) are the complete elliptic integrals of the first and second kind, and mu₀=4*pi*10^-7 H/m

We should be able to calculate the mutual inductance of two aligned circular filaments.

We know that mutual inductance decreases with an increase of the distance between the two coils, however the following code:

Code: [Select]

x = 0:0.01:1;
k = sqrt( (4.*dT1.*dR1)./((dT1+dR1).^2+x.^2));
y = (1/2) .* mu_0 .* sqrt(dT1*dR1) .* ( (2./k-k).*ellipticK(k) - (2./k).*ellipticE(k));
plot(x,y)
yields the following plot:

[The Electrician]

I see a more recent paper by Babic which appears to be freely available online:

http://www.wseas.us/e-library/conferences/2006istanbul/papers/520-115.pdf

If you could use this paper as your reference then those who would help you could see the full paper under discussion.

This paper considers several coil configurations, including the thin disc coils configuration.

You haven't shown the physical disposition of your coils.  It would be nice if you provided a diagram such as in figure 3 of the referenced paper.

[tommydyhr]

I see a more recent paper by Babic which appears to be freely available online:

http://www.wseas.us/e-library/conferences/2006istanbul/papers/520-115.pdf

If you could use this paper as your reference then those who would help you could see the full paper under discussion.

This paper considers several coil configurations, including the thin disc coils configuration.

You haven't shown the physical disposition of your coils.  It would be nice if you provided a diagram such as in figure 3 of the referenced paper.

Great, thanks! My most recent post with the plot is a plot of equation (1) from that paper, which behaves quite opposite of what I'd expect, with M increasing with increasing distance.

[The Electrician]

In reply #9 above, you didn't say what values you used for dT and dR.

I've noticed that modern mathematical software sometimes expects the modulus of elliptic integral functions to be the square, or square root, of what you might expect.  Compare what you get for EllipticE(.5) and EllipticK(.5) with what Mathematica gives me:



Also, here's what I get when I plot Maxwell's elliptic formula for radii 10 and 10 versus separation:



You should download this classic paper for your reference:

http://www.g3ynh.info/zdocs/refs/NBS/Sci169noerr.pdf

[The Electrician]

You might also try plotting just the elliptic functions to make sure that is working OK:

[tommydyhr]

In reply #9 above, you didn't say what values you used for dT and dR.

I've noticed that modern mathematical software sometimes expects the modulus of elliptic integral functions to be the square, or square root, of what you might expect.  Compare what you get for EllipticE(.5) and EllipticK(.5) with what Mathematica gives me:

...

Also, here's what I get when I plot Maxwell's elliptic formula for radii 10 and 10 versus separation:

...

You should download this classic paper for your reference:

http://www.g3ynh.info/zdocs/refs/NBS/Sci169noerr.pdf

Thank you. I get the same results (and indeed the same plot) as you for radii 10 and 10. My example in #9 was for radii 17 mm and 3.95 mm.

[IanB]

Thank you. I get the same results (and indeed the same plot) as you for radii 10 and 10. My example in #9 was for radii 17 mm and 3.95 mm.

In that example from post #9, does the mutual inductance continue to increase with increasing separation beyond 1, or does it start to decrease?

[tommydyhr]

Thank you. I get the same results (and indeed the same plot) as you for radii 10 and 10. My example in #9 was for radii 17 mm and 3.95 mm.

In that example from post #9, does the mutual inductance continue to increase with increasing separation beyond 1, or does it start to decrease?

It converges at y = 1.6E-8:

[IanB]

It converges at y = 1.6E-8

That result doesn't follow from the equation, though. Consider:



When h increases, the arguments to each of the elliptic integrals approach zero, and both integrals become equal at zero:



Furthermore the multiplying term in front of the first integral approaches 1, meaning that the difference term inside the parentheses should approach zero. As a result the whole expression should tend to zero with increasing h.

[tommydyhr]

It converges at y = 1.6E-8

That result doesn't follow from the equation, though. Consider:

...

When h increases, the arguments to each of the elliptic integrals approach zero, and both integrals become equal at zero:

...

Furthermore the multiplying term in front of the first integral approaches 1, meaning that the difference term inside the parentheses should approach zero. As a result the whole expression should tend to zero with increasing h.

You're absolutely right. Further examination of all of the equations (including plotting every single part separately) showed me that in order for the equations to match the results from the papers, the argument for the elliptic functions shall be squared. Now all my previous work actually matches the papers' results as well.

This should feel like a success, really, but after this many hours... I can't believe it was just that.

Examining Matlab's definition of elliptic integrals does indeed confirm this fact, though:
Matlab's definition:


Wikipedia:


Thanks for all the help, guys, I can't begin to describe how much it's appreciated. I might even be able to sleep tonight!

[The Electrician]

You're absolutely right. Further examination of all of the equations (including plotting every single part separately) showed me that in order for the equations to match the results from the papers, the argument for the elliptic functions shall be squared. Now all my previous work actually matches the papers' results as well.

This should feel like a success, really, but after this many hours... I can't believe it was just that.

I found this out myself years ago after many hours of trying to figure out why my results didn't match those from some paper. 

That's why I warned you about it in reply #12, and suggested you check the values you were getting for EllipticE[.5) and EllipticK[.5].

I had a feeling that was your problem.   But now it's more like:   

[tommydyhr]

You're absolutely right. Further examination of all of the equations (including plotting every single part separately) showed me that in order for the equations to match the results from the papers, the argument for the elliptic functions shall be squared. Now all my previous work actually matches the papers' results as well.

This should feel like a success, really, but after this many hours... I can't believe it was just that.

I found this out myself years ago after many hours of trying to figure out why my results didn't match those from some paper. 

That's why I warned you about it in reply #12, and suggested you check the values you were getting for EllipticE[.5) and EllipticK[.5].

I had a feeling that was your problem.   But now it's more like:   

I actually did compare my results to the results from Mathematica, which is why I dismissed the idea in the first place.
Then again, let's look at Mathematica's definition of the first elliptic integral:


Ah well, you live and you learn! Thanks again!