sum of a sum of a sum of a sum

Hi, I'm working on implementing a set of equations to find how the mutual inductance changes for multilayer misaligned coils.

After implementing all equations I try to complete the last step:

but i seem to be getting a different result than the one shown in the paper (M = -1.42256038 µH). In my case the values of S,N,m,n have a massive impact of the final result as well, which doesn't seem to be the case in the paper.

Here is my code for this part:

M = @(a,b,c,xc,yc,zc,q,g,p,h,phi,K,E) N1*N2 .* M0(a,b,c,xc,yc,zc,q,g,p,h,phi,K,E)./(2*n+1)./(2*N+1)./(2*S+1)./(2*m+1);
M_temp=0;
for i = -n:1:n
    q = i;
    for ii = -m:1:m
        p = ii;
        for iii = -N:1:N
            h = iii;
            for iv = -S:1:S
                g = iv;
                
                k = sqrt(abs(k_squared(a,b,c,xc,yc,zc,g,p,h,phi)));
                [K,E] = ellipke(k);
                M_temp = M_temp + M(a,b,c,xc,yc,zc,q,g,p,h,phi, K, E);
            end
            
            M_temp = M_temp + M(a,b,c,xc,yc,zc,q,g,p,h,phi, K, E);
        end  
        M_temp = M_temp + M(a,b,c,xc,yc,zc,q,g,p,h,phi, K, E);
    end 
    M_temp = M_temp + M(a,b,c,xc,yc,zc,q,g,p,h,phi, K, E);
end  

[Marta Gorecka]

Here is the latest version of the entire code, as suggested by @David Goodmanson. The value still isn't the same sadly, but seems much closer to the corect result.

Any advice would be much appreciated!

clc; clear all;
%% self inductance of the coils
%https://ieeexplore.ieee.org/stamp/stamp.jsp?arnumber=1669896
%all dimensions in inches
a = 1.476;
b = 0.315;
c = 0.945;
N12 = 27; %total number of turns
mu0 = 4*pi*10^(-7); %permeability of free space in [H/m]
%Wheeler's formula for multi-layer coil
%in microhenries
L_Wheeler = 0.8*a.^2*N12^2/(6*a+9*b+10*c);
disp('L_Wheeler = '); disp(L_Wheeler);

%% Mutual inductance calculations based on https://www.researchgate.net/publication/265332426_Mutual_inductance_calculation_between_misalignment_coils_for_wireless_power_transfer_of_energy
%everything in centimeters

R1 = 7.8232;
R2 = 11.7729;
a_axial = 14.2748;
b_axial = 2.413;
hp = 1.397;
hs = 4.1529;
zc = 7.366;
yc = 30.988;
xc = 0;
N1 = 1142;
N2 = 516;
N = 2;
n = N;
S = 2;
m = S;
h = 0;
q = h;
p = 0;
g = p;
Rp = @(h) R1 + sum(hp*h/(2*N+1));
Rs = @(q) R2 + sum(hs*q/(2*n+1));
alfa = @(q,h) Rs(q)/ Rp(h);

%initial coefficients of the plane of the secondary coil
a = 0;
b = 1;
c = 0.001;
l = @(a, c) sqrt(a^2+c^2);
L = @(a, b, c) sqrt(b^2 + l(a, c)^2);

x = @(p, xc)  xc + b_axial*a*p/(2*m+1);
y = @(p, yc) yc + b_axial*b.*p./(2.*m+1);
z = @(zc, g, p) zc + a_axial.*g./(2.*S+1) + b_axial*c.*p./(2.*m+1);


beta = @(p, xc, h) x(p,xc)./Rp(h);
gamma = @(p, yc,h) y(p, yc)./Rp(h);
delta = @(zc, g, p, h) z(zc, g, p)./Rp(h);
phi = 0;

p1 = @(p, yc, h, a, c) gamma(p,yc, h) .*c/l(a, c);
p2 = @(p, xc, yc, h, a, b, c) (beta(p, xc, h).*l(a, c)^2+gamma(p, yc, h)*a*b)./(l(a,c)*L(a, b, c));
p3 = @(a, b, c) a*c/L(a, b, c);
p4 = @(a,b,c,xc,yc,zc,g,p,h) (gamma(p,yc,h)*l(a,c)^2-beta(p, xc, h)*a-delta(zc, g, p, h)*b*c)/(l(a,c)*L(a,b,c));
p5 = @(a,b,c,xc,zc,g,p,h) (delta(zc,g,p,h)*a-beta(p, xc, h)*c)/l(a,c);

