Computer modelling

 

Choice of language

The choice of language is mainly governed on the processes required and the hardware the end product is to be run from. Currently, there are a number of languages (or mathematics packages) available, such as Mathematica and MatLab which would also be considered.

While the current trend is to use PCs, this is rather myopic as the trend for processing power is away from the PC range of machines (or at least, the Pentium style of processor) and onto much faster 64 and 128 bit processors, dominated by Alpha and ARM. For this reason, any software solely available for Windows is counted out.

Basing software on a particular operating system is not a good idea as it means when the operating system dies, so does the software. It also precludes the use on other machines without that particular operating system.

For this reason, from the outset FORTRAN 77 (FORmula TRANslation, a language developed in the early 60s for maths processing) was chosen. There is a F77 compiler available for most processors. For my work, I am using Acorn Desktop F77 from Intelligent Interfaces on a combination of StrongARM RiscPC and RiscStation machines.

While F77 is a very powerful maths based language, it has a number of drawbacks. Firstly is the age of the language. When originally developed, F77 (or F66 as it was) had all of the program stored on punch cards which a computer read in and executed. This meant that the only data out would be either saved as a file or teleprinted. The lack of graphical output requires the use of other software to interpret the results.

The language is also not very friendly to use (it is very strict as to where code must begin on a page and how constructs are formed). When the idiosyncrasies of the language are learned though, programming is fairly simple and rapid.

Other languages such as PASCAL and C were considered, but neither have the mathematical power of F77.

 

Software.

The software produced for this project falls into two distinct categories; calculation of the liquid bridges and calculation of the scattering, the listings for which are held in the appendix.

 

a. Liquid bridges

Calculation of the liquid bridge is performed by the use of elliptic integrals of the first and second kind as described by Carlson [9],[10] to produce the theoretical shape of a liquid bridge by the variation of vapour pressure on two perfect spheres (which to an approximation, is a good representation of SiO2 crystal).

Two separate programs are used to calculate the liquid bridges - one gives the full information (such as the fractional areas wet, volume of liquid surface of the sphere covered), the other produces a graphical representation for the shape of the bridge (in 1 dimension).

 

aa. Liquid bridges

This is primarily performed using FindEE9 and FindEE10. The full derivation for the calculations used are covered elsewhere..

The second program used is Britest.

Britest operates in a different way to FindEE in that the existence of the bridge is found and parameters returned to a routine to generate the shape of the bridge for plotting in a spreadsheet or graph plotting utility. As it does not require the to output the same amount of data, the program is smaller. The mathematics followed though is the same as for FindEE.

While the representation of the bridge in 1D is of use in showing the shape, a 3D representation gives a far better indication of shape. There are problems with the conversion from 1 to 3D.

Any two dimensional array can be simply converted to a three dimensional array by taking a third vector along the array (z). By creating this third vector from the first two vectors (x and y).

The third vector is found by the use of (52). NZ = number points along the z axis, S = inter-distance between two spheres, Rts = sphere radius, L = loop from 1 to NZ. For ease, NZ2 = NZ/2, φ = an angle at which the line is drawn along.

equation equation 52

52

The x and y vectors are now to be reprocessed to obtain the 3D image from the flat image. This is simply done by (53) and (54) where R(Z) is the radius of the sphere as a function of z.

equation 53

53

equation 54

54

φ is increased at a rate which will determine the resolution until φ=360 °.

Here in lies a major problem - plotting the data!

While to us, the data outputted would be simple to interpret (our minds work on a relational methodology whereby a point in space is related to another point and from that, we obtain a picture), a computer can only work on a "join the dots" system. Take for instance a simple 3D box.

3D Box

Here we can see what a simple 3D cube would seem to us. The 3 axis, x, y and z have been drawn in for ease of explanation. If the co-ordinates for this were to be shown on the graph we would have a simple x,y,z data sheet.

3D coords

The figures are set out as x, y, z co-ordinates with 0,0,0 being x = 0, y = 0 and z = 0 and 1,1,1 meaning x = 1, y = 1 and z = 1. Given this, you would assume that feeding these numbers into a spreadsheet or data plotting program would result in this 3D box being redrawn and sadly, you would be wrong. The real result is below.

3D real plot

Obviously, this means that any output produced will need either a dedicated program which can sort out this problem, or a simpler solution to adjust the way the data is being produced so that a matrix whereby a sort of "road map" is formed and a plotting program can follow.

This road map is a simple procedure to think of, but not to program!

Lined box

If the box is split into it's component parts, the formation of the new output can be clearly seen.

If the box is split into it's component parts, the formation of the new output can be clearly seen.

Take the red box first and construct from the bottom left around anticlockwise back to the starting position. Next, from the origin of the red square, we create a point at the end of the blue line 1. This is the staring point for the second box (green). Again, this is followed round.

Now comes the problem!. We have two boxes connected by blue line 1. How do we get lines 2 to 4?

Simple - we cheat!

It's not actually cheating, just repetition of the points. The two squares are connected by the blue line 1. If we go around the red square in the same way as before except to join up at the points were the blue lines 2, 3 and 4 should be, the full box is drawn. This does lead to a larger data file (and therefore needing more storage space and memory to process it in), but the representation is complete.

This simple approach can be performed on complex objects of revolution (such as the liquid bridge). Instead of a single box, we have a series of small boxes, but effectively, the same transformation takes place, but the size of data file will not be as large as would be assumed!

Consider the data array produced in 2D. If it were considered as a set of lines, it would look like this:

lines

The matrix would begin at the bottom left hand side and follow the path along. The red lines simply denote that this would be were there would be a join line.

The matrix would begin at the bottom left hand side and follow the path along. The red lines simply denote that this would be were there would be a join line.

