Thread Subject:

how to integrate a vector over a discrete surface?

I have a surface S in discrete from, i.g., i only have [x,y,z]data, but not a function. This surface is under air pressure P, which is perpendicular to the surface everywhere with a constant value. I want to know the force vector F on the surface caused by the air pressure, i.e. integrate P over S in discrete form. Could anybody help? How about surface S is in polar format? Thanks a lot.

"yl " <shanshicn@hotmail.com> wrote in message <il3k7n$dqu$1@fred.mathworks.com>...
> I have a surface S in discrete from, i.g., i only have [x,y,z]data, but not a function. This surface is under air pressure P, which is perpendicular to the surface everywhere with a constant value. I want to know the force vector F on the surface caused by the air pressure, i.e. integrate P over S in discrete form. Could anybody help? How about surface S is in polar format? Thanks a lot.
- - - - - - - - - - - -
  With constant pressure it is not necessary to integrate over the surface to get the total force. That force depends only on the boundary of the surface and the constant pressure P. This principle is well-known in rocket propulsion theory.

  Just take the cross product line integral

 1/2*int r cross dr

around the surface boundary and then multiply the result by P. The 'r' here is a variable vector pointing from some fixed point such as the coordinate origin to points along the boundary curve. The line integration must proceed forward with the surface always off to the left. (Of course if you have a closed surface this force will be zero.)

  For example suppose you have a balloon as your surface and the balloon surface boundary (the open end) is a circle in the x-y plane below the balloon with unit radius. Then r = cos(t)*i + sin(t)*j where i, j and k are the usual x, y, and z unit vectors and t is the angle with respect to the x-axis. The line integral as t varies through a full circle would be

 1/2*int (r cross dr) =
 1/2*int (cos(t)*i+sin(t)*j) cross (-sin(t)*i+cos(t)*j) dt =
 1/2*int (cos(t)^2+sin(t)^2)*k dt = 1/2*int k dt = pi*k

When multiplied by P you would get pi*P*k which is a force along the z-axis of pi*P. This is just the circle's area times P. That is, the force is the same as if the balloon were just a flat diaphragm stretched across the unit circle.

Roger Stafford

Roger,Thanks a lot. That is very helpful. Just to understand better about the principle, any recommend readings? thanks.

"Roger Stafford" wrote in message <il459e$opn$1@fred.mathworks.com>...
> "yl " <shanshicn@hotmail.com> wrote in message <il3k7n$dqu$1@fred.mathworks.com>...
> > I have a surface S in discrete from, i.g., i only have [x,y,z]data, but not a function. This surface is under air pressure P, which is perpendicular to the surface everywhere with a constant value. I want to know the force vector F on the surface caused by the air pressure, i.e. integrate P over S in discrete form. Could anybody help? How about surface S is in polar format? Thanks a lot.
> - - - - - - - - - - - -
> With constant pressure it is not necessary to integrate over the surface to get the total force. That force depends only on the boundary of the surface and the constant pressure P. This principle is well-known in rocket propulsion theory.
>
> Just take the cross product line integral
>
> 1/2*int r cross dr
>
> around the surface boundary and then multiply the result by P. The 'r' here is a variable vector pointing from some fixed point such as the coordinate origin to points along the boundary curve. The line integration must proceed forward with the surface always off to the left. (Of course if you have a closed surface this force will be zero.)
>
> For example suppose you have a balloon as your surface and the balloon surface boundary (the open end) is a circle in the x-y plane below the balloon with unit radius. Then r = cos(t)*i + sin(t)*j where i, j and k are the usual x, y, and z unit vectors and t is the angle with respect to the x-axis. The line integral as t varies through a full circle would be
>
> 1/2*int (r cross dr) =
> 1/2*int (cos(t)*i+sin(t)*j) cross (-sin(t)*i+cos(t)*j) dt =
> 1/2*int (cos(t)^2+sin(t)^2)*k dt = 1/2*int k dt = pi*k
>
> When multiplied by P you would get pi*P*k which is a force along the z-axis of pi*P. This is just the circle's area times P. That is, the force is the same as if the balloon were just a flat diaphragm stretched across the unit circle.
>
> Roger Stafford

