da D4lF4zZI0 » 21/05/2019, 07:45
Infatti, lo scopo è estendere l'algoritmo anche a geometria non semplici.
Comunque, in questi giorni ho buttato giù un pò di codice ( che credo che funzioni, ma una consulenza e un consiglio sarebbe gradito ) che incollo qua:
clf
clear all
clc
%% Dominio
% Numero di punti del dominio
ND = 10;
% Estensione del dominio
Dom = [-1 1;
-1 1;
0.1 6];
%% Geometria di analisi
% Parametri
I = 20;
mu = 4*pi*10^-7;
gamma = mu*I/(4*pi);
% Passo di integrazione
ds = 1e-4;
% Geometria filo rettilineo ( supposto disposto lungo l'asse z )
L = [0 0 0;
0 0 5];
Nl = numel(L)/3;
%% Inizializzazione variabili
% Componenti campo di induzione magnetica B = (Bx, By, Bz);
Bx = zeros(ND, ND, ND);
By = zeros(ND, ND, ND);
Bz = zeros(ND, ND, ND);
% Volume Mesh
[X, Y, Z] = meshgrid(linspace(Dom(1,1), Dom(1,2), ND), ...
linspace(Dom(2,1), Dom(2,2), ND), ...
linspace(Dom(3,1), Dom(3,2), ND));
%% Integrazione legge di Biot-Savart
for i = 1:ND
for j = 1:ND
for k = 1:ND
% Punto in cui calcolare il campo
pTest = [X(i,j,k) Y(i,j,k) Z(i,j,k)];
% La curva viene discretizzata in Nl punti e l'argoritmo
% viene iterato Nl-1 volte ( Infatti sono Nl-1 i segmenti in
% cui viene divisa la curva.
% Ogni segmento, a sua volta, viene suddiviso in passi elementi elementari
% di lunghezza ds ( passo di integrazione ) e per ognuno degli
% elementi viene calcolato l'incremento del campo dB
for pCurv = 1:Nl-1
% Lunghezza dl
len = norm(L(pCurv,:) - L(pCurv+1,:));
% Numero di punti dl
% Viene effettuato un contollo della discretizzazione della
% curva; infatti, se dl < ds, allora non avrebbe senso la
% discretizzazione.
Npi = ceil(len/ds);
if Npi < 3
error('Il passo di integrazione è troppo grande')
end
% Discretizzazione curva
Lx = linspace(L(pCurv,1), L(pCurv+1,1), Npi);
Ly = linspace(L(pCurv,2), L(pCurv+1,2), Npi);
Lz = linspace(L(pCurv,3), L(pCurv+1,3), Npi);
% Integrazione
for s = 1:Npi-1
Rx = Lx(s) - pTest(1);
Ry = Ly(s) - pTest(2);
Rz = Lz(s) - pTest(3);
dLx = Lx(s+1) - Lx(s);
dLy = Ly(s+1) - Ly(s);
dLz = Lz(s+1) - Lz(s);
% Distanza
R = sqrt(Rx^2 + Ry^2 + Rz^2);
% Biot-Savart
dL = [dLx dLy dLz];
r = [Rx Ry Rz];
dB = gamma*cross(dL,r)/R^3;
% Somma ( integrale ) del campo
Bx(i,j,k) = Bx(i,j,k) + dB(1);
By(i,j,k) = By(i,j,k) + dB(2);
Bz(i,j,k) = Bz(i,j,k) + dB(3);
end
end
end
end
end
%% Modulo del campo di induzione magnetica
M = sqrt(Bx.^2+By.^2+Bz.^2);
%% Grafico
quiver3(X,Y,Z,Bx,By,Bz);
hold on
plot3(L(:,1),L(:,2),L(:,3),'r-','linewidth',3);
hold on
xlabel('x');
ylabel('y');
zlabel('z');
title('Campo vettoriale induzione magnetica');