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:
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_result
Where I've defined the Heuman Lambda function in a separate file:
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!
![Smiley :)](https://www.eevblog.com/forum/Smileys/default/xsmiley.gif.pagespeed.ic.R8GFI-pF6f.png)
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.
![](http://i.imgur.com/SoQzgRZ.png)
I also tried solving the examples with the simple summation equation (4), which was realized with the following:
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!