l1 = @(a, b, c) 1 - (b^2*c^2)/(l(a,c)^2*L(a,b,c)^2);
l2 = @(a,b,c) c^2/l(a, c)^2;
l3 = @(a,b,c) a*b*c/(l(a,c)^2*L(a,b,c));

q1 = @(a,b,c,xc,yc,p,h) (gamma(p,yc,h)*l(a,c)^2-beta(p,xc,h)*a*b)/(l(a,c)*L(a,b,c));
q2 = @(a,c,xc,p,h) (-beta(p,xc,h)*c)/l(a,c);

A0 = @(a,b,c,xc,yc,zc,q,g,p,h,phi) 1 + alfa(q,h)^2 + ...
    beta(p,xc,h).^2 + gamma(p,yc,h).^2 + delta(zc,g,p,h).^2 + ... 
    2*alfa(q,h)*(p4(a,b,c,xc,yc,zc,g,p,h).*cos(phi) + p5(a,b,c,xc,zc,g,p,h).*sin(phi));

V0 = @(a,b,c,xc,yc,zc,q,g,p,h,phi) sqrt(beta(p,xc,h).^2 + gamma(p,yc,h).^2 + ... 
    alfa(q,h)^2*(l1(a,b,c)*cos(phi).^2+l2(a,b,c)*sin(phi).^2 + l3(a,b,c)*sin(2*phi)) + ...
    2*alfa(q,h)*(q1(a,b,c,xc,yc,p,h)*cos(phi) + q2(a,c,xc,p,h)*sin(phi)));

k_squared = @(a,b,c,xc,yc,zc,q,g,p,h,phi) 4*V0(a,b,c,xc,yc,zc,q,g,p,h,phi)./(A0(a,b,c,xc,yc,zc,q,g,p,h,phi)+2*V0(a,b,c,xc,yc,zc,q,g,p,h,phi));
k =  sqrt(k_squared(a,b,c,xc,yc,zc,q,g,p,h,phi));

[K,E] = ellipke(k);

Psi = @(a,b,c,xc,yc,zc,q,g,p,h,phi, K, E) (1 - k_squared(a,b,c,xc,yc,zc,q,g,p,h,phi)./2).*K-E;

int_M0 = @(a,b,c,xc,yc,zc,q,g,p,h,phi, K, E) (p1(p, yc, h, a, c).*cos(phi) + ...
    p2(p, xc, yc, h, a, b, c).*sin(phi) + p3(a, b, c)).*Psi(a,b,c,xc,yc,zc,q,g,p,h,phi, K, E)./(k.*sqrt(V0(a,b,c,xc,yc,zc,q,g,p,h,phi).^3));

M0 = @(a,b,c,xc,yc,zc,q,g,p,h,phi, K, E) mu0*Rs(q)./pi .* integral(@(phi) int_M0(a,b,c,xc,yc,zc,q,g,p,h,phi, K, E), 0, 2*pi, 'ArrayValued', true);
disp('M0 = '); disp(M0(a,b,c,xc,yc,zc,q,g,p,h,phi, K, E));

%% 1 test on example 1 from the paper

M = @(a,b,c,xc,yc,zc,q,g,p,h,phi,K,E) N1*N2 .* M0(a,b,c,xc,yc,zc,q,g,p,h,phi,K,E)./(2*n+1)./(2*N+1)./(2*S+1)./(2*m+1);
M_temp=0;
for i = -n:1:n
    q = i;

    for ii = -m:1:m
        p = ii;

       for iii = -N:1:N
            h = iii;
    
            for iv = -S:1:S
                g = iv;
                k = sqrt(abs(k_squared(a,b,c,xc,yc,zc,q,g,p,h,phi)));
                [K,E] = ellipke(k);
                M_temp = M_temp + M(a,b,c,xc,yc,zc,q,g,p,h,phi, K, E);
            end
           
       end  
      
    end 
    
end 
   disp('M_temp = '); disp(M_temp);
   disp('M = '); disp(M(a,b,c,xc,yc,zc,q,g,p,h,phi,K,E));

Answer 1

[David Goodmanson]

Hi Marta,

The code below does not address the algebra in the paper, but I am posting this because it is another way of doing the problem and could be used for comparison purposes. I apologize for yet another system of notation but I think it is fairly clear. (The units are cm, but that is not an indication that the answer is in cgs electromagnetic units, The final answer is still SI, Henries).

