Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

MATLAB Mutual inductance of two disk coils (MATLAB)

  1. Oct 12, 2014 #1

    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 (Text):

    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;

    M = ((mu_0*N1*N2)/(3*(R2-R1)*(R4-R3)))*sum_result
    Where I've defined the Heuman Lambda function in a separate file:
    Code (Text):
    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.
  2. jcsd
  3. Oct 17, 2014 #2
    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?
  4. Oct 18, 2014 #3
    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! :rolleyes:
Share this great discussion with others via Reddit, Google+, Twitter, or Facebook