surface integral of discrete data

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

Answer 1

[Richard Brown]

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

[Peter]

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]

Yeah, your question is actually easier than that - see my other answer

Answer 2

[Richard Brown]

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

[Richard Brown]

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]

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.