MATLAB 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:

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! :)

Edit: 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.

 

[Greg Bernhardt]

Thanks for the post! Sorry you aren't generating responses at the moment. Do you have any further information, come to any new conclusions or is it possible to reword the post?

 

[tommydyhr]

Thanks Greg. I found a completely alternative solution, which I can't really elaborate on here. I tried to delete this topic a couple of days ago, but I couldn't find the button for it!