The solution is general, including rotation and translation of the second coil compared to the first one. It does use random points ro represent the coils. The method will not work very well if parts of the coils get excessively close to each other compared to their size. If each coil is represented by 10,000 points, then at the end of the calculation there are two matrices whose sum total is 1.6 GBytes of RAM. The memory requirement is proportional to the product of the number of points in each coil. Whatever memory limitations there may be, one can always make several runs and average the results. I didn't multiply by the product of the number of turns in each coil but that is easy enough to do.

mu0 = 4*pi*1e-7;
mu0cm = 4*pi*1e-9;   % distances are in cm

% coil 1 sits in the xy plane (z=0), center at the origin. 
a1 = 2;      % inner radius
b1 = 2.2;    % outer radiius
c1 = .3;     % thickness
meanr1 = mean([a1 b1]);

% coil 2
a2 = 1;      % inner radius
b2 = 1.2;    % outer radius
c2 = .2;     % thickness
meanr2 = mean([a2 b2]);
theta = 30;  % rotation of coil 2 about its center point.
             % rotation is in DEGREES, axis of rotation is the x axis.
x0 = 1;      % center point of coil 2 is translated to x0,y0,z0.
y0 = 1;      % z0 is the height above the plane.
z0 = 4;

n1 = 1e4;    % number of points describing each coil
n2 = 1e4;

% first coil
z1 = (c1/2)*(2*rand(1,n1)-1);
r1 = sqrt(a1^2+(b1^2-a1^2)*rand(1,n1));
phi1 = 2*pi*rand(1,n1);
x1 = r1.*cos(phi1);
y1 = r1.*sin(phi1);
ux1 = -sin(phi1);    % components of unit vector in the direction of current flow
uy1 =  cos(phi1);

% second coil
z2 = (c2/2)*(2*rand(1,n2)-1);
r2 = sqrt(a2^2+(b2^2-a2^2)*rand(1,n2));
phi2 = 2*pi*rand(1,n2);
x2 = r2.*cos(phi2);
y2 = r2.*sin(phi2);
ux2 = -sin(phi2);
uy2 =  cos(phi2);
% second coil is rotated, then translated
x3 = x2 + x0;
y3 = y2*cosd(theta) + z2*sind(theta) + y0;
z3 = z2*cosd(theta) - y2*sind(theta) + z0;
ux3 = -sin(phi2);
uy3 =  cos(phi2)*cosd(theta);
% uz3 = -cos(phi2)*sind(theta);

R = sqrt((x1'-x3).^2 + (y1'-y3).^2 + (z1'-z3).^2);
udot = ux1'*ux3 + uy1'*uy3;
M = (mu0cm/(4*pi))*sum(sum(udot./R))*(2*pi*meanr1/n1)*(2*pi*meanr2/n2)

plot3(x1,y1,z1,'.')
hold on
plot3(x3,y3,z3,'.')
plot3(x0,y0,0,'o')
hold off
xlabel('X')
ylabel('Y')
zlabel('Z')
axis equal

% analytic expression should match when x0=y0=0, theta=0
% and the coils have zero thickness, a1=b1 a2=b2 c1=c2=0 
f1 = meanr1;
f2 = meanr2;
g = z0;
k = 2*sqrt(f1*f2)/sqrt((f1+f2)^2+g^2);
[K E] = ellipke(k^2);
Ma = -mu0cm*sqrt(f1*f2)*((k-2/k)*K + (2/k)*E)
format long
M/Ma
format

ORIGINAL ANSWER

Hi Marta,

You don't show how M is calculated so I will take your word on that, plus the fact that k depends on g.p.h but not q. But one thing that definitely needs to be changed is the M_temp = Mtemp+M statements. You only need the first one, contained in all four for loops, so the code should end

            M_temp = M_temp + M(a,b,c,xc,yc,zc,q,g,p,h,phi, K, E);
         end
      end  
   end 
end  

[Marta Gorecka]

Thanks a lot, @David Goodmanson! I've been trying to sort out my approach, but yours looks solid! Really appreciate your help. Do you have any reference/paper to your code as well?

[David Goodmanson]

Hi Marta,

I don't have a reference for this but it's just a discretization of the Neumann formula

M = mu0/(4*pi) * Int Int (u1.u2) / |x1 - x2| ds1 ds2

where 3d vectors x1 and x2 describe the loops, u1 and u2 are unit vectors in the direction of the current, ds1 and ds2 are arc lengths along the loops and the current density is constant across the cross section of each loop.

In the line that determines M I don't have real good justification for using the arithmetic mean of the inner and outer radii, but experimentally it seems to work.

[Marta Gorecka]

That's a massive help, thank you!