Roger, one more question. which integration command in matlab for the integration on discrete boundary? thanks a lot.

"Roger Stafford" wrote in message <il459e$opn$1@fred.mathworks.com>...
> "yl " <shanshicn@hotmail.com> wrote in message <il3k7n$dqu$1@fred.mathworks.com>...
> > I have a surface S in discrete from, i.g., i only have [x,y,z]data, but not a function. This surface is under air pressure P, which is perpendicular to the surface everywhere with a constant value. I want to know the force vector F on the surface caused by the air pressure, i.e. integrate P over S in discrete form. Could anybody help? How about surface S is in polar format? Thanks a lot.
> - - - - - - - - - - - -
> With constant pressure it is not necessary to integrate over the surface to get the total force. That force depends only on the boundary of the surface and the constant pressure P. This principle is well-known in rocket propulsion theory.
>
> Just take the cross product line integral
>
> 1/2*int r cross dr
>
> around the surface boundary and then multiply the result by P. The 'r' here is a variable vector pointing from some fixed point such as the coordinate origin to points along the boundary curve. The line integration must proceed forward with the surface always off to the left. (Of course if you have a closed surface this force will be zero.)
>
> For example suppose you have a balloon as your surface and the balloon surface boundary (the open end) is a circle in the x-y plane below the balloon with unit radius. Then r = cos(t)*i + sin(t)*j where i, j and k are the usual x, y, and z unit vectors and t is the angle with respect to the x-axis. The line integral as t varies through a full circle would be
>
> 1/2*int (r cross dr) =
> 1/2*int (cos(t)*i+sin(t)*j) cross (-sin(t)*i+cos(t)*j) dt =
> 1/2*int (cos(t)^2+sin(t)^2)*k dt = 1/2*int k dt = pi*k
>
> When multiplied by P you would get pi*P*k which is a force along the z-axis of pi*P. This is just the circle's area times P. That is, the force is the same as if the balloon were just a flat diaphragm stretched across the unit circle.
>
> Roger Stafford

"yl " <shanshicn@hotmail.com> wrote in message <il5gld$2oh$1@fred.mathworks.com>...
> Roger,Thanks a lot. That is very helpful. Just to understand better
> about the principle, any recommend readings? thanks.

  I'm not sure what reading to recommend to you. (My own sources are very ancient.) I'll show you an intuitive derivation of the formula I gave and that should at least be an indication of where in more modern literature to look.

  The line integral formula "F = P/2*int (r cross dr)" is a consequence of the famous Green's Theorem. See:

 http://en.wikipedia.org/wiki/Green's_theorem

  If S is the surface you are integrating the constant pressure, P, over, let ds be an infinitesimal normal vector area element. Then the orthogonal projection of P*ds onto the x-y plane would just be the infinitesimal z-component of the force exerted on that surface element. Therefore the total z-component of the force would simply be the projected area of S onto the x-y plane times P which is bounded by the similar projection of the S boundary there. (If the projected surface bulges past this projected boundary there will always be an appropriate cancellation of that portion, so we only need consider the interior of the projected boundary.)

  The z-component of the formula's "P/2*int (r cross dr)" will be

 P/2*int (x*dy-y*dx)

integrated along the projected boundary. By Green's Theorem using their L = -y and M = x (see the above website) this would be equal to

 P/2*int (dM/dx-dL/dy) dxdy = P/2*int (1+1) dxdy = P*int 1 dxdy

integrated over the projected x-y area, which will give us just the projected (interior) area times P, and therefore the z-component of the total force.

  Applied to all three x, y, and z components, this justifies the above formula.

> Roger, one more question. which integration command in matlab for the
> integration on discrete boundary? thanks a lot.

  That question is not so easy to answer. In your original x,y,z data, to do a proper surface integration you would have had to know accurately where the boundary of the surface lay. With your data at scattered interior discrete points, it would not be evident just where this was.

  Presuming you somehow obtain a sequence of sufficiently closely-spaced successive points lying right along this required boundary, you could use the 'trapz' function for each of the three components of P/2*int (r cross dr) to find the three force components. The remaining data points of the interior would not be needed at all. If such boundary points were not sufficiently close together to achieve the desired accuracy, you might have to resort to some interpolation to provide more closely-spaced points for 'trapz' or perhaps you could use one of the higher order integration schemes for discrete points present on the file exchange. I wrote one at:

 http://www.mathworks.com/matlabcentral/fileexchange/19152

