surface integral of discrete data

41 次查看（过去 30 天）
Peter
Peter 2012-4-22
评论： Sindhu 2014-3-7
I would like to compute the circulation of a velocity field. I think that the best way would be to compute the vorticity and then calculate the surface integral. At the moment I have computed vorticity using curl(X,Y,U,V) Where X,Y,U,V are all 2D matrices. Now that I have vorticity, how can I calculate the surface integral of vorticity? I found solutions if the velocity field can be defined by a function, but not if it is a set of descrete points.
thanks

请先登录，再回答此问题。

回答（2 个）

Richard Brown
Richard Brown 2012-4-22
This came up recently -- in this thread there is code for computing the surface area based on dividing a regular mesh into triangles:
http://www.mathworks.com/matlabcentral/answers/35501-surface-area-from-a-z-matrix
  2 个评论
Peter
Peter 2012-4-23
Thanks for the pointer. I looked over at the other answer and I tried implementing it into my code but i think i'm missing something.
If I consider the a rectangular domain 0<x<3 and 0<y<1 with u = x^2*y and v = x*y I should get -7.5 for the circulation around the rectangle
However I'm getting a value of ~5.8. If you could help clarify I would greatly appreciate it. my code is below
function green_test
clear all
close all
clc
dx = 0.01;
dy = 0.01;
x = 0:dx:3;
y = 0:dy:1;
[X,Y] = meshgrid(x,y);
U = X.^2.*Y;
V = X.*Y;
vort = curl(X,Y,U,V);
a = surface_int(dx,dy,vort,-5)
function [surfaceArea] = surface_int(dx,dy,Z,cutoff)
[m, n] = size(Z);
areas = 0.5*sqrt((dx*dy)^2 + (dx*(Z(1:m-1,2:n) - Z(1:m-1,1:n-1))).^2 + ...
(dy*(Z(2:m,1:n-1) - Z(1:m-1,1:n-1))).^2) + ...
0.5*sqrt((dx*dy)^2 + (dx*(Z(1:m-1,2:n) - Z(2:m,2:n))).^2 + ...
(dy*(Z(2:m,1:n-1) - Z(2:m,2:n))).^2);
zMean = 0.25 * (Z(1:m-1,1:n-1) + Z(1:m-1,2:n) + Z(2:m,1:n-1) + Z(2:m,2:n));
areas(zMean <= cutoff) = 0;
surfaceArea = sum(areas(:));
sprintf('Total surface area is %2.4f\n', surfaceArea)
Richard Brown
Richard Brown 2012-4-23
Yeah, your question is actually easier than that - see my other answer

请先登录，再进行评论。

Richard Brown
Richard Brown 2012-4-23
Sorry, I read your question too fast and pointed you to code for surface area, not surface integral. In this case it's even easier, as your "surface" that you're integrating over is just the xy plane, and the vorticity is all oriented directly upwards. So the surface integral is simply
dx = 0.01;
dy = 0.01;
x = 0:dx:3;
y = 0:dy:1;
[X,Y] = meshgrid(x,y);
U = X.^2.*Y;
V = X.*Y;
[W,~] = curl(X,Y,U,V);
[m,n] = size(W);
vort = 0.25 * (W(1:m-1,1:n-1) + W(1:m-1,2:n) + W(2:m,1:n-1) + W(2:m,2:n));
circulation = sum(vort(:))*dx*dy
A couple of notes:
  1. Make sure you call curl with two output arguments, otherwise it returns the angular velocity only
  2. The vorticity that you're integrating over is evaluated at the centre of each rectangle, hence the averaging expression on the vort = line
  2 个评论
Richard Brown
Richard Brown 2012-4-23
And you could do it without the averaging by defining your grid slightly differently
x = (dx/2):dx:(3-dx/2);
y = (dy/2):dy:(1-dy/2);
circulation = sum(W(:)) * dx * dy
Sindhu
Sindhu 2014-3-7
I also have the same problem statement. "To find Circulation in the yz plane". I understand the coding till the vorticity part. But why do you take an averaged value for vorticity? This is hopefully first order accurate and what should I do If I need second order accuracy? Please suggests me some books or sites to understand the calculations. I went through few lectures on surface integral and parameterization of surfaces. It seems to be confusing and I could not figure out what is the math behind it.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Fluid Dynamics 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by