MATLAB Вычисление Интеграла Линий в векторном поле - PullRequest
Новая
       15

MATLAB Вычисление Интеграла Линий в векторном поле

0 голосов
/ 18 октября 2018

ниже приведен код, который создает векторное поле на основе скорости свободного потока и силы вихря в начале координат.

Затем я пытаюсь вычислить циркуляцию путем суммирования площади линейных интегралов,

Example

Подход:

1) Изолировать вектор от его поля скоростей.Пример: line1_2 - горизонтальный вектор из матрицы u, line2 - вертикальный вектор из матрицы u

2) использовать trapz с получением области на линии и с разрешением в качестве фактора. подход

3) Суммируйте все области, и он должен быть равен числителю C для любой замкнутой фигуры, в данном случае 1000. Мы можем поиграть с признаками области 1/2/3/ 4, если это делает его = 1000.

Мой метод не дает правильных значений для суммирования с тем, что должно быть, любые мысли будут очень признательны.

По существу нужно сделать строкуинтеграл для замкнутой формы в поле скоростей или интеграл площади.(теорема Стокса)

Спасибо 

U_i = 10;      % free stream velocity
C = 1000/(2*3.14);% vortex strength
meshfactor = 100;
[x,y] = meshgrid(-U_i:U_i/meshfactor:U_i);

% streamline= -U_i*y + G/(4*3.14)*log(y.^2 + x .^2)
% potential = -U_i*x - G/(2*pi)*atan(y/x)

% u and v components % note: substituted gamma/2i into the variable (line 2)

u = U_i - C*y./(2*(x.^2 + y.^2)); %%define the velocity component in the x
v = C*x./(2*(x.^2 + y.^2)); %%define the velocity component in the y

u(isnan(u)) = 0;
v(isnan(v)) = 0;

%prompt user for input of point coordinates in vector form, ie "[3,3]"
prompt1= 'pt1';
pt1 = input(prompt1);

prompt2= 'pt2';
pt2 = input(prompt2);

prompt3= 'pt3';
pt3 = input(prompt3);

prompt4= 'pt4';
pt4 = input(prompt4);

%% define line 1-2 in the horizontal, goes from left to right
line_1 = pt1(1):U_i/meshfactor:1;pt2(1);
yvec=u(pt1(2),pt1(1):pt2(1));
area_1 = U_i/meshfactor*trapz(yvec);

%% define line 2-3 in the vertical, goes from left to right
line_2 = pt2(2):U_i/meshfactor:1;pt3(2);
yvec2=v(pt2(2):pt3(2),pt2(1));
area_2 = U_i/meshfactor*trapz(yvec2);

%% define line 3-4 in the horizontal, goes from right to left
line_3 = pt4(1):U_i/meshfactor:1;pt3(1);
yvec3=u(pt4(2),pt4(1):pt3(1));
area_3 = U_i/meshfactor*trapz(yvec3);

%% define line 4-1 in the vertical, goes from left to right
line_4 = pt1(2):U_i/meshfactor:1;pt4(2);
yvec4=v(pt1(2):pt4(2),pt1(1));
area_4 = U_i/meshfactor*trapz(yvec4);

circulation = area_1 + area_2 + area_3 + area_4

 figure
 hold on
    title((sprintf(' Unifrom flow with vortex at origin \n Vortex Strength = %0.0f, Free Stream Velocity = %0.0f',C,U_i)));
    quiver(x,y,u,v);
 hold off
...