which does third-order cumulative integration of a sequence of discrete points, and there are a number of others to be found there.

Roger Stafford

"Roger Stafford" wrote in message <il6873$mhm$1@fred.mathworks.com>...
> ......
> The line integral formula "F = P/2*int (r cross dr)" is a consequence of the famous Green's Theorem. See:
> .....
 - - - - - - - -
  There are a couple of things to mention about using 'trapz' to integrate around the boundary of your surface.

  First, the reasoning using Green's theorem I went through earlier also can be used to show that

 P*[int y*dz, int z*dx, int x*dy]

can be used in place of

 P/2*int r cross dr =
 P/2*[int y*dz-z*dy, int z*dx-x*dz, int x*dy-y*dx]

which I described before, and this reduces the necessary computation by half.

  Second, it is necessary that in the vectors used by 'trapz' for this purpose the first element must be the same as the last element, so that the path of integration is closed. If this is not the case, you must augment the vectors by repeating the first elements at the vectors' ends.

  The expression for the force using this shorter form would then be:

 F = P*[trapz(z,y),trapz(x,z),trapz(y,x)]

where x, y, and z are such "closed" vectors around the boundary path (advancing with the surface to the left.)

Roger Stafford

Roger,

Thanks a lot for your detailed explanation. I got the boundary data. instead of trapez, i multiply force component with segment length between two neighbor points and add them together. I repeated the integration 3 times for the 3 directions. I may need to check my method using one simple case. Anyway, thanks a lot for remind me the Green's theorem, which i totally forgot.

If the balloon holding the pressure is not uniform, say, some part is thicker than others, the Green's theory and your analysis should still hold as long as the pressure is constant. Am I right? thanks.

"Roger Stafford" wrote in message <il952f$kmm$1@fred.mathworks.com>...
> "Roger Stafford" wrote in message <il6873$mhm$1@fred.mathworks.com>...
> > ......
> > The line integral formula "F = P/2*int (r cross dr)" is a consequence of the famous Green's Theorem. See:
> > .....
> - - - - - - - -
> There are a couple of things to mention about using 'trapz' to integrate around the boundary of your surface.
>
> First, the reasoning using Green's theorem I went through earlier also can be used to show that
>
> P*[int y*dz, int z*dx, int x*dy]
>
> can be used in place of
>
> P/2*int r cross dr =
> P/2*[int y*dz-z*dy, int z*dx-x*dz, int x*dy-y*dx]
>
> which I described before, and this reduces the necessary computation by half.
>
> Second, it is necessary that in the vectors used by 'trapz' for this purpose the first element must be the same as the last element, so that the path of integration is closed. If this is not the case, you must augment the vectors by repeating the first elements at the vectors' ends.
>
> The expression for the force using this shorter form would then be:
>
> F = P*[trapz(z,y),trapz(x,z),trapz(y,x)]
>
> where x, y, and z are such "closed" vectors around the boundary path (advancing with the surface to the left.)
>
> Roger Stafford

"yl " <shanshicn@hotmail.com> wrote in message <iloh37$72n$1@ginger.mathworks.com>...
> Roger,
>
> Thanks a lot for your detailed explanation. I got the boundary data. instead of trapez, i multiply force component with segment length between two neighbor points and add them together. I repeated the integration 3 times for the 3 directions. I may need to check my method using one simple case. Anyway, thanks a lot for remind me the Green's theorem, which i totally forgot.
>
> If the balloon holding the pressure is not uniform, say, some part is thicker than others, the Green's theory and your analysis should still hold as long as the pressure is constant. Am I right? thanks.
- - - - - - - - -
  Your description "multiply force component with segment length" worries me a bit, yl. The x-component line integral quantities

 int 1/2*(y*dz-z*dy)

or

 int y*dz

