MENG21712 10 Credits
Modelling with computers TB1 Dr. A.G.W. Lawrie
Contents
1 Introduction 2
1.1 ComputingforMechanicalEngineers ……………. 2
1.2 EducationalAims……………………… 3
1.3 ComputingFacilities ……………………. 4
1.4 AttendanceandStudy …………………… 4
1.5 TB1Assessment………………………. 5
1.6 Rubricandintendedlearningoutcomes . . . . . . . . . . . . . . . 6
1.7 Referencematerial …………………….. 9
1.8 Documentorganisation…………………… 9
2 Solution techniques 10
2.1 Analyticalsolution …………………….. 11 2.2 Numericalsolution …………………….. 12 2.3 Algorithmandimplementation ………………. 14 2.4 Twodimensions………………………. 16 2.5 Underlyingprinciples……………………. 18
3 Problem Descriptions 20
3.1 A.Thermodynamics ……………………. 20 3.2 B.Fluiddynamics …………………….. 21 3.3 C.StructuralMechanics ………………….. 23
4 Optimising your code 25
4.1 Therationale ……………………….. 25 4.2 Efficiencyandscaling……………………. 26
1
5 Final Submission 28
5.1 Reportsubmissionguidance ………………… 28 5.2 Automatedmarking ……………………. 29 5.3 Codesubmissionspecification ……………….. 29 5.4 CallingScript ……………………….. 32 5.5 RopeBridgeWeightedMembrane exemplar . . . . . . . . . . . . 35 5.6 Finalchecklist……………………….. 37 5.7 Matlabfunctionlist……………………. 38
A Physics of heat conduction 41
B Stream function descriptions of incompressible flow 43
C Derivation of thin membrane mechanics 46
D Numerical convergence 49
1 Introduction
The aim of Modelling 2 is to introduce students to the art of using a range of modelling techniques to understand and find solutions to engineering problems. The part taken in Teaching Block 1 TB1 this year is a short project where the task is to model three closely related physical problems using Matlab. The followon component in Teaching Block 2 involves the designbuildtest of an amphibious vehicle. The following course notes pertain to TB1.
1.1 Computing for Mechanical Engineers
In ComputerBased Modelling MENG11511 in Year 1 you have been trained in the use of Matlab to encode a list of tasks into a procedural algorithm that can be executed on computer, using the software package Matlab. Matlab is a convenient Swiss army knife spanning everything from calculations you might perform on a pocket calculator to complex numerical problemsolving. The TB1 component of this course MENG21712 will teach the principles of solving real world physical problems whose structure is governed by partial differential equa tions. Our aim will be to develop the skills, tools and programming techniques to express engineering problems in a form that a computer can implement for us, then you will designbuildtest a method that with small adjustments can solve a broad range of such physical problems.
2
1.2 Educational Aims
Nowadays especially, enabling a computer to work for you is a core skill in any practicing Engineer, and more broadly across all numerate, technical and problem solving professions. Some students retain apprehensions regarding the use of com puters as a tool under their control to solve their problems. This is likely to be detrimental to their careers, whether they be in Engineering or elsewhere. This course aims to give you experience of using computers to produce credible Engi neering solutions to physical problems, and should allow you to grow comfortable with the tools, algorithms, implementation strategies and some level of technical background. Engaging with this course should help you overcome any remaining misgivings.
Any form of systematic procedure to solve a problem is called an algorithm, and it can be broken down into several small steps. Individually each step may be simple and easily understood, but as the length of the algorithm increases then the potential for unexpected behaviour quickly grows. Programming is a skill that develops through practice because it takes repeated exposure in order to restructure your own thinking so that it matches the rigorously systematic execution of a computer. Although devising an implementation of an algorithm is often tricky, understanding errors discovered when testing it is the real challenge in programming. With practice and experience, both things this course aims to provide, your errors will become less serious and less frequent, and your confidence will grow.
On this course we aim to provide a supportive environment for the development of programming skills. To manage your workload, a stagegate in week 5 will ensure that you have produced code that meets the passstandard of the course by the halfway point. During the final 5 weeks you will be encouraged to improve your codes performance, and pointed towards some techniques and algorithms that can help you achieve that. The final submission in week 10 to avoid collision with other deadlines will comprise your code and a report justifying your use of algorithms, implementation choices and presenting the results of your performance optimisation.
The university assesses coursework on a 21point scale that matches marks to descriptors of achievement, and firstclass 70 marks are awarded for going beyond what is explicitly taught in the course. To allow you to demonstrate this achievement, in the second half of the course support will still be available but you will be given guidance pointing you towards ways of improving your code rather than offering explicit teaching. As a further incentive, a prize is awarded for the bestperforming code!
3
1.3 Computing Facilities
The teaching computer resources for the Faculty of Engineering are located in the Merchant Venturers Building MVB on Woodland Rd. Before the start of this course all students should have a login and password for the Faculty computers. If not you will need to obtain one from IT services. Laboratory sessions will take place in both of the PC labs, rooms 1.07 and 1.08 in the MVB, and all students enrolled on the course will be given a 2hour slot in one of these rooms. In the laboratory there will be teaching assistants who can assist you wherever you have difficulty in completing your tasks. You must attend the laboratory slot you are given in your timetable because we have calculated your allocation to either group 1 or group 2 on the available capacity of each room, and distributed teaching assistants accordingly. Swapping session without permission from the unit director is not permitted. The University also provides nocost offcampus access to Matlab in a variety of ways: firstly the software itself can be downloaded onto your personal laptop using a licence key purchased by the University and available through IT services, www.bristol.ac.uksoftware and secondly the University provides Matlab as one of the installed packages on the Remote Desktop service www.bristol.ac.ukstudentdesktop. The various editions of Matlab all operate in slightly different ways and code developed on one version cannot be guaranteed to work on another. A common reason for this is that certain Matlab toolboxes may be installed in the version you have used to develop your code, but not on the version used to assess your code. This becomes a problem should you have called a function that is contained only in this toolbox. To avoid such problems, a list of admissible Matlab functions has been provided in 5.7. Before submitting your final work, it is your responsibility to ensure that your code only uses functions from the prescribed lists.
1.4 Attendance and Study
Attendance at laboratory sessions is the best means of obtaining assistance and feedback from teaching assistants, and is also compulsory. Attendance lists will be recorded for each laboratory session, and in total you will be expected to spend approximately 70 hours on this part of the course, including contact time and independent study. The introductory lectures will cover the mathematical background you need to understand the physical problem and the approach we will use to convert the problem into a form the computer can easily execute. The laboratory sessions will be performed in groups of 2 or 3 students that you will form in week 1. Each group will be required to hand in a group report.
4
1.5 TB1 Assessment
Although the work is performed in groups, submissions are made per group, and each package of submissions is assessed as a unit, marks are awarded to indi vidual students. In line with uuniversity policy on group work, confidential peer assessment is used, where appropriate, to redistribute marks, though it is anticipated that in most cases there will be no need to do so.
There is a stagegate that must be reached by your laboratory session on Monday 28th October 2019, and this ensures that your groups code meets the passstandard for this element of the course well in advance of the final dead line. The assessment will involve demonstrating your groups code to teaching assistant during the laboratory and they may also ask a few questions to test your understanding. For the purposes of this assessment, your code should produce simple graphical output, showing a final solution and a splot showing the path taken to converge to it. If satisfactory, the teaching assistant will sign each in dividuals stagegate assessment form and assign your group a unique identifier number. It is your responsiblity to return the completed form to the unit director.
For the final deadline, 5pm, Friday 6th December 2019, your submission must include a code a single Matlab .m file containing ASCII text only and a report a single .pdf file of no more than 10 pages produced using the LATEX doc ument preparation system, submitted electronically on Blackboard. Both will pass through a sequence of manual and electronic honestychecks to guard against plagiarism. We may also perform random spotchecks, and request constituent components of your report, eg. figures, .tex .bib files and so forth. Your report will be marked conventionally against the 21point scale, and your code will be tested exhaustively in an automated script to assess its performance on a range of testcases. For the purposes of this assessment your code must not produce graphical output or require user input.
Marks for each part of the course will be accumulated and an individual total compiled at the end of the year. The work in TB1 contributes in total 60 of the mark for this course, and this is subdivided into a 10 weighting towards the week 5 stagegate, 50 towards your code and 40 to the report.
There is no timed examination for this course. Note that no additional credit is awarded for the sophistication of graphical output: graphical user interfaces are elegant but timeconsuming to develop, and in this course the focus is on algorithmic efficiency rather than aesthetics of output. In recent years Matlab has tried to improve its underlying execution speed which is slow relative to other coding platforms and now offers a variety of tools to mitigate its disadvantages. Compilation to binary of any form is strictly forbidden since there is no way to verify authorship or quality of content, and the only admissible submission format is a single .m ASCII text file.
5
1.6 Rubric and intended learning outcomes
At the end of TB1 students should be able to:
represent realworld physical problem in mathematical notation using partial derivatives
transform partial differential equations into a sequence of instructions a computer can execute
devise systematic ways of testing instruction sequences code for errors
repackage code to conform to a specification
graphically present results obtained by using your code
use the LATEX document preparation system
summarise their work in a written report A firstclass student will also be able to:
optimise the implementation of an algorithm for efficiency
optimise the choice of algorithm to improve scaling performance
critically evaluate their choices, articulate their reasoning in some depth
show evidence of independent study of alternative algorithms
evaluate their codes performance over a suitably chosen range of testcases
6
Your code will primarily be evaluated by automated script that will perform a timing test on a range of test cases. This will test your code for
successful execution
returning correct results
correct error handling
scaling performance
In the event of unsuccessful execution, incorrect results, or incorrect error han dling, the automated script will continue on to the next case in the test suite. In the event of successful, correct results, a time threshold will decide if the scaling performance exhibited by your code warrants running a more complex test case. For consistency, codes from all groups will be tested on a common machine. All code will be evaluated against a benchmark implementation and graded accord ingly.
The report will be marked according to the Universitys 21point scale for coursework, and the relevant interpretation of the scale for this task is tabled beside the descriptor in figure 1.
There is no LATEXtemplate provided for this course, you should use your own choice of formatting. In general terms, your marker will be seeking to see:
clarity of structure
lucidity of prose: brevity without castration
sections and paragraphs used as intended
correct typesetting in LATEX
figures with correctly labelled axes
figures with captions
a relevant bibliography, cited in the body of text
Outline the physics of the problems, but focus on the algorithms developed in your code, key optimisations, and evaluate why and how well they work. One thing I strongly encourage is to interleave results and discussion together in a single section. Undiscussed results have little meaning, discussion detached from results involve a lot of tedious scrolling…
Bear in mind that your final mark for TB1 has three contributions, 10 from the week 5 stagegate, 50 from the final code submission and 40 from the report, so we are placing emphasis on quality and behaviour of your code over other elements of assessment.
7
1820: Report reaches publication standard in terms of content for a reputable academic journal, with little other than format adjustment to fit the style requirements of the journal.
1517: Clear evidence of originality, ingenuity and depth of understanding apparent from report content, with lucid written communication, and with the message intended of diagrams and figures unambiguous at a first reading.
1214: Clear evidence in the report of successful execution of the required tasks, expressed with written clarity and good support from diagrams and figures, but lacking in originality when seeking solutions to problems.
911: Evidence of some thought on the problems at hand from a thorough reading of the report, but the message is delivered without adequate clarity in text, diagrams or figures, and there is little or no evidence of originality or ingenuity.
68: Evidence of negligible effort in execution of the project andor preparation of the report, widespread lack of clarity in communicating a message and clear evidence of insufficient thought on the problems at hand.
05: Partial, invalid or nonsubmission of report, clear evidence of dysfunctional execution of the project andor preparation of the report.
Figure 1: Interpreting the University 21point marking scale for coursework.
8
1.7 Reference material
The main reference text which applies to this course is 1. For the programming aspects of the course, see also 2 and 3. In addition there are many good web based resources which include Matlab tutorials. These can be found by searching the web using a good search engine, as well as searching within the help function of Matlab. For the modelling aspects of the course, a further useful reference is 4. All necessary materials are provided in this set of course notes, but reading more broadly is one indicator of going beyond what is explicitly taught. In your report references in LATEX should be cited using citeHahn in the following manner,
beginthebibliography99
bibitemHahn B.D.Hahn, 2007. it Essential MATLAB
for engineers and scientists. 3rd Edition, Elsevier.
bibitemAUSCH M.Austin and D.Chancogne, 1999.
it Engineering programming: C, Matlab, Java, John Wiley.
bibitemHUNT B.R.Huntit et al, 2001. it A guide to
MATLAB : for beginners and experienced users, Cambridge University Press.
bibitemFAUS L.V.Fausett, 1999. it Applied numerical
analysis using MATLAB, Prentice Hall.
endthebibliography
resulting in formatted output as follows:
References
1 B.D. Hahn, 2007. Essential MATLAB for engineers and scientists. 3rd Edi tion, Elsevier.
2 M. Austin and D. Chancogne, 1999. Engineering programming: C, Matlab, Java, John Wiley.
3 B.R. Hunt et al, 2001. A guide to MATLAB : for beginners and experienced users, Cambridge University Press.
4 L.V. Fausett, 1999. Applied numerical analysis using MATLAB, Prentice Hall.
1.8 Document organisation
In this part of the course we will look at three physical problems that arise in different areas of engineering: Thermodynamics, Fluid dynamics and Structural
9
Mechanics, but are nonetheless related. The objective is to solve each of the given problems using the same solution technique and develop a single common code able to solve all three problems. This is possible because the differences between them are simply in the physical context, the underyling partial differential equation shares the same form and so they can all be solved using a common method.
The remainder of this document is organised as follows: The section on so lution techniques, 2 walks you through an approach to computermodelling of the problems posed. We demonstrate in 2.1 how to take an equation derived from physical principles and solve it on a simple onedimensional problem using analytical methods. Then in 2.2 we show how to discretise the equation, and in 2.3 illustrate how a numerical solution strategy can be implemented in Matlab. Thereafter, in 3 the specific remit for each problem is discussed in detail. In section 4 there are suggestions for ways to optimise your code for greater effi ciency, and ways to maintain efficiency as the problem resolution grows. Finally 5 contains important information about the final submission specification and format. In the appendices, we provide background information on how several different physical principles result in the same form of partial differential equa tion, and additional technical details on why numerical methods get less efficient as the problem size grows. The information in the appendices will inform your report writing and finalisation of the code towards the end of your project, but it is not necessary to read this material to begin work on your code.
The first step to make progress towards a solution is to copy the code for the rope bridge demonstrated in lectures into a text file that Matlab can run. A good approach in computer programming is to keep the programming running but slowly migrate it slowly from something simple that you already have working, to something more complex that does more of what you want. A sensible approach would be to modify the code in some simple ways so that instead of solving a 1D UDL the rope bridge problem it solves the analogous problem of a 2D UDL on a thin membrane. This fairly simple step opens the door to tackling all of the three projects listed the section on problem descriptions. The next step is then to modify this 2D code so that it addresses each specfic problem.
2 Solution techniques
This section introduces the common approach to solving all the problems posed. Most engineering problems can be decomposed into simple rules that govern be haviour. In many familiar cases these rules are based on Newtons laws coupled with a constitutive relation that comes from physics, for example, Hookes law to relate stress and strain, specific heat capacity to relate heat energy and tem perature, or the perfect gas law to relate pressures and densities in a gas.
10
Figure 2: Uniformly distributed load on a rope.
Here, well look at an inuitive problem to get started. A rope bridge is a centuriesold engineering solution to crossing a gorge, long before our local hero, Isambard Kingdom Brunel, designed the first suspension bridge in the modern style to span the Avon gorge at Clifton. Supported only at each end, and only able to resist axial stress and does not resist stress due to bending a rope bridge hangs under its own weight, as shown in diagram 2. The rule governing the vertical displacement, , of the system is as follows:
k2 mg, 1 x2
where m is the mass per unit length, g is gravity and x is the ordinate in the spanwise direction.
2.1 Analytical solution
Equation 1 is a differential equation we can solve directly by integrating twice. The first integration gives,
mgxC1, 2 x k
where C1 is a constant of integration. Integrating again we have,
x 1mgx2 C1xC2. 3
2k
We can set the constants C1 and C2 by applying boundary conditions. In the
trampoline problem, and as is common in many beam deformation problems, the boundaries have Dirichelet boundary conditions, the prescribed displacement 0. Wecansolveequation3atbothboundariesx0andxL:
x0,0:0 1mg02 0C1 C2 4 2k
0 C2 5 Similarly, and using the fact that C2 0,
xL,0:0 1mgL2C1L0 6 2k
C1 1mgL. 7 2k
11
So the deformation can be described as,
x 1mgx2 1mgLx0 8
2k 2k
mgxxL, 9
2k
which is consistent with the parabolic sag shown in diagram 2.
2.2 Numerical solution
As engineering problems become more realistic they become more complex and it becomes progressively more difficult to find analytical solutions. Instead we try to approach problems in a way that can handle more general cases, and we resort to numerical methods of the type we will study here. Analytical calculations still have a place in our understanding, particularly here in checking that our computer software is providing satisfactory answers. To use a computer to model an engineering rule like equation 1, we can discretise gradients using finite differences, essentially using the original definition on which Calculus is based,
fx fx h fx, 10 h
but instead of allowing h 0 we let h x, some small nonzero value. If we discretise the x axis into small segments, like links in a chain rather than a continuous rope, then the point where two links hook together we can in general call position i, the previous join is position i1 and the next join is position i1. Thus a gradient at a location halfway between i and i 1 is given by,
i1 i. 11 xi1 x
2
Often, as in this case, the rules of engineering problems involve second derivatives, the derivative of the derivative. Thus,
2
xi1 xi1
x i2 indices i,i1 backwards to i1,i we have,
i i1. xi1 x
2
x i2
2 2. 12 x 2 i x
. We already have the value for 1 , but we still need 1 . Shifting the
13
12
x , i1 i1
x
x , i1 i1
x , ii
x 2 2i1
x xi x xi x i1 1 i 1
Figure 3: A typical outofplane displacement in space.
as can be inferred from the typical representation of displacement against x shown in diagram 3. So now the second derivative is given by,
2 i1i x
ii1
Simplifying, we have,
x . 14
x2 i
2 i1 ii i1 i1 2i i1. 15
x2 i x2 x2 So the discretisation of our rule, equation 1, is as follows:
ki1 2i i1 mg. 16 x2
Written like this, it is expressed as a statement of fact: if we know all the correct values on the lefthandside then the expression must balance the righthandside. However in the absence of knowing the correct values, we are stuck. Instead we seek a procedure that takes some initial guess at what the values should be and update them until these guesses progressively become more acccurate. This process is known as convergence on the correct values. To implement such an approach we need a procedure that moves us from one guess to another, a number machine that takes in our first guess and returns our second guess, takes in this second guess and returns our third guess, and so on, subject to some parameters that come from the problem. The information we need is encoded in equation 16, but we need to turn this equation into an update procedure, a formula acting on the initial guess.
x
13
As a first step we rewrite our discrete equation 16, rearranging in terms of the value of at our general position i and obtain,
i1 2i i1 x2 mg 17 k
i i1i1 x2mg. 18 2 2k
This statement now explicitly defines a value of i that is dependent on its neigh
bours i1 and i1 and details of the problem, in fact we take the average of
surrounding values and account for the local contribution of the uniformly dis
tributed load as an additional constant term x2mg that increases with the load 2k
and decreases with stiffer ropes.
Equation 18 holds true when the elastic displacement responses i1, i and
i1 exactly balance the applied force, resulting in a static equilibrium. However it can also form the basis of a numerical method that starts marching from an initial guess and converges towards this exact equilibrium solution. If the marching index is n, then,
n1 ni1 ni1 x2 mg . 19 i 2 2k
This update formula only exactly corresponds to the original problem statement 16 when there is zero difference between n1 and n and so the update indices n and n1 are interchangeable and therefore can be safely ignored. We take this to be the definition of convergence. Prior to convergence, the current values of will not balance equation 16. By careful construction of the update formula and behind this lies an area of mathematics known as numerical analysis the de gree of imbalance of the equation is guaranteed always to reduce until it becomes asymptotically close to zero. In practice, we have finite patience and tend to trun cate our marching index when the degree of imbalance is within some acceptable userdefined and problemdependent limit.
More detail on how and why these methods work is provided below, but in the meantime one could code this update formula in Matlab.
2.3 Algorithm and implementation
The numerical method could be described in words as follows:
1. set physical constants, eg. gravity, spatial dimensions
2. set userspecified parameters for gridspacing, tolerances
3. create arrays to hold old and new values of the solution variable
14
4. set old array to starting guess
5. set values on the various boundaries according to the specified problem
6. calculate the updated value of at all points in the domain, using equation 19
7. estimate how far we are from convergence
8. copy updated new solution to old array in preparation for next iteration
9. if we are further from convergence than given tolerance, return to step 5, else continue
10. plot the function x vs. x
This stepbystep algorithm can be transformed into the following Matlab code.
set constants
nx100;
mass10;
gravity9.8;
stiff0.1;
length1;
phinewzerosnx;
phioldzerosnx;
dxlengthnx;
force0.5dxdxmassgravitystiff;
set up while loop
errc1;
erroerrc;
tol1e6;
it0;
while errcerrotol
enforce boundary displacement
phiold10;
phioldnx0;
calculate updated displacements
for ix2:nx1
phinewix0.5phioldix10.5phioldix1force;
end
calculate gradient of the error
errc0;
for ix1:nx
15
errcerrcabsphinewixphioldix;
end
errcerrcnx;
if it1
erroerrc;
end
swap old and new for next iteration after calculating error
for ix1:nx
phioldixphinewix;
end
itit1; end
plot0:dx:lengthdx,phinew;
Note in particular how the convergence tolerance is enforced. The tolerance tol is a userdefined constant, which in all the exercises is set to 1e6. The error between the current guess and the exact solution is never computed since we dont a priori know the exact answer. Instead we compute the rate of change of the error with respect to the iteration count. We do this by taking the absolute value in each grid cell of the difference between the current iteration and the previous one, then summing the contributions from each cell. To obtain a measure that is independent of the grid resolution nx, then we normalise this total error by the number of cells so that the error is expressed as a spatial average. We want to provide a relative measure of convergence that is common to all problems. The final converged value is measured with respect to the value at the end of the first iteration, given by the ratio errcerro. We are thus defining convergence to be obtained when the rate of change of error per iteration drops from its initial value to a value at least 6 orders of magnitude smaller.
2.4 Two dimensions
In more realistic engineering problems we often encounter more than just one spatial dimension, and constraints that cannot be resolved analytically. In the following, we consider problems in two dimensions where, like the rope bridge above, the physics requires discretisation of second derivatives. We generalise the idea of discretising like links in a chain, to discretising like a fishing net with square holes. We have to calculate gradients in both the x and y directions, so we introduce a second index j to identify the grid position in this new direction. A rectangular mesh with spacing x, y is shown in figure 4 In many cases the physics established for a one dimensional problem can be very simply translated
16
i,j1
i1,j i,j i1,j i,j1
Figure 4: Stencil on a rectangular mesh.
into two dimensions. Wherever derivatives 2 appear, these translate to x2
2 2
x2 y2 . 20
The equation is balanced by a term that represents the physics of the particular
problem, and in general we call it a forcing term fx,y. In the above one
dimensional example this forcing would be f x mg . The resulting equation is k
known as Poissons equation,
2 2 2
x2 y2 fx,y. 21
You are required to solve an equation of this form for each physical problem, but initially we will consider the special case of fx,y 0, which is known as Laplaces equation. We already know how the second derivative in x is discre tised:
2 i1 2i i1 , 22 x2 i x2
and the y derivative follows suit:
2 j1 2j j1 . 23
y2 j y2
When we put them together, we need both i and j indices simultaneously to clarify where exactly in space were retrieving values of . Our approximation to Laplaces equation becomes,
2220i1,j2i,ji1,j i,j12i,ji,j1. 24 x2 y2 x2 y2
17
Laplaces equation is a statement about the physics of a problem, not in itself a solution technique. Written in this form we require to already know values of x,y and we can then verify that the answer is zero, whereas we really want to establish a numerical method for finding those values, based on some initial guess. Following the suggested approach from the onedimensional example above, we rearrange the physical statement into a procedure that we can implement iteratively and believe will converge progressively from our initial guess to the correct answer. Our first step is to note that the spatially central value i,j appears more often than any other, it is said to be dominant in the solution, so we should perhaps treat it as special. We can rearrange our statement of physics so that it is described by a function g of neighbouring values:
i,j gi1,j,i1,j,i,j1,i,j1,x,y. 25 and then repeatedly compute new values of i,j until they converge. Multiplying
equation 24 by x2y2 to get a common denominator, we have,
0 y2 i1,j 2i,j i1,j x2 i,j1 2i,j i,j1 26
and teasing out the i,j terms, we have,
2y2i,j 2x2i,j y2i1,j i1,jx2i,j1 i,j1. 27
We really want the equation in terms of i,j alone, so we pull it out as a common factor and divide both sides by the other factor, ie.
i,j y2 i1,j i1,j x2 i,j1 i,j1. 28 2×2 2y2
Of course this becomes a much neater expression if y x, as on a regular square grid, and in this case would reduce to calculating the average of four
numbers:
i,j i1,j i1,j i,j1 i,j1 . 29 4
Coding this twodimensional update scheme equation in Matlab is a straight forward extension of the script provided above. Further details are provided in sections specific to each problem.
2.5 Underlying principles
So what does this equation mean? It says that the centre value i,j is the average of the surrounding values in x ad y. It turns out that this averaging process mirrors the behaviour of the physics of problems that contain second derivatives,
18
such as heat flow, fluid flow, and stress distribution. It implies that there can be no sudden jumps in the value of , since to raise one value egregiously with respect to its neighbours would mean that it was not the average of its surroundings. So how can we use it? Well, were aiming to develop an iterative scheme to walk towards the correct answer from some arbitrary guess. Until we get to the correct answer, or at least close enough, this equation probably wont balance. The amount by which it is unbalanced is called the residual the residue of error that prevents the physical law from being satisfied. This equation forms the basis of a technique called relaxation that works by taking an initial guess and repeatedly smoothing, averaging and relaxing values until the equation balances. So how does this work? We take equation 29 and although it solves a steadystate problem, we treat it in the same way wed treat a timeintegration in some sort of pseudotime as we march towards convergence. We say the RHS is evaluated at the old time level n, say, and the LHS is the new updated value at timelevel n 1, giving,
n1 ni1,j ni1,j ni,j1 ni,j1 . 30 i,j 4
Of course, when n1 n then weve converged and therefore the equation must be in balance and the steadystate correctly represented. Until convergence, it describes a marching scheme, which can be written in the following form for 1,
n1 n ni1,j ni1,j ni,j1 ni,j1 4ni,j . 31 i,j i,j 4
Rearranging, it is easier to see that for an appropriate timescale hidden in the change of timescale from t to is an implied constitutive relation of some sort the update scheme is effectively a discrete version of a time derivative.
n1 n n n n
i,j i,j i1,j i1,j i,j 1
4 2 2.
x2 y2
n 4n
i,j 1 c 32
33
So effectively our relaxation technique obtains a steadystate solution by evolving our arbitrary initial guess towards a steady state through a contrived sort of time. This observation is often used to appropriately modify the update formula to speed up convergence to a steady state. If we take larger steps through time, then we could reach a steady state faster. Of course as we converge, then rates of change of values reduce and by definition these time or pseudotime derivatives tend towards zero. When they reach zero, we satisfy the rule that was obtained from physics.
19
3 Problem Descriptions 3.1 A. Thermodynamics
This problem concerns the calculation of temperatures inside a plate of size 0.8m 0.4m as shown in figure 5. Temperature distributions T x, y satisfy Laplaces equation in stationary equilibrium,
2 2T 2T T x2 y2 0.
34
The plate is heated along parts of its top and left edges and cooled along parts of its lower and right edges. The faces and remainder of the edges are insulated. There is also a circular hole in the plate. The hole contains a highly conducting fluid, such that all points around the boundary of the hole have equal temperature. Your Matlab program will meet the following minimum requirements for this problem:
Uses a relaxation method to calculate the steadystate temperature distri bution in the plate specified below.
The user should be able to specify x and y resolution number of grid points, a temperature convergence criterion, and be able to change the boundary temperatures and the length of the boundary held hot or cold, and position of the circular hole.
For the stagegate assessment and preparing your report, the code should be able to provide clear graphical representation of a converged temperature distribution T x, y and a convergence plot of T vs. iteration count. This should be turned off for the code performance assessment.
The plates dimensions are 0.8m0.4m. By default the length along the vertical edge of the heated and cooled sections are 0.2m and 0.1m along the horizontal, and maintained respectively at 100K and 0K. By default the circular hole has radius 0.05m and is positioned at 0.3m to the right and 0.3m above the lower left corner of the plate.
In locations where boundaries have a prescribed temperature, the value at points on the outer edge is held fixed. Boundaries that are described as insulating indicate that there is no heat flux through that part of the surface. Where there is no heat flux, there can be no temperature gradient driving one, so the edge value and the nextinnermost value must be held to be identical though with no constraint on the specific value. On the edge of the circular hole, the values of temperature must all hold a common value, though again there is no constraint on the specific value, simply that they must hold the same value.
20
Figure 5: Boundary conditions for heated plate. Figure not to scale.
3.2 B. Fluid dynamics
This problem concerns the calculation of steadystate incompressible flow patterns around obstacles in a twodimensional channel. Specifically, the purpose is to estimate the directions of streamlines and local velocities in a 4m long section of a 2m wide channel partially obstructed by several cylinders circles, as depicted below. In steady state, the velocities U and V can be described by a single variable called a stream function x, y according to the following relationships,
x
y
2 2 2
x2 y2 0. 36
V
U.
Under appropriate conditions detailed in appendix B, this stream function sat
isfies Laplaces equation,
In addition to the requirements for problem A, your Matlab program must also satisfy the following minimum requirements for problem B:
Uses a relaxation method to calculate the steadystate velocity field in the 2D channel specified below.
35
21
!
Figure 6: Typical arrangement of Streamlines around cylinders. Figure not to scale.
The user should be able to specify x and y resolution number of grid points, a convergence criterion for stream function, and be able to change the inflow velocity, and the positionradius of the cylinders.
For the stagegate assessment and preparing your report, the code should be able to provide clear graphical representation of a converged velocity
Ux,y
field V x, y , plotted as streamlines contours of stream function, and
a convergence plot of vs. iteration count. This should be turned off for the code performance assessment.
The quantities and dimensions required to construct the geometry are tabu lated for your convenience:
Quantity
Channel length Channel width Upstream flow speed Upstream flow speed
Symbol Dimension
x 4m
y 2m
U 1ms1
V 0ms1
Cylinder locations and radii for the default problem are:
22
Cyl xc m yc m radius
1 0.5
2 1.0
3 1.5
4 1.9
5 2.5
6 2.3
0.7 0.2 1.4 0.3 0.5 0.3 1.1 0.2 0.4 0.2 1.6 0.2
On bounding walls and on cylinders in the interior, there is no penetration of flow through these boundaries. In such cases, the stream function which is an indirect measure of mass flux per unit distance takes a constant value. There is no direct way of specifying what particular value that should be, just that all values on the same boundary should have a common value. To provide a baseline value, set 0 in the bottom left corner. Thus all of the bottom wall will have 0. To be clear, the top wall and the cylinders will all have different, but individually constant values. For inflow we use the formula,
y Uy 37
ie. a linear variation of with y whose gradient is the prescribed upstream flow
speed. It follows that the stream function value on the wall at y 2m has a
constant value 2U. At the outflow the value of must vary smoothly from
0 to 2U, however disturbances from the cylinders may still affect the
distribution of velocity in the wake. To account for this, we do not enforce a
specific velocity profile, but simply that there is no V component to the outflow
velocity. From equation 35 this implies 0 and we implement this by holding x
the edge value and the nextinnermost value identical though with no constraint on the specific value.
3.3 C. Structural Mechanics
This problem concerns the displacement of a trampoline when subject to the weight of a gymnast standing offcentre as depicted in figure 7. Stretched mem branes can only resist inplane stress, and do not support bending moment. There is, unfortunately, a rip in the trampoline skin. Stress cannot be borne across the rip; stresses can only be supported in the direction parallel to the rip.
In steady state, the displacement u satisfies Poissons equation, 2 2u 2u
u x2 y2 fx,y, 38 where the forcing is provided by the weight of the gymnast distributed over two
23
Figure 7: Skin of a trampoline with numerical grid superimposed, showing typical gymnast contact patches, a typical rip in the membrane and the zerodisplacement boundary around the edge. Figure not to scale.
small circles around their feet. In those two circles, the forcing is given by fx,y mg 39
Ak
where m is the mass of the gymnast, g is gravity, A is the total area of contact of the gymnasts feet, and k is the stiffness of the elastic material. In addition to the requirements for problems A and B, your Matlab program must also satisfy the following minimum requirements for problem C:
Uses a relaxation method to calculate the steadystate displacement of the trampoline skin specified below.
The user should be able to specify x and y resolution number of grid points, a convergence criterion for displacement, and be able to change the position of the rip, and the positionweight of the gymnasts feet, and the extensional stiffness of the material.
For the stagegate assessment and preparing your report the code should be able to provide clear graphical representation of a converged displacement ux,y and a convergence plot of u vs. iteration count. This should be turned off for the code performance assessment.
The trampoline has a length of 4m and a width of 2m and is fixed around the edge. Taking the bottomleft corner as an origin, the default location of the
24
rip runs from 1.5m, 0.5m to 1.5m, 0.8m. The default weight of the gymnast is 56kg, the contact area of their feet is approximated by two circles of radius 20cm. By default the gymnasts feet lie at 3m, 0.75m and 3m, 1.25m. The stiffness value you shoudl choose to scale your computation is largely determined by the precise details of how you define it whether by length, area, normalised strain etc. for our purposes here it is not important how it is defined, so choose an appropriate value for the materials extensional stiffness such that the maximum displacement is around 45 50cm, typical of the deflections found on a real trampoline.
All of the edges of the trampoline are fixed, so they have displacement u 0 everywhere on the edge. The rip is not able to support stress in the direction across the rip ie. in the xdirection, but it can support stress in the direction parallel to the rip just as well as everywhere else. Since outofplane displacement gradients are an elastic response to inplane stress, then there can be no gradient of displacement in the xdirection anywhere adjacent to the rip on either side. Thus values on the very edge of the rip must be equal to those on the nextinner most row of grid points, and the equivalent must be true on the opposite side of the rip. Note that the displacements at both ends of the rip must be equal, so the average difference in stress between the end points must be equal. In this case we can make the convenient modelling approximation that the rip lies exactly coincident with a mesh boundary, but the vertical displacements on either side of the rip may take different values.
4 Optimising your code 4.1 The rationale
Now that you have passed the Week 5 stagegate you are expected to have pro grams that solve each of the three problems, thermodynamics, fluid dynamics, and structural mechanics. You will have noticed that the time taken to converge to a solution is dependent both on the problem youre solving and how you write the code that solves it. Optimisation is the process of carefully rewriting a work ing code to make it faster. Good programming practice dictates that you first get a program working before you attempt to optimise it: only once you have full oversight of what the code needs to do in order to work should you consider taking shortcuts that may or may not make the code more efficient. The Week 5 stagegate has now given you this oversight and only now should you proceed to optimise your code. Code performance on the sorts of physical problems you study in this unit is not simply a matter of timing code on a particular case: the aim is to have a general purpose tool that solves a range of related physical problems to a userspecified accuracy in terms of convergence tolerance and grid resolution.
25
4.2 Efficiency and scaling
Timing information from one particular case can be misleading because there are, broadly speaking, two components to the evaluation of a Matlab command: parsing and execution. A command is parsed from a string of characters that you ask Matlab to evaluate, and only then can it be executed as an operation acting on data held in memory. The parsing step has a considerable overhead cost, which is dependent on the length and complexity of the character string that it receives. The execution cost is the bit that gets the job done, and as the problem size increases, the cost of parsing doesnt change but the relative fraction of total time spent doing useful work on data held in memory increases. Thus at very low grid resolution the cost of parsing dominates, but the overhead decreases proportionally as grid resolution increases.
Matlab have tried very hard to reduce the overhead costs, and the vectorised notation : for scanning over an array helps with this, since operations that would naturally be expressed over several lines as a loop operation can be expressed in a single line:
phinew1,:0.0;
has the same effect as
for iy1:ny
phinew1,iy0.0;
end
but requires fewer parsing operations to interpret. Internally the operation is executed by a for loop in a compiled language usually C or Fortran in which the parsing step is performed in advance rather than onthefly and the for operation is encoded as lowlevel binary instructions that execute much more rapidly. Execution speed in Matlab has long been a drawback of the package, especially with computationally intensive calculations.
The use of for and while shoudl be minimised in Matlab scripts. Careful coding of array operations can yield equivalent outcomes to explicitly coded loops, but are more efficiently executed. Consider the following code:
x 0:dx:gridlength;
y 0:dy:gridheight;
nx lengthx; ny lengthy;
xarr, yarr meshgridx, y;
phi zerosny, nx;
This generates arrays x arr and y arr having appropriate entries note that a sensible array ordering and standard matrix rowcolumn ordering are not always synonymous. For instance, this allows statements such as:
26
dirichelet1 findyarr gridheight xarr gridlength bcportion;
dirichelet2 findyarr gridheight bcportion xarr 0;
diricheletboundary uniondirichelet1,dirichelet2;
where Dirichelet is the type of boundary where the value is fixed. These state ments identify the grid points that can then be set to the desired value:
phidiricheletboundary boundaryvalue;
Similar constructions can be used to identify other boundaries. Naturally there are other methods to achieve the same results, but they are generally less easy to interpret than the above. As an example, in the special case x y, the update step can be performed extremely efficiently:
innerx 2:nx1;
innery 2:ny1;
xdiff phiinnery, innerx1 phiinnery, innerx1;
ydiff phiinnery1, innerx phiinnery1, innerx;
oldinner tempinnery, innerx;
phiinnery, innerx xdiff ydiff4;
though it should be noted that in the general case x y and your code should have the capability to handle the general case
There are other ways of improving the efficiency of your program; the gap between your starting guess and the final answer plays a big role in the time taken, so how would you improve your starting guess? You could start by solving a simpler relevant problem, and using the answer to this simpler problem as a starting guess. If the target resolution is, say, an 80 40 grid, you could instead start with a 10 5 grid, a much larger x, y, and a much shorter time to converge on this approximate problem. Then interpolate this coarsely resolved solution as a starting guess of a finer mesh, maybe 20 10, solve this, then repeat on 40 20, then finally on an 80 40. This method is a form of multigrid, sonamed for obvious reasons, and it one of the fastest and most flexible ways of solving Poissonequation problems. The copying and interpolation process from a coarse grid to a finer one takes a little care to get the indexing correct, because only every second value in each direction will have a directly supplied value. You will need to perform a simple linear interpolation to fill in the gaps. One hint for designing your multigrid: refinement by factors of 2 require interpolations that are straightforward to implement, so it makes sense to choose the number of gridpoints in x and y to be a power of two. You may well be asked to perform the calculation on a grid with a nonpowerof2 on either side, so you should consider carefully your strategy for such cases.
27
In all cases you must check that your modifications actually make your code more efficient… From Matlab 2015a onwards there is a new button Run and time on the top bar, which will give you performance statistics about your script. In general multigrid will actually run less efficiently on very small problems, but will get progressively more efficient as the number of x and y grid points increases. This is a favourable scaling property. The basic iterative method gets very slow as the problem size increases, but this may be mitigated using a wellcrafted multigrid method. You should plot the performance of your code at a variety of resolutions. It is often instructive to consider the time taken to converge divided by the total number of grid points as a measure of efficiency. If the resulting graph were horizontal, then the time taken scales linearly with the problem size, and this is the ideal case for an efficient code.
The second way in which examining just one particular case can be misleading is that time to convergence doesnt scale linearly with grid resolution. You might expect that the cost per grid point is a constant, independent of the size of the problem, but in fact the cost per grid point increases. The reason for this is that physical situations that give rise to Poisson or Laplace equations have one key thing in common. A change of property at any location in space affects every other point, so a loose analogy might be the handshakes problem 10 people in a room, 45 handshakes needed every point is connected to every other and so it takes more iterations before everything equilibrates and the current solution balances the original equation to within an acceptable tolerance. Some of the linear algebra behind this is detailed in the appendix, for those interested.
5 Final Submission
5.1 Report submission guidance
Only one person from your group should submit a copy of the report. The only acceptable format is .pdf do not try to upload .zip files or any other format. The filename of the uploaded .pdf should be group.pdf where is your group number expressed as three digits, and lowercase g is used. For example group 7 would submit a file group007.pdf. No other filenames will be accepted.
The only other stipulation on format is that the .pdf contains less than 10 pages, and note that there is no LATEXtemplate: you have a free format within those 10 pages, but do take note of the marking rubric given in 1.6, which suggests how you might organise your submission.
28
5.2 Automated marking
The need to test how your code performs over a range of resolutions is imperative, so we need to invoke your program multiple times for each problem at different res olutions and convergence tolerances. We will then plot the performance statistics we need to evaluate correctness, error handling and scaling performance of your code and assign a mark accordingly. It is unnecessary to do this datagathering by hand, provided a common format is used to invoke your program, and the solution is returned in a common format. The specification in 5.3 details exactly the arguments we will supply to invoke your code, and the solution we expect in return. You will need to build a small wrapper around your code to accept our arguments and return the data in the correct format. To ensure you can prepare your code accordingly, a simplified assessment script is provided in 5.4, and af ter the week 5 stage gate a representative exemplar submission solving the 2D rope bridge problem will be issued, showing how to integrate your code with the automated script.
5.3 Code submission specification
The code submission is 50 of the computational part of Modelling 2. You should submit a single .m file containing all your code that you wish to be assessed. The filename of the uploaded .m should be group.m where is your group number expressed as three digits, and lowercase g is used. For example group 7 would submit a file group007.m. Inside this file the name of the first function must correspond directly to the filename. It must also accept input and output arguments in the following format, illustrated for group 7:
function phiend,converge,st1,st2,st3,error…
group007phistart,tol,prob,probargs;
put your code in here
end
The input arguments are specified as follows:
phistart is the array on which you are expected to return a solution. You will need both the values and its its dimensions number of elements in each row and column. Note carefully that the grid will be interpreted as vertex centred, so elements on the last row or column of the array represent values that are on the boundary.
tol represents the convergence tolerance, and this is defined to behave in the same way as the template provided: it measures the gradient of error at the beginning comparing the original phistart with the values after one
29
iteration with respect to the gradient of error at the current iteration com paring current with immediately previous, and your iterative loop should finish when the ratio of starting to current gradient is less than the specified value.
prob is a single character specifier for the problem you will solve and set ting it to a should solve the thermodynamics problem, b should solve the fluid dynamics problem and c should solve the structural mechanics problem. In the case shown below, d is used for the rope bridge exemplar.
probargs is a problemspecific parameter array, whose format is detailed next.
The input argument probargs will in each case just be an array of numbers cod ifying the setup of the problem. For the thermodynamics case a, the argument variables will determine only the location and size of the circle, and are specified in the following order:
probargs1x; x coordinate of centre
probargs2y; y coordinate of centre
probargs3r; radius of circle
For the fluid dynamics case b, the argument variables will determine only the location and sizes of the circles, and are specified in the following order:
probargs1x1; x coordinate of centre of circle 1
probargs2y1; y coordinate of centre of circle 1
probargs3r1; radius of circle 1
probargs4x2; x coordinate of centre of circle 2
probargs5y2; y coordinate of centre of circle 2
probargs6r2; radius of circle 2
…
and there may be more than 6 circles specified as input so you will have to check the size of the incoming array to be sure how many there will be. For the structural mechanics case c, the argument variables will determine the position of the rip and the position of the circles but not their loading which should be preset along with an appropriate stiffness for the membrane to achieve the specified maximum displacement.
probargs1x1; x coordinate of centre of circle 1
probargs2y1; y coordinate of centre of circle 1
probargs3r1; radius of circle 1
30
probargs4x2; x coordinate of centre of circle 2
probargs5y2; y coordinate of centre of circle 2
probargs6r2; radius of circle 2
probargs7xrip; x coordinate of the rip
probargs8ybot; y coordinate of the bottom of the rip
probargs9ytop; y coordinate of the top of the rip
Note that the rip will always be oriented as per the description surrounding figure 7. All dimensions and coordinates will be in the units of the original problem specification.
All information not provided by these arguments must be hard coded into the function.
The output arguments are specified as follows:
phiend is the array which contains the solution. This must have exactly as many elements in each row and column as phistart. Note carefully that the grid will be interpreted as vertex centred, so elements on the last row or column of the array represent values that are on the boundary. You should comment in your report about any consequences for accuracy that this vertexcentred representation introduces. Be careful to ensure that numelsizephiend is identical to numelsizephistart. The safest approach is to set phiendphistart and copy values into an unchanged array.
converge is a two column data array supplying data for a 2D plot of con vergence gradient ratio vs iteration count. The orientation of the array that the calling script expects is given by convergezerosnrows,2;.
st1, st2 and st3 are character strings for the first part of the email addresses for each member of the group, and these should be set as follows:
st1jb007;
st2cm16039;
st3null; only 2 in this group
Failing to provide these in the format specified here will mean your code will not contribute to your overall mark.
error should return 0 if your code has reached convergence and you are confident that the solution is a good response to the input parameters. Instead return 1 if your code recognises that there is an unsuitable input specification, or 2 if a failure to converge has occurred during execution.
31
5.4 Calling Script
modelling 2 201920
simplifed script for automated assessement
clear all
close all
ngroups61;
timeout30;
nx100;
ny100;
default settings for problems A. B. and C.
thermodynamics 0.8m x 0.4m domain
probazeros3,1;
proba1,10.3; x coordinate of centre
proba2,10.3; y coordinate of centre
proba3,10.05; radius of circle
fluid dynamics 4m x 2m domain
probbzeros18,1;
probb1,10.5; circle 1
probb2,10.7;
probb3,10.2;
probb4,11.0; circle 2
probb5,11.4;
probb6,10.3;
probb7,11.5; circle 3
probb8,10.5;
probb9,10.3;
probb10,11.9; circle 4
probb11,11.1;
probb12,10.2;
probb13,12.5; circle 5
probb14,10.4;
probb15,10.2;
probb16,12.3; circle 6
probb17,11.6;
probb18,10.2;
structural mechanics 4m x 2m domain
32
probczeros9,1;
probc1,13.0; circle 1
probc2,10.9;
probc3,10.2;
probc4,13.0; circle 2
probc5,11.1;
probc6,10.2;
probc7,11.5; rip xcoordinate
probc8,10.5; rip ystart
probc9,10.8; rip yfinish
rope bridge
probdzeros5,1;
probd110; mass
probd29.8; gravity
probd30.1; stiff
probd41; length
set up space for output statistics
durationAzerosngroups,1;
durationBzerosngroups,1;
durationCzerosngroups,1;
durationDzerosngroups,1;
userszerosngroups,3,10;
crashzerosngroups,3;
figure1;
for ig7:7 1:ngroups
groupsprintfgroup03d,ig;
dispgroup;
test sequence
probnumd;
probcurprobd;
tol1e6;
phistartzerosnx,ny;
try
tstartcputime;
the line which calls your code is below:
33
phiend,converge,st1,st2,st3,error…
fevalgroup,phistart,tol,probnum,probcur;
durationDig,1cputimetstart;
usersig,1,1:numelst1st1:;
usersig,2,1:numelst2st2:;
usersig,3,1:numelst3st3:;
if error0 expecting no error on this testcase
errsumsumabsphiendphistart;
if err1e4
disptest 1 failed to return a valid solution;
crashig,11; flag unrecognised internal error
else
verify solution convergence by further iteration…
display outputs
dispdurationDig,1;
subplot3,2,1;
plotconverge:,1,converge:,2;
subplot3,2,2;
surfphiend;
consider recursively doubling resolution
if durationAig,1timeout
close all
clear all
phistartzerosnx2,ny2;
try…
catch…
end
end end
else
disptest 1 failed to recognise incorrect input;
crashig,12;
end catch
disptest 1 crashed internally;
crashig,13; flag unrecoverable internal error
end
34
end
5.5 RopeBridgeWeightedMembrane exemplar
function phiend,converge,st1,st2,st3,error…
group007phistart,tol,prob,probargs
prepare userids
st1jb007;
st2cm16039;
st3null; only 2 in this group
if probd
obtain resolution and initial guess from phistart
szsizephistart;
nxsz1;
nysz2;
phinewphistart;
phioldphistart;
phiendphistart; ensure output array has right size
create some space for convergence plot
convergezeros30000,2;
set constants
massprobargs1;
gravityprobargs2;
stiffprobargs3;
lengthprobargs4;
dxlengthnx;
force0.5dxdxmassgravitystiff;
set up while loop
errc1;
erroerrc;
it1;
while errcerrotol
enforce boundary displacement
phiold1,:0;
phiold:,10;
35
phioldnx,:0;
phiold:,ny0;
calculate updated displacements
for iy2:ny1
for ix2:nx1
phinewix,iy0.25phioldix1,iy0.25phioldix1,iy…
0.25phioldix,iy10.25phioldix,iy1force;
end end
calculate error gradient
errc0;
for iy1:ny
for ix1:nx
errcerrcabsphinewix,iyphioldix,iy;
end end
errcerrcnxny;
if it1
erroerrc;
end
swap old and new for next iteration after calculating error
for iy1:ny
for ix1:nx
phioldix,iyphinewix,iy;
end end
convergeit,1it;
convergeit,2errc;
itit1;
end
phiend1:nx,1:nyphinew1:nx,1:ny; copy into output array
error0; rash assumption: dont do this at home!
end end
36
5.6 Final checklist
Never trust user input. You have no control over what the automated script will ask your code to do, you cant be certain that youve interpreted the units in the same way as the user and it is your responsibility to ensure that the code returns gracefully. Dont assume that because the automated script will have lots of code to mark that it will not attempt to pass inconsistent or corrupted data. It will. Matlab has two keywords, try and catch that enable us to test your code in a safe playpen and if it crashes the marking script will simply continue to the next invocation attempt.
Ensureyourusernamesareprovidedwithinthefileandreturnedasargument in the correct manner, otherwise your marks will go to someone else, or to nobody!
Ensure your code returns a nonzero error if the combination of input arguments dont make sense, or your code cant solve the problem posed by these input arguments.
Ensure your code measures the initial rate of change of error from the phis tart array exactly as given to you as an input argument: the input may not be a zero initial condition. Measure the ratio of your current rate of change of error to this initial rate of change of error to decide whether you have met the convergence tolerance before returning phiend as an argument. This calculation must be done at the resolution given by the phistart input argument.
Ensure that the array size of your return argument phiend exactly matches the array size of the input argument phistart, and that the requested convergence tolerance tol as defined above is met at this grid resolution, not just on any other internal grids you may use.
Ensure that the physical problem the boundary conditions, geometric di mensions, positions of circles etc. you solve is independent excepting dis cretisation errors of the dimensions of the phistart array; there is no guar antee that the array dimensions will match the aspect ratio of the physical dimensions.
Ensure that your code does not require user intervention GUIs, waiting for input values etc. because no such intervention will be forthcoming, and the waiting time will be treated as execution time and count against you.
Ensure that your final submission does not produce any plots. 37
Ensure that your code only uses functions in the prescribed lists given in 5.7.
5.7 Matlab function list
The following lists show several available functions and arithemticalgebraic oper ators that are commonly used in Matlab. Your program, when finally submitted, is likely to make use of commands from the first two lists: Programming, logic and arithmetic and Common Functions. There is unlikely to be need to use commands from the subsequent lists in your final submission, however the plotting commands will be especially useful for debugging and prepring plots to present in your report. To obtain a good mark for your code, it is essential that you remove any commands that will require user input, pauses execution for any reason, or displays information to screen, because this will slow your code down and possibly halt it altogether, resulting in a performance penalty that you can easily avoid.
Programming, logic and arithmetic operators
:
. .
for
if, else, elseif, end while
Perform operation on an entire array rowcolumn Arithmetic on scalars
Algebraic operations on matrices
Gaussian elimination
Elementwise operations on arrays
Logical operator for and
Logical operator for or
Relational operator for less than
Relational operator for greater than Relational operator for greater than or equal Relational operator for less than or equal Relational operator for equal
Relational operator for not Open a for loop Conditional statements Open a while loop
38
Common Functions
abs Absolute value
ceil Round up to nearest integer values
cos Return cosine of input variable
max Maximum value in a rowcolumn of an array mean Mean value of a given array
min Minimum value in a rowcolumn of an array round Round valuesto nearest integer
sin Return sine of input variable
sign Returns sign of variable
size Return size of an array
numel Number of elements in an array
pi 3.141… to machine precision
tan Return tangent of input variable
zeros Array of zeros
ones Array of ones
Plotting commands
bar Plot a bar chart
compass Compass plot
get List attributes of plots etc.
grid Show grid
gtext Put text on the plot using mouse
hist Plot a histogram
hold Maintain current plot while others are plotted mesh Mesh surface plot
plot 2d plot function
plot3 3d plot
set Change attributes of plots etc.
stairs Stairstep plot
stem Plot a stem plot
surf 3d coloured surface plot
text Put text on the plot
title Plot title
xlabel Label for the x axis
ylabel Label for the y axis
39
Matrix algebra, complex number and statistical commands
inv Matrix inverse
pinv Matrix pseudoinverse
tril Lower triangular part of a matrix triu Upper triangular part of a matrix
det Determinant of a matrix
eye Identity matrix
mldivide Left matrix divide
conj Complex conjugate
real Real part of complex number
imag Imaginary part of complex number angle Phase angle or complex argument rand Generate a random number
randn Normally distributed random numbers cov Covariance
std Standard deviation of a given array
40
Miscellaneous commands
cd
clear
clock
date
diary
dir, ls
edit
eps
format
help
helpwin
load
lookfor
path
quit
save
who
whos
input
pause
disp
clc
fprintf
num2str
meshgrid
Change directory
Clear all variables
Current date and time as date vector Display current date
Start recording a diary
List all files in the directory
Start the editor
Floating point relative accuracy
Set output format for displaying numbers Display Matlab help related to command Display the help window
Load file into Matlab workspace
Look for commands or .m files
Getset search path
Exit from Matlab
Save workspace variables to disk
List current variables
List current variables showing type and size Prompt user for input
Wait for user response of length of time Display input variable
Clear command window
Write data output to a file
Convert a number to a string
Create an indexed array
A Physics of heat conduction
Fouriers Law of thermal conduction says that the local heat flux density, q, is equal to the product of materials thermal conductivity, k, and the negative local temperature gradient, T. The heat flux density is the amount of energy that flows through a unit area per unit time. The reason for the negative sign is that the heat flows away from hot things towards cold things, downhill and therefore against the temperature gradient as its vector direction is defined. Thus the heat flux is given by,
or in one dimension, simply,
q kT, 40 q kT . 41
x
41
To describe how temperature or any other quantity diffuses in time through a medium, we invoke Ficks 2nd Law of diffusion, which in one dimension says that
T q. 42 t x
To see why this makes physical sense, we integrate the LHS in space over an arbi trary length 2L. L T dx cant be solved immediately, but we recognise that in
L t
most physical situations we can swap the order of integration and differentiation,
and equivalently write L T dx, which is now much easier to interpret as the t L
time rate of change of the total amount of T in the length 2L. Doing the same integration on the RHS, we can say that L q dx qL qL qL is the
L x L
difference in the amount of q at both extremes of length 2L at this timeinstant.
Remember, q is a flux of heat, and so heat could be flowing very quickly through the length 2L without the total amount contained in the length 2L changing at all. The temperature T will only increase in time if there is more heat flowing in through, say, the L boundary than, say, is flowing out through the L boundary of our arbitrary length.
Substituting the heat flux definition 41 into equation 42, we realise we have a mixed partial differential equation in one dimension,
T T
tx kx . 43
Often the thermal conductivity k is not a spatial variable, so it can be pulled outside the operator. Also, we often as in this problem consider only
x
steadystate behaviour, where the time derivatives become negligible: T 0. In t
44
45
such cases, we have the much simpler onedimensional form 0 k2T
x2
In the more general case of two dimensions, this becomes,
2T 2T 0kx2y2 ,
which appears in a wide range of physical contexts. Here, it describes the varia tion in temperature T in a plate with no internal generation of heat. Since the conductivity k is normally constant and nonzero, then the second term must be zero. The vector form of this simpler equation is generally known as Laplaces equation,
2T 2T 2T 0. 46 x2 y2
42
B Stream function descriptions of incompress ible flow
Incompressible fluid is a fluid whose volume is a conserved quantity, you cannot
squeeze the fluid into a smaller space think water!. This means that in one
dimension think pipe! you can measure the flux of volume same as the mass
flux for a density of 1 at one end of an inelastic pipe, and it must exactly match
the flux of volume at the other. Flux of volume Q and average flow rate U are
connected, Q U A, where A is the crosssectional area of the pipe. If our pipe
has length 2L between L and L, then we can say that QL QL 0. Now
think about the solution to this integral: L Q dx. So what were really saying L x
by an incompressible onedimensional flow is that,
L Qdx0. 47
L x
In two dimensions it is only a little more complicated, but instead of thinking of a pipe, we think of a square box where flow can only pass through its faces at rightangles to them. Again incompressibility means that the total amount of stuff in the box is constant, but for every new bit that flows in, an old bit must flow out. This time, however, there is no constraint on the direction that fluid can flow in or flow out, just that the total volume mustnt change. We can extend the integral description above to two dimensions by noting that we need to define a vector version of Q. We could use U and V for the x and y velocities respectively, and we should clarify whether we are talking about face areas in the vertical or the horizontal, so Q U Ax V Ay . In this context we can choose the shape of our box, so lets make it square, Ax Ay A const., and with sides 2L as before. Thus we could write,
L Q
L Q
y dy 0. 48 weve defined through our box, Q doesnt vary in y and Q doesnt vary in x, so
x dx
Theres not much we can do with this equation until we realise that with the flow
L
L
x y
adding terms like L Qdy, simply adds zero to the equation. However it helps
L x
us group the terms meaningfully,
L Q Q
L Q Q
xy dy0. 49 Following from the fundamental theorem of calculus, if we let our box tend to zero
size without changing its shape, then we could disregard the integrals altogether, 43
L
xy dx
L
Figure 8: Illustrative plot of stream function and corresponding streamlines in a uniform flow from the left.
and write down a differential form,
Q Q 0. 50
x y
From the way weve defined the flow in our box, Q AU and Q AV . The
area A is a common factor which tends towards but never reaches zero, so we could say that in two dimensions we satisfy the following equation,
U V 0, 51 x y
which is called the Continuity Equation for incompressible flow.
We need to introduce another concept that of the stream function. Steady incompressible flow in two dimensions has a special property that its velocity, a vector field quantity, can be described completely by a separate scalar field, a single number whose value varies throughout the domain. Returning to our square box, this single number is a relative measure of the volume flux through a surface or a bit of a surface. It needs a datum somewhere, usually a bottom lefthand corner is set to zero. For every unit of distance up the left side of our box from the bottom corner to a variable point Y , it contributes a small amount to the
stream function :
Y
If we let Y L, then we recover the volume flux through the whole left side of the
L
for a uniform flow is illustrated in figure 8. Another contribution comes if we now
44
x y
Udy 52 box, L U dy U Ax as introduced earlier. The accumulation of stream function
Y L
move from top left to top right, and to establish a consistent sign convention, we say postive volume flux is V 0 flow out the positive yaxis direction just as postive volume flux is U 0 flow in the postive xaxis direction. So we have,
L L
Udy V dx. 53
L L
One of the interesting features of the stream function is that provided the flow is incompressible, it doesnt matter how one moves from top left to bottom right, it can be any curve so long as we do this without forming a loop or a knot; the volume flux is still how much stuff moves from one side of the line to the other and the difference in values of at bottom left and top right reflect this. These values are not dependent on the curve drawn between them.
Letting our box once more tend to zero size without changing shape, we obtain the differential form
Uy V x. 54 If we differentiate the stream function, then since xdifferentiation has no effect
on y and viceversa, we have,
x
y
we dont like the sign convention? We could invent another single field that is
U
everywhere oriented at 90 degrees to : the vector V is always perpendicular
V 55
U, 56 thus a single field can describe both components of velocity. But what if
to
V
U . Lets call our new field , defined as follows:
x
U 57
V, 58 So what constraints are there on ? We need to ensure incompressibility, so we
y
substitute expressions for U and V into the continuity equation 51,
UV 0 59 x y
xxyy 0 60
2 2 0. 61 x2 y2
45
The vector form of this simpler equation is generally known as Laplaces equation,
2 0, 62
that crops up very often in all sorts of physical problems. Here, is in fact the symbol traditionally used for this quanitity, and is known as the velocity potential. We can show that the stream function also satisfies incompressibility, by substituting as follows, and recognising that for smooth fields the order of the
derivatives in x and y can be swapped:
U V
xyx y y x 0 63 If we substitute U and V for , then,
2
x x y y xVyU, 64
and by comparing with the vorticity vector in two dimensions, we have the fol lowing:
i j k 0
u 0 0 , 65
xy VU UV0 x y
so we can reduce this to the following simple equation 2
66
For the condition 2 0 the equation we are actually solving to be true, we require that 0 everywhere as an additional constraint. This is called the condition of irrotationality and only holds exactly true for straight flow in free space. However we tell ourselves that this holds approximately everywhere because it makes obtaining solutions a bit easier, and we cross our fingers and hope that the error isnt too large. In this problem we will aim to solve 2 0, but be aware that the flows it produces are not always physically meaningful. You should identify in your report where you believe the flows are representative, and where they are not.
C Derivation of thin membrane mechanics
To a good approximation, thin elastic membranes can be modelled as supporting tensile stress but have no resistance to compression and consequently no nega tive strain, and also cannot resist bending moment or shear. If any one part of
46
Figure 9: Discrete analog of an elastic string.
the elastic material locally undergoes only small displacements as a result of in plane tensile stress, then a linear elastic model is appropriate, even if globally, displacements from the unstressed state are large. If initially we consider a one dimensional elastic material think elastic band! as a slightly simpler analogue, then we can derive an equation of motion relating displacements u to axial forces F, and material stiffness k. To make it easier to conceptualise, we pretend that the mass of the elastic is lumped at regular intervals, and in between are exten sional springs, as depicted in figure 9. We derive the equation of motion for the mass at a general position x h. Newtons second law states that,
F mu , 67 and in the absence of local applied forcing, the force applied to the mass is due to
tensile stress in the elastic springs, ie.,
mu F FrightFleft 68 kux2huxhkuxuxh 69 kux2h2uxhux, 70
noting that extension of the right spring leads to a rightward force on the mass at x h, but extension of the left spring leads to a leftward negative force. Now we need to account for contributions from all the point masses, and the cumulative stiffness and length of all the springs put together. Thus we say the total length L is given by L Nh, where N is the total number of point masses. The total mass
47
M Nm, and the total stiffness K kN because as we introduce more and more springs, any displacement is distributed amongst the springs, so the ratio of applied forcing to strain response, ie. the stiffness, is reduced. We know that the local dynamical system for the mass at x h will oscillate in simple harmonic motion with a natural frequency given by,
2 k , 71 m
but we need to derive a new frequency to describe the whole system. By rear ranging the expressions for M, L and K, we have,
k NK 2K L2K
m M NM h M. 72
N
Thus by dividing equation 70 by m we can write,
u L2 K ux 2h 2ux h ux 73 M h2
and following the fundamental theorem of calculus we can let h 0. We know
that the numerator takes the form of a second derivative of u evaluated at location
xh,soifcLK wecansay, M
2u c2 2u, 74 t2 x2
for onedimensional inplane elastic motion with small strain. The quantity c has dimensions of velocity and is interpreted as a wave speed, with K a timescale
M
that represents the natural frequency of the system as a whole. This equation is also valid for displacements in the direction normal to the plane that result in small angular deviation of the elastic eg. a drum skin, or in this case a trampoline, since the local behaviour is still dominated by small strain. In the more general case of two dimensions, this becomes
2u 2 2u 2u
t2 c x2y2 , 75
a form of wave equation that appears in a wide range of physical contexts. Wave
equations of this sort describe the vibration of membranes such as drum skins,
or in this case, of a trampoline. In this problem we will be interested in the
equilibrium stress on the torn trampoline, when 2u 0 in a steady state. In t2
this case, then since c is not normally zero, equation 75 reduces to Laplaces equation,
2u 2u 2u 0. 76 x2 y2
48
D Numerical convergence
To understand why larger grid sizes lead to slower convergence, we need to express the numerical method as operations in linear algebra. Taking the 1D equation for a rope bridge for simplicity,
2 mg, 77 x2 k
then we treat the unknown values of on the grid as a vector , one element
for each grid point. We also treat the loading as a vector of the same length,
where in this case every element is equal. The operator 2 can be discretised and x2
represented as a square matrix with the same number of rows and columns as the vector has elements. Using our normal discretisation,
2 i1 2i i1 78 x2 x2
we can separate out the coefficients in front of from itself, and put those coefficients in the matrix that is otherwise zero. Effectively we are saying that the action of performing a doublederivative on any vector of values can be written down as a matrix in the following way:
.. .. ..
2… 1 2 1
2 2 2 . 79 x2 x x x
.. .. .. …
Thus, the linear algebra problem that is approximately equivalent to the physical statement for the rope bridge is given by,
. .
i1 .. .. ..i1
.. .. ..
mg
. . .
121 mkg
x2 x2 x2 i mkg , 80
which we might also write in a more compact form as A f. Direct inversion, A1f is in principle possible provided we check that A has no zero eigen values and is therefore invertible, but using a direct method to compute A1, then for every N grid points it takes ON3 arithmetic operations to compute the exact solution. Faster methods tend to use iteration to convergence and are less
…. k. . .
49
exact, but the very best methods can reduce the cost to ON, so the benefit is worthwhile on larger problems relative to the loss of accuracy. The pathway to these fast methods is described in 4.2. Below we discuss the principle behind all iterative approaches to solving this and many other problems.
The iteration process can also be interpreted as a matrix operation on vectors. This time, we need to think of the matrix as a tool for moving a vector from point ing in one particular direction towards another, which is what matrices do when operating on vectors. We start with an initial guess, and move this initial guess around in space until it converges. The elementwise operation being performed during an iteration of the rope bridge case is the following,
old old mg
new i1 i1 x2 , 81
i 2 2k and expressed for all elements we have,
. .
. . .
.oldnew . i1i1
old new mgx2 i i
.. .. 1
1
0 old new , 82 22 2k
.. .. . i1 i1 …..
. . . 1 1
which can be thought of as a sequential scaling and rotation of the vector the
modified vector composed of with a 1 added on the end by the matrix that
we might call B. We need the extra column on the matrix and the extra element
on the vector to account for the loading mg x2. Repeatedly applying B to a 2k
vector has the following structure,
1 B 0
2 B 1
3 B 2
4 B 3
5 B 4 .
83
The general case is k1 Bk, and when convergence is reached then the solution changes by an imperceptible amount, so by definition k1 k. This
50
is just a special case of an eigenvalue problem where the eigenvalue happens to be 1:
1. B. 84
so the converged solution is simply an eigenvector associated with the matrix B constructed from the numerical method. In fact the only way to compute eigenvalues and eigenvectors in large matrices with 6 or greater rows and columns is to iteratively change a vector from a starting location to convergence. The only additional trick is that when the eigenvalue is not 1, the vector length needs to be rescaled, and the amount by which it needs to be rescaled is the eigenvalue associated with that eigenvector. This is called a power iteration. The choice of initial vector is purely arbitrary, provided no eigenvalues of B are larger than 1 then the numerical method will converge to the principal eigenvector, and this is the solution of the problem.
What speeds up or slows down a numerical method is how quickly the principal eigenvector can be isolated from the others, and this is related to the relation ship between its eigenvalue and all the other eigenvalues. The largest eigenvalue tends to dominate, so when the ratio of max, say, to the next most important eigenvalue is large, isolating the dominant eigenvector happens quickly. When the eigenvalues are similar and the ratio is small, the solution takes much longer to converge. The solution will never converge if the dominant eigenvalue exceeds 1.
As matrices get larger in rows and columns, even if their structure doesnt change, their coefficients do if they represent the same physical problem. Many coefficients are dependent on x and these coefficients become smaller as x 0, changing the relative importance of the others. It turns out that the eigenvalues tend to cluster closer together, and closer to 1 in magnitude, so it becomes pro gressively more difficult to converge as it takes longer to separate out the dominant eigenvector. Consequently, as grid resolution increases, the cost of converging to the same tolerance grows as ON2 in the number of points, and not ON as one might have originally expected.
51