sum of a sum of a sum of a sum

8 views (last 30 days)
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
  1 Comment
Marta Gorecka
Marta Gorecka on 28 Feb 2022
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));

Sign in to comment.

Accepted Answer

David Goodmanson
David Goodmanson on 3 Mar 2022
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
  3 Comments
Marta Gorecka
Marta Gorecka on 21 Mar 2022
That's a massive help, thank you!

Sign in to comment.

More Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!