around the boundary are precisely what you would use if you projected the boundary orthogonally onto the y-z plane and then proceeded to find the area in that plane within the resulting curve (making due allowances for the curve possibly crossing itself.) The same is true of the other two line integral components.

  In terms of the quantity

 int 1/2*(r cross dr)

you should be taking the vector cross product of the vector r pointing from the origin by a vector dr along an infinitesimal line segment of the boundary, and the result would be orthogonal to both these two vectors. The vector sum of all these vector cross products times P will give you the required net force.

  If your boundary is actually a polygonal series of straight line segments, this can be expressed precisely as the finite sum of such cross products and there would be no need for an integration routine.

  The pressure P is merely a scalar factor in all of this.

  Is what you are doing equivalent to the above?

  In answer to your second paragraph question, the net force you are calculating has nothing whatever to do with the container (balloon) or its thickness as long as neither it nor the air inside are moving. When the air begins to escape out the orifice and creates a pressure gradient inside, the situation will of course change.

Roger Stafford

Roger, thanks for your reminding. i may not describe correctly last time and here is what I did:

1) define vector R (x,y,z) for each point(x,y,z)
2) define dR=diff(R)
3) Fi=cross(R,dR) and sum together for all points to get F
4) F has three components in x, y, z directions respectively, i.e., force components

    Rv0=[Rd.*cos(Td),Rd.*sin(Td),Zd];%transfer from polar system to (x,y,z)
    Rv=Rv0;
    dR=diff(Rv);
    Rv(1,:)=[];
    s=0;
    for i=1:length(dT)
        s=s+cross(Rv(i,:),dR(i,:));
    end
s is the force vector include components in x,y,z direcsions. Please let me know if I was wrong. thanks.



"Roger Stafford" wrote in message <ilon09$o2h$1@ginger.mathworks.com>...
> "yl " <shanshicn@hotmail.com> wrote in message <iloh37$72n$1@ginger.mathworks.com>...
> > Roger,
> >
> > Thanks a lot for your detailed explanation. I got the boundary data. instead of trapez, i multiply force component with segment length between two neighbor points and add them together. I repeated the integration 3 times for the 3 directions. I may need to check my method using one simple case. Anyway, thanks a lot for remind me the Green's theorem, which i totally forgot.
> >
> > If the balloon holding the pressure is not uniform, say, some part is thicker than others, the Green's theory and your analysis should still hold as long as the pressure is constant. Am I right? thanks.
> - - - - - - - - -
> Your description "multiply force component with segment length" worries me a bit, yl. The x-component line integral quantities
>
> int 1/2*(y*dz-z*dy)
>
> or
>
> int y*dz
>
> around the boundary are precisely what you would use if you projected the boundary orthogonally onto the y-z plane and then proceeded to find the area in that plane within the resulting curve (making due allowances for the curve possibly crossing itself.) The same is true of the other two line integral components.
>
> In terms of the quantity
>
> int 1/2*(r cross dr)
>
> you should be taking the vector cross product of the vector r pointing from the origin by a vector dr along an infinitesimal line segment of the boundary, and the result would be orthogonal to both these two vectors. The vector sum of all these vector cross products times P will give you the required net force.
>
> If your boundary is actually a polygonal series of straight line segments, this can be expressed precisely as the finite sum of such cross products and there would be no need for an integration routine.
>
> The pressure P is merely a scalar factor in all of this.
>
> Is what you are doing equivalent to the above?
>
> In answer to your second paragraph question, the net force you are calculating has nothing whatever to do with the container (balloon) or its thickness as long as neither it nor the air inside are moving. When the air begins to escape out the orifice and creates a pressure gradient inside, the situation will of course change.
>
> Roger Stafford

maybe the code should be:
    Rv0=[Rd.*cos(Td),Rd.*sin(Td),Zd];
    Rv=Rv0;
    dR=diff(Rv);
    Rv(1,:)=[];
    dT=diff(Td);%angle
    s=0;
    for i=1:length(dT)
        s=s+cross(Rv(i,:),dR(i,:)).*dT(i);
    end