When the line gets to the bottom right, the side connection lines need to be produced. This is simply done by connecting the bottom right to bottom left. The process restarts, but goes up to the a point on the y axis (x = 0) and joins to the point where x = n (far right)

This is carried on until the 2D matrix has all the points adjoined.

The 3D system is performed in the same way, but with the additional axis. The 2D block is produced, line draw back to 0,0,1 and the matrix reproduced in the z=1 area. The box system for data sets can now be produced. Overall, the data size will be larger, but not x9 larger (it's about x4 larger).

For a liquid bridge, we have a nodoid. This can be treated in exactly the same but using circles and ellipses instead of straight lines and boxes.

A simple bridge between two spheres in 2D gives a CSV file around 19K. The 3D profile is around 56Mb. It's apparent that the files produced are large!.

We already have a program which will generate the profile of the liquid bridge between two spheres of a given radius. The program is easily adjusted for a varied contact angle of the liquid to the sphere face. This is all well for a 1 dimensional system consisting of two spheres and a vapour / liquid. The aim of the project is to show a 3D representation of the FCC lattice and model both this and the SANs for the same system.

FCC lattice

This shows the arrangement from the front of an FCC lattice system showing the six nearest neighbours. It can be seen that there is the potential for 6 bridges to occur - one between the red central sphere and it's neighbour with the bottom edges of the bridge possibly being covered by the next sphere.

FCC flat

The view from the top and bottom will show and give rise to another 4 bridges - a total of 10 bridges possible.

The bridge contour can be generated for the 1D model and by manipulation, can be made to render the 3D model.

If a series of steps are taken up to the no bridge condition, a whole series of bridge contour generations can be constructed. This would be achieved by simply moving around the sphere by a small number of Radians and calculating the distance between the new position and the equivalent position on the opposing sphere. As the distances increase, the deformation due to the bridge elasticity will increase until no bridge is formed. This procedure can then be carried on for the remaining 9.5 bridges. The original will only be a half bridge contour map due to not revolving completely around the central sphere until it is met at the end.

Consideration in this way does lead to a problem. Until saturation is achieved, there will be an increasing degree of interaction between adjoining bridges. This can be best seen in the following graph :

Bridges

As the vap. pressure increase, the area of the bridges increase and also the amount of gap between the bridges decreases. This degree of interaction can also be found by variation of the vapour pressure from between 30% to around 98% (the maximum for water / silica before saturation)

The method of modelling has been greatly used in the video arts industry for a few years. It is known as contour mapping. The algorithm is quite simple.

We have a series of contour maps built up from the method outlined above. This will give the deformation and maps of the bridges between two spheres. If these are taken and at suitable points, lines drawn across them, a 3D contour map can be seen :

Contour maps

(obviously, it looks better than this!. These pictures were constructed to demonstrate a point rather than give the true representation).

The entire of these plots will form a deformed "hour glass" shape with the thinnest point at the centre.

As it stands, I have only considered 2 spheres from the 10. If all the spheres are of the same size and there is an even diffusion of vapour throughout the system then it can be assumed that all the bridges will be the same shape.

 

Drawbacks

The main drawbacks of performing these tasks are that huge amounts of data will be generated (each single profile as opposed to complete profile generates upto 2000 points. If a 0.5° angle was taken each time, then there is the potential to generate 720 ´ 2000 points = 1440000. If I take it that each point occupies a byte of memory (and taking 1K = 1024 bytes) then this will require 1406.25K (about 1.4Mb). There are 10 of these profiles meaning 14Mb of memory are needed for storing the 10 bridges. If I take 100 points on the central bridge and draw around the profiles, that make another 1000 points. There are 10 bridges making 10,000 points in total. Add these together and a singles 3D display of a single FCC lattice will require around 25Mb of memory (or equivalent amount of disc space). This is a large amount of space in anyones books!.

The second drawback is the speed of determination. Obviously, the full 25Mb of points would not be stored in memory - each profile would be dumped to disc and the next bridge determined. The problem occurs in the two sets of figure generation. While the bridge generations are very fast due to FFT techniques, the plotting of the contours will not be so simple. It should be possible to create a wire frame diagram of the liquid bridge, however, without FFT this may be a long haul. Even on a P233 or my SA233, I doubt that there will be much change out of 10 hours for single bridge relief map to be generated.

 

Advantages

The biggest advantage is that once a single profile has been generated and the position of the start and finishing points on the arc of the sphere for each bridge are determined, the full 3D representation can be determined.

The profile program is already working in F77, it will just required a DO loop adding in which would do something like this

        1 DO 2 LL=1, 1000
        2 (profile routine for a point on the arc)
        IF (BRIDGE.NE.1.) GOTO 10
        BRIDGEARC=BRIDGEARC+0.5D0
        CONTINUE
        10 (no bridge is made. store the data to disc)
        11 BRIDGEARC=BRIDGEARC+0.1D0
        CALL CHECKBRIDGE(BRIDGEARC)
        IF (BRIDGE.NE.0.) THEN
        GOTO 11
        ELSE
        IF (BRIDGE.EQ.1.AND.BRIDGEARC.GE.360) GOTO 20
        GOTO 1
        ENDIF

The transformation to a 3D representation the liquid bridge is simpler than would first thought of. By taking a point on the sphere where the liquid bridge comes into contact and then "wrapping" around the data which the program has generated, the following graph can be produced.

The liquid bridge

The translation from standard x, y and z data to the above picture is relatively simple and given in the formulae below

equation 55a

55a

equation 55b

55b

the value for the y axis depends on the absolute (positive value) of the z axis being less the the z axis curvature of the sphere. If it is less then (55c) applies otherwise (55d) applies

equation 55c

55c

equation 55d

55d

equation

While the above representation gives how the spheres and the bridge looks, it is simpler to consider just the bridge.