To represent a surface in three space we paramaterize each of the three coordinates with two arbitrary parameters. For a given paramaterization i,j there is a corresponding point in three space P(xi,j,yi,j,zi,j).
On a computer we cannot store the infinite number of points that make up a surface. As a consequence we chose a finite set of discrete values. If we are describing a sphere in spherical coordinates, we may choose values of theta that are linearly spaced between 0 and 2*pi and values of phi that are linearly spaced between 0 and pi. (Note r(theta,phi)=Ro for the unstrained shpere).
The Cartesian product:
{theta}x{phi}x{r}={(theta,phi,r)| theta=(i*2*pi/n), phi=(j*pi/n), r=rsphere, i,j=1,2,.....n}
defines all points on the surface and each point is paramterized by the set of indices i,j. To represent all these points, we place each coordinate in a matrix where the value of coordanate that corresponds to i=i1 and j=j1 is found in row i1 and column j1.
Figure 1. This is Matrix representing a single coordinate. To describe the points in three space three of these matrices are used (one for each coordinate). The coordantes for a single point which is paramaterized by i,j are found in the ith row and jth column of each matricies.
It is often convenient when plotting these points to transform them into rectangular coordinates. If phi is the angle between z and the position vector and theta is the rotation of the position vector about the z axis from the x-z plane then we can use the transformations:
(1)
0<2 0
When the computer draws the surface it draws a mesh of three dimensional curves. One way to draw the curves is with line segments. When the computer draws the curve with line segments it connects P(xi,j,yi,j,zi,j) to P(xi+1,j,yi+1,j,zi+1,j) and P(xi,j,yi,j,zi,j) to P(xi,j+1,yi,j+1,zi,j+1). In the matlab code bellow a sphere is defined in spherical coordinates, it is converted to rectangular coordinates and then plotted. (For the program we will describe later the surface was created in a function. This function can be found in appendix A..1)
n=n-1
for i=1:n+1
for j=1:n+1
theta(i,j)=(2*pi*(i-1))/(n);
phi(i,j)=(pi*(j-1))/(n);
end
end
n=n+1
%Transform to Cartesian
x=r.*sin(phi).*cos(theta);
y=r.*sin(phi).*sin(theta);
z=r.*cos(phi);
% Plot Surface
mesh(x,y,z);
% Label axies
xlabel('x axis');
ylabel('y axis');
zlabel('z axis');
title('Figure 3. A Mesh Plot of a Sphere');
Closing Curves constant in theta
In the above code we were carful to close the curves. Curves wich are constant in theta (e.g lines of latitude) can be closed by connecting curves that differ in theta by 180 degrees. This was done by ensuring we included the end points (i.e. phi=0 and phi=pi) for each curve. Connecting curves that are constant in phi can be cone by connecting the end point to the beginning point. This was done by having the first row of the coordinate matricies contain points where theta=0 and the last row contain points where theta=2*pi.
The area Vectors.
The mesh divides the sphere into sections that are approximately rectangular. In the limiting case where the number of mesh points approaches infinity the average deviation from a square shape of these sections will approaches zero. The cross product of a two distinct vectors going from the same mesh point to an adjacent mesh point will give a vector whose area is approximately the area of on of an adjacent section and the direction of the vector will be approximately perpendicular to the surface. The average deviation in the direction of these vectors from a normal vector to the surface will approach zero as the number of mesh points approaches zero.
Examining the surface and using the right-hand rule it is clear that one cross product that will give a vector perpendicular to and pointing away from the surface is {P(xi,j,yi,j,zi,j) to P(xi,j+1,yi,j+1,zi,j+1)} cross {P(xi,j,yi,j,zi,j) to P(xi+1,j,yi+1,j,zi+1,j)}
The matlab function bellow takes a surface given by the three nxn matrices x,y,z and returns the area vectors for the sections gridded off by lines constant in i and lines constant in j.
We are able to find the volume of the balloon by integrating in spherical coordinates without actually converting to spherical coordinates. The surface of the sphere is divided up into a mesh. We can divide the volume of the sphere up into pyramid's where the apex of the pyramid is at the origin and the base of the pyramid has a mesh point at each corner. Choosing as the base the sections that are on average are approximately square we can find the area of the base as was discussed above. The height of the pyramid is the shortest distance from the apex to a point in the plane of the pyramid's base. This will be the magnitude of: the scalar projection of a vector that runs from the apex of the pyramid to one of the corners of the base of the pyramid onto a vector that is perpendicular to the plane of the base of the pyramid.
If u and v are vectors that run along the base of the pyramid and rx,y,z is a vector that runs from the origin to one of the corners of the base then we can use the following relations:
Volume=(1/3)A*h (2)
A=|uxv| (3)
h=|projuxv
rx,y,z| (4)
When uxv points away from the center of the sphere we can use the following relation:
h=compuxv rx,y,z (5)
combining (2) (3) and (5) we get:
(6)
(7)
The following Matlab code computes the volume of each pyramid then sums the volume of all pyramids together to get the total volume.
With the matlab code above and n=20 The approximation yields a value of: 4.0848 and the actual volume of a sphere with radius 1 is: 4.1887.
{P(xi,j,yi,j,zi,j) to P(xi,j+1,yi,j+1,zi,j+1)} cross {P(xi,j,yi,j,zi,j) to P(xi+1,j,yi+1,j,zi+1,j)}
We recall that this produces an area vector that points away from the surface. The average pressure force in this piece of area will be equal to the area vector multiplied by: the pressure inside the balloon minus the pressure outside the balloon.
(8)
It is more convenient for us to know the pressure force of a section of area for wich a mesh point is at the center. This section will have approximately the same area. As a consequence the magnitude of the average force acting on this section of area will be approximately the same.
However when we consider its direction we must be careful. Because in our matrices the area vector we assigned to the mesh points was from the section in the direction of increasing i and j. If kept the force acing in that direction but changed the line of action so it acted through the mesh point we would end up with a net moment in our system about the center of mass. This would cause an angular acceleration of our sphere. We clearly would not want our approximation of internal forces to result in this side effect.
So the idea is to move the force while keeping the moment about the center of mass the same and the magnitude of the force the same. Let:
ro be the vector from the center of mass to the point
where the force on the adjacent area acts.
Fo be the magnitude of the net pressure force acting on the
adjacent piece of area.
r1 be the vector from the center of mass to the mesh
point.
F1 be the average force of the section of area for wich the
mesh point is at the center..
Are Condition that the moments are the same gives
ro x Fo = r1 x F1
<=>
r1 x ( ro x Fo)
= r1 x ( r1 x F1) =
(r1 o F1) r1
- |r1|2 >F1 <=>
F1=|r1|-2 [r1 x (ro x Fo)]
+ |r1|-2(r1 o
F1) r1 <=>
F1=|r1|-2 [r1
x (ro x Fo)] + [(| F1|/|r1|)
cos ()] r1 <=>
F1=|r1|-2 [r1
x (ro x Fo)] + [(| Fo|/|r1|)
cos ()] r1 <=>
F1=|r1|-2 [(r1
o Fo) ro
- (r1 o ro) Fo]
+ [(| Fo|/|r1|) cos ()] r1
(9)
Equation (9) with the condition |Fo|=|F1| will not always be solvable in wich case we look for the value a value alpha that will make |Fo| as close to the magnitude of |F1| as possible. We want the solution that points away from the surface. Knowing that vector Fo points away from the surface and F1 acts close by, alpha will be close in value to the angle between ro and Fo. Thus, we look for solutions close to Fo that satisfy (9) and minimize | |Fo|-|F1| |.
In a sphere alpha will be exactly the same. For now we will use this approximation and it is assumed that the nature of the error will be such that any net moment induced by errors in our estimation of the infernal forces will cancel out.
(10)
P=pressure
z=altitude
Po=Initial Pressure
Zo=Initial Altitude
g=acceleration due to gravity
R=the gas constant (for air R=0.2870 kJ/Kg- K for helium R=2.07702
KJ/Kg-k) (2)
T=temperature
The initial pressure inside the balloon is determined by the ideal gas equation.
(11)
rho=the density of the gas
(The other variables are defined above)
In the code bellow the external pressure on bottom of the balloon is
given and does not change.
The code computes the pressure force on each mesh point.
%Input -[(x,y,z)_{i,j}]nxn the paramaterized
surface of the balloon
% -g the acceleration due to gravity
% -n the number of mesh points in each angular
direction
%Output -[(Pox,Poy,Poz)_{i,j}]n-1xn-1
% the net preasure force at [(x,y,z)_{i,j}]n-1xn-1
%Determine the area vectors
[Ax Ay Az]=surface_Area(x,y,z,n);
%Find the midpoint of the aproxmatly square
sections
[rmx,rmy,rmz]=surface_GridMidPoints(x,y,z,n);
%Compute The Volume
Vp=(1/3)*(Ax.*rmx+Ay.*rmy+Az.*rmz);
V=sum(sum(Vp));
%Get the properties of the gass inside and
outide the balloon
[Rb,Pb,Tb,Rair,Pair,Tair,m]=balloon_getIntialState_gass(V);
%Calculate the internal preasure ForcesP=mRT/V
minZ=min(min(rmz));
Pint=Pb*exp(-(rmz-minZ)*g/(Rb*Tb)); %The interal preasure force
on the aprox.
%square sections
%Calculate the external preasure Force
Pex=Pair*exp(-(rmz-minZ)*g/(Rair*Tair));
%Find the average force over the aproximately
square sections
Pmx=Ax.*(Pint-Pex);
Pmy=Ay.*(Pint-Pex);
Pmz=Az.*(Pint-Pex);
%Compute the angle beteen the preasure vecots
and the postin vector
%of the point of action
cosalpha=cosvec(rmx,rmy,rmz,Pmx,Pmy,Pmz);
%Find where an equivlanet force through the
lattic points would act
rox=x(1:length(x(:,1))-1,1:length(x(1,:))-1);
roy=y(1:length(y(:,1))-1,1:length(y(1,:))-1);
roz=z(1:length(z(:,1))-1,1:length(z(1,:))-1);
[Pox Poy Poz]=force1(Pmx,Pmy,Pmz,rmx,rmy,rmz,rox,roy,roz,cosalpha);
The code to find the spring constant is shown bellow:
Fkx=Fxi+Fxj;
Fky=Fyi+Fyj;
Fkz=Fzi+Fzj;
dxi=abs(x(1:n-1,1:n-1)-x(2:n,1:n-1));
dyi=abs(y(1:n-1,1:n-1)-y(2:n,1:n-1));
dzi=abs(z(1:n-1,1:n-1)-z(2:n,1:n-1));
dxj=abs(x(1:n-1,1:n-1)-x(1:n-1,2:n));
dyj=abs(y(1:n-1,1:n-1)-y(1:n-1,2:n));
dzj=abs(z(1:n-1,1:n-1)-z(1:n-1,2:n));
Di=sqrt(dxi.^2+dyi.^2+dzi.^2);
Dj=sqrt(dxj.^2+dyj.^2+dzj.^2);
The net force was obtained by simply adding up the force vectors. Since some of the matrices elements with different indices contained the same force it was necessary to make sure all forces that acted on these points were summed up. Finally the size of the matrix was expanded to nxn by assuming the forces at the end points were similar to their neighbors.
The initial velocity is defined. The mesh points are translated parallel to the velocity vector for at time delta t. The coordinates are then recalculated so the center of mass stays at the origin. The velocity of the lattice points is changed over a time delta t by the force. The forces are then recalculated for the balloon and the loop is repeated. The output is shown if figure 6.
%Get the surface (paramaters i,j coordanates
(x,y,z)
n=40; %Number of points in each angular dirction
Ro=1; %The Intial Radius
t=0; %Intial Time
delta_t=1e-5; %The incremental change in
time (s)
%Get the mesh points representing the balloons
surface.
[x y z]=balloon_newSurface(n,Ro);
g=9.8;
%Find the mass of the sections for wich thoese
preasure forces act
ml=balloon_getMassMeshPoints(x,y,z,n);
F_ml=-ml*g; %The force due to mass
ml(:,1)=ml(:,2);
ml=[ml ml(:,1)]; ml=[ml; ml(1,:)];
%Calculate The Intial Displacement Between
lattice points
[Doi,Doj]=Surface_Distance(x,y,z,n);
%Calculate The Spring Constants
[Ki,Kj]=balloon_getSpringConstants(Doi,Doj,n);
%Find The force acting on the mesh points
[Fx,Fy,Fz]=balloon_getForce(x,y,z,Doi,Doj,Ki,Kj,n,g,Ro,F_ml);
%Define Intial Velocity
Vx=zeros(n); Vy=zeros(n); Vz=zeros(n);
%Calculate New time
DELTA_T=8e-5; %The total change of time
t_final=t+DELTA_T;
while t<t_final
t=t+delta_t;
%Calculate New Postion
x=x+Vx.*delta_t; y=y+Vy.*delta_t; z=z+Vz.*delta_t;
%Make the orgin the center of mass
[x y z]=balloon_resetAxis(x,y,z,ml);
%Calculate New Velociety
Vx=Fx./ml*delta_t; Vy=Fy./ml*delta_t; Vz=Fz./ml*delta_t;
%Find the forces on the center of mass
[Fx,Fy,Fz]=balloon_getForce(x,y,z,Doi,Doj,Ki,Kj,n,g,Ro,F_ml);
%[x,y,z]=balloon_postionTouchUp(x,y,z,n);
mesh(x,y,z)
hold on
end
hold off
It is not clear why in the simulation it was unnecessary to introduce damping forces. It is likely in a lot of applications this will not be the case. Even for static analysis it would be desirable to able to watch how the balloon dynamically reaches equilibrium. This might make it easer to spot any fundamental modes of oscillation that may occur in a structure.
To make dynamic observations of the transition to equilibrium possible techniques could be introduced to select the mesh points to apply the laws of dynamics to at random. This would allow for a more gradual transition that could be observed through animation.
Further approximation techniques could be applied to allow the algorithm to iteratively converge on the correct solution. This could be done by further partitioning the surface after an initial approximation in made. To do this a new data structure will be required. This will be dealt with in future papers.
n=20; %Number of points in each
angular dirction
Ro=1; %The Intial Radius
theta=zeros(n); %Intialize
theta as an nxn array of zeros
phi=zeros(n); %intialize phi
r=Ro*ones(n); %intialize r
as an nxn array of Ro
n=n-1;
for i=1:n+1
for j=1:n+1
theta(i,j)=(2*pi*(i-1))/(n);
phi(i,j)=(pi*(j-1))/(n);
end
end
n=n+1;
%Transform to Cartesian
x=r.*sin(phi).*cos(theta);
y=r.*sin(phi).*sin(theta);
z=r.*cos(phi);