"yl " <shanshicn@hotmail.com> wrote in message <ilt2vo$7ms$1@ginger.mathworks.com>...
> Roger, thanks for your reminding. i may not describe correctly last time and here is what I did:
>
> 1) define vector R (x,y,z) for each point(x,y,z)
> 2) define dR=diff(R)
> 3) Fi=cross(R,dR) and sum together for all points to get F
> 4) F has three components in x, y, z directions respectively, i.e., force components
>
> Rv0=[Rd.*cos(Td),Rd.*sin(Td),Zd];%transfer from polar system to (x,y,z)
> Rv=Rv0;
> dR=diff(Rv);
> Rv(1,:)=[];
> s=0;
> for i=1:length(dT)
> s=s+cross(Rv(i,:),dR(i,:));
> end
> s is the force vector include components in x,y,z direcsions. Please let me know if I was wrong. thanks.
>
>
>
> "Roger Stafford" wrote in message <ilon09$o2h$1@ginger.mathworks.com>...
> > "yl " <shanshicn@hotmail.com> wrote in message <iloh37$72n$1@ginger.mathworks.com>...
> > > Roger,
> > >
> > > Thanks a lot for your detailed explanation. I got the boundary data. instead of trapez, i multiply force component with segment length between two neighbor points and add them together. I repeated the integration 3 times for the 3 directions. I may need to check my method using one simple case. Anyway, thanks a lot for remind me the Green's theorem, which i totally forgot.
> > >
> > > If the balloon holding the pressure is not uniform, say, some part is thicker than others, the Green's theory and your analysis should still hold as long as the pressure is constant. Am I right? thanks.
> > - - - - - - - - -
> > Your description "multiply force component with segment length" worries me a bit, yl. The x-component line integral quantities
> >
> > int 1/2*(y*dz-z*dy)
> >
> > or
> >
> > int y*dz
> >
> > around the boundary are precisely what you would use if you projected the boundary orthogonally onto the y-z plane and then proceeded to find the area in that plane within the resulting curve (making due allowances for the curve possibly crossing itself.) The same is true of the other two line integral components.
> >
> > In terms of the quantity
> >
> > int 1/2*(r cross dr)
> >
> > you should be taking the vector cross product of the vector r pointing from the origin by a vector dr along an infinitesimal line segment of the boundary, and the result would be orthogonal to both these two vectors. The vector sum of all these vector cross products times P will give you the required net force.
> >
> > If your boundary is actually a polygonal series of straight line segments, this can be expressed precisely as the finite sum of such cross products and there would be no need for an integration routine.
> >
> > The pressure P is merely a scalar factor in all of this.
> >
> > Is what you are doing equivalent to the above?
> >
> > In answer to your second paragraph question, the net force you are calculating has nothing whatever to do with the container (balloon) or its thickness as long as neither it nor the air inside are moving. When the air begins to escape out the orifice and creates a pressure gradient inside, the situation will of course change.
> >
> > Roger Stafford

"yl " <shanshicn@hotmail.com> wrote in message <ilt2vo$7ms$1@ginger.mathworks.com>...
> Roger, thanks for your reminding. i may not describe correctly last time and here is what I did:
>
> 1) define vector R (x,y,z) for each point(x,y,z)
> 2) define dR=diff(R)
> 3) Fi=cross(R,dR) and sum together for all points to get F
> 4) F has three components in x, y, z directions respectively, i.e., force components
>
> Rv0=[Rd.*cos(Td),Rd.*sin(Td),Zd];%transfer from polar system to (x,y,z)
> Rv=Rv0;
> dR=diff(Rv);
> Rv(1,:)=[];
> s=0;
> for i=1:length(dT)
> s=s+cross(Rv(i,:),dR(i,:));
> end
>
- - - - - - - - - - - -
  Yes that looks right, provided that the first and last values in Rd, Td, and Zd represent the same point and provided that length(dT) is the same as length(dR).

  You could just write:

 s = sum(cross(Rv,dR,2),1);

instead of the for-loop.

  Your second version isn't correct. You don't want dT(i) in there. That would make the answer a scalar and moreover vanishingly small with both dR and dT multiplied together.

  You can test your answer by running the boundary around a flat circle and see if your answer is (close to) the area of the circle.

Roger Stafford

Roger,

You are right and I have tested the codes on a flat circle. One question: does Green's theory or the code requires the boundary in one plane? Looks like it is true according to the assumptions of GREEN'S THEOREM. (1827) If F(x; y) = (P(x; y);Q(x; y)) is a vector field in the plane and R is a region
in the plane which has as a boundary a piecewise smooth closed curve
 traversed in the direction so that the
region R is "to the left"...

For 3d case, if the boundary is not on one plane, will the same conclusion holds? thanks.

"Roger Stafford" wrote in message <iltbhd$726$1@ginger.mathworks.com>...
> "yl " <shanshicn@hotmail.com> wrote in message <ilt2vo$7ms$1@ginger.mathworks.com>...
> > Roger, thanks for your reminding. i may not describe correctly last time and here is what I did:
> >
> > 1) define vector R (x,y,z) for each point(x,y,z)
> > 2) define dR=diff(R)
> > 3) Fi=cross(R,dR) and sum together for all points to get F
> > 4) F has three components in x, y, z directions respectively, i.e., force components
> >
> > Rv0=[Rd.*cos(Td),Rd.*sin(Td),Zd];%transfer from polar system to (x,y,z)
> > Rv=Rv0;
> > dR=diff(Rv);
> > Rv(1,:)=[];
> > s=0;
> > for i=1:length(dT)
> > s=s+cross(Rv(i,:),dR(i,:));
> > end
> >
> - - - - - - - - - - - -
> Yes that looks right, provided that the first and last values in Rd, Td, and Zd represent the same point and provided that length(dT) is the same as length(dR).
>
> You could just write:
>
> s = sum(cross(Rv,dR,2),1);
>
> instead of the for-loop.
>
> Your second version isn't correct. You don't want dT(i) in there. That would make the answer a scalar and moreover vanishingly small with both dR and dT multiplied together.
>
> You can test your answer by running the boundary around a flat circle and see if your answer is (close to) the area of the circle.
>
> Roger Stafford

ok, for higher dimensions, there is an expansion of Green's theorem: Stokes' Theorem. So the codes should work well for 3d cases.

"yl " <shanshicn@hotmail.com> wrote in message <in2i95$8fp$1@fred.mathworks.com>...
> Roger,
>
> You are right and I have tested the codes on a flat circle. One question: does Green's theory or the code requires the boundary in one plane? Looks like it is true according to the assumptions of GREEN'S THEOREM. (1827) If F(x; y) = (P(x; y);Q(x; y)) is a vector field in the plane and R is a region
> in the plane which has as a boundary a piecewise smooth closed curve
> traversed in the direction so that the
> region R is "to the left"...
>
> For 3d case, if the boundary is not on one plane, will the same conclusion holds? thanks.
>
> "Roger Stafford" wrote in message <iltbhd$726$1@ginger.mathworks.com>...
> > "yl " <shanshicn@hotmail.com> wrote in message <ilt2vo$7ms$1@ginger.mathworks.com>...
> > > Roger, thanks for your reminding. i may not describe correctly last time and here is what I did:
> > >
> > > 1) define vector R (x,y,z) for each point(x,y,z)
> > > 2) define dR=diff(R)
> > > 3) Fi=cross(R,dR) and sum together for all points to get F
> > > 4) F has three components in x, y, z directions respectively, i.e., force components
> > >
> > > Rv0=[Rd.*cos(Td),Rd.*sin(Td),Zd];%transfer from polar system to (x,y,z)
> > > Rv=Rv0;
> > > dR=diff(Rv);
> > > Rv(1,:)=[];
> > > s=0;
> > > for i=1:length(dT)
> > > s=s+cross(Rv(i,:),dR(i,:));
> > > end
> > >
> > - - - - - - - - - - - -
> > Yes that looks right, provided that the first and last values in Rd, Td, and Zd represent the same point and provided that length(dT) is the same as length(dR).
> >
> > You could just write:
> >
> > s = sum(cross(Rv,dR,2),1);
> >
> > instead of the for-loop.
> >
> > Your second version isn't correct. You don't want dT(i) in there. That would make the answer a scalar and moreover vanishingly small with both dR and dT multiplied together.
> >
> > You can test your answer by running the boundary around a flat circle and see if your answer is (close to) the area of the circle.
> >
> > Roger Stafford