A Simple Finite-Difference Approach

<BODY LANG="EN"> <HR><center> <table border=1> <TR> <td><A HREF="node11.shtml"><img height=36 width=36 alt="PREVIOUS" src="/UCESIcons/tools/previous.gif"></A></td> <td><A HREF="node11.shtml"><img height=36 width=36 alt="UP" src="/UCESIcons/tools/up.gif"></A></td> <td><A HREF="node13.shtml"><img height=36 width=36 alt="NEXT" src="/UCESIcons/tools/next.gif"></A></td> <td><a href="contents-node12.html"><img height=36 width=36 alt="CONTENTS" src="/UCESIcons/tools/toc.gif"></a></td> <td><a href="reset.hamlet"><img height=36 width=36 alt="RESET" src="/UCESIcons/tools/stop.gif"></a></td> <td><a href="copyright.html"><img height=36 width=36 alt="COPYRIGHT" src="/UCESIcons/tools/copyright.gif"></a></td> <td><a href="mailto:uces_info@krellinst.org"><img height=36 width=36 alt="FEEDBACK" src="/UCESIcons/tools/feedback.gif"></a></td> <td><a href="http://www.cs.utah.edu/%7Ehamlet/icons.html"><img height=36 width=36 alt="ABOUT" src="/UCESIcons/tools/about.gif"></a></td> </TR> </table> </center> <HR> <!DOCTYPE HTML PUBLIC "-//IETF//DTD HTML 2.2//EN">  <HTML> <HEAD> <TITLE>A Simple Finite-Difference Approach</TITLE> <BASE TARGET="_parent"> <META NAME="description" CONTENT="A Simple Finite-Difference Approach"> <META NAME="keywords" CONTENT="cone"> <META NAME="resource-type" CONTENT="document"> <META NAME="distribution" CONTENT="global"> <LINK REL=STYLESHEET HREF="cone.css"> </HEAD> <BODY LANG="EN"> <H3>A Simple Finite-Difference Approach</H3> <P> Let us, instead, seek a <EM>numerical</EM> solution by replacing the derivatives with <EM>finite differences</EM>. In a way, this is a step backwards in mathematics, since we are going to replace an equation from <EM>calculus</EM> with one which is purely algebraic. Recall from your first course in calculus that you crossed the border between algebra and calculus by <EM>taking a limit</EM> to define a derivative. That is, the derivative <nobr><I>df</I>/<I>dx</I></nobr> was defined as: <P> <IMG WIDTH=361 HEIGHT=35 ALIGN=BOTTOM ALT="[limit definition]" SRC="img56.gif"> <P> So, we might approximate <I>dy</I>/<I>dt</I> in our equation as: <P> <IMG WIDTH=350 HEIGHT=35 ALIGN=BOTTOM ALT="[limit approx.]" SRC="img57.gif"> <P> where <IMG WIDTH=11 HEIGHT=12 ALIGN=BOTTOM ALT="delta_t" SRC="img58.gif"> is `small.' Rearranging our <EM>differential equation</EM> slightly gives: <P> <IMG WIDTH=529 HEIGHT=43 ALIGN=BOTTOM ALT="[Eq. 20]" SRC="img59.gif"> <P> Again the time-dependent major radius of the cone at the fluid level is given by: <P> <IMG WIDTH=339 HEIGHT=38 ALIGN=BOTTOM ALT="[Eq. for R(y)]" SRC="img45.gif"> <P> If we are willing to approximate <nobr><I>dy</I>/<I>dt</I></nobr> with the <EM>finite differences</EM> <IMG WIDTH=36 HEIGHT=24 ALIGN=MIDDLE ALT="delta_y/delta_t" SRC="img60.gif"> , then <P> <IMG WIDTH=560 HEIGHT=43 ALIGN=BOTTOM ALT="[Eq. 21]" SRC="img61.gif"> <P> is the equation we need to solve on the computer. Notice that Equation (21) is <EM>non-linear</EM>; that is, the right-hand side contains the solution <I>y</I> and the square of the solution  <nobr><I>y</I><sup><small>2 </small></sup>.</nobr> Note also that just as the vessel is about to empty, the radical <IMG WIDTH=109 HEIGHT=31 ALIGN=MIDDLE ALT="[sqrt expression]" SRC="img63.gif"> will go to zero if we are not careful to require that <IMG WIDTH=42 HEIGHT=18 ALIGN=MIDDLE ALT="r > r_e" SRC="img64.gif"> , since <IMG WIDTH=65 HEIGHT=23 ALIGN=MIDDLE ALT="R(y) -> r" SRC="img65.gif"> as <IMG WIDTH=40 HEIGHT=22 ALIGN=MIDDLE ALT="y -> 0" SRC="img66.gif"> . <P> Now there is clearly some ambiguity about what value of the fluid height <I>y</I> we should use in the right-hand side of Equation (21). Two approaches are obvious. <P> <UL> <LI> We could evaluate <IMG WIDTH=59 HEIGHT=17 ALIGN=MIDDLE ALT="y(t+delta_t)" SRC="img68.gif"> completely in terms of <nobr><I>y</I>(<I>t</I>),</nobr> which is always known from the last time step. This includes evaluating the fluid surface radius <nobr><I>R</I>(<I>y</I>)</nobr> using the values computed at the previous time step.<P> <LI> We could <EM>iterate</EM> at each time level. That is, we could first solve for <IMG WIDTH=59 HEIGHT=17 ALIGN=MIDDLE ALT="y(t+delta_t)" SRC="img68.gif"> using the value of <I>y</I> from the previous time level. Then we could repeat the calculation several times without advancing the time level and presumably get better and better agreement between <IMG WIDTH=59 HEIGHT=17 ALIGN=MIDDLE ALT="y(t+delta_t)" SRC="img68.gif"> on the left-hand side and the value of <I>y</I> in the terms on the right-hand side. </UL> <P> This second approach is the one we should take. No matter what programming language one might choose to use, a program to obtain numerical solutions to Equation (21) should take roughly the following form: <OL> <LI> Specify physical parameters <IMG WIDTH=59 HEIGHT=22 ALIGN=MIDDLE ALT="R_O, r, y_0," SRC="img69.gif"> and <I>g</I> in a consistent set of units. <LI> Specify a time step size <IMG WIDTH=11 HEIGHT=12 ALIGN=BOTTOM ALT="delta_t" SRC="img58.gif"> . <LI> Set initial values for the dependent variable <nobr><I>y</I>(<I>t</I>)</nobr> and <nobr><I>R</I>(<I>y</I>).</nobr> <LI> Solve the <EM>finite-difference</EM> Equation (21) for <IMG WIDTH=59 HEIGHT=17 ALIGN=MIDDLE ALT="y(t+delta_t)" SRC="img68.gif"> , using values of <nobr><I>y</I>(<I>t</I>)</nobr> in the right-hand side. <LI> Repeat the previous step using the values of <IMG WIDTH=59 HEIGHT=17 ALIGN=MIDDLE ALT="y(t+delta_t)" SRC="img68.gif"> just computed to evaluate the right-hand side. <LI> <EM>Iterate</EM> by repeating this last step until the value of <IMG WIDTH=59 HEIGHT=17 ALIGN=MIDDLE ALT="y(t+delta_t)" SRC="img68.gif"> stops changing significantly. <LI> Advance the solution another time step by beginning with Step 3. <LI> Stop and write out the solution when <IMG WIDTH=59 HEIGHT=17 ALIGN=MIDDLE ALT="y(t+delta_t)" SRC="img68.gif"> is sufficiently close to zero, or whatever the final value is. </OL> <P> Below is a simple program written in Fortran which solves our equation in the manner outlined above. <P> <PRE> program cone parameter (imax=1500) ! Set the maximum dimensions. dimension y(imax),t(imax),r(imax),dydt(imax) open(1,file='cone1.plt',form='formatted')! Open a file for plotting. c Below we compute 'machine epsilon,' the smallest number your c computer can add to one and get something different from one. c This is very useful for testing for 'machine zero', etc. eps=1. 100 eps=eps/2. if(eps +1. .gt. 1.)goto 100 print *,'Machine epsilon: ',eps c Read constant parameters. c Here they have been set at values sufficient for evaluating c the behavior of the system. print *, 'Enter R0,Y0,ro,re, and g in that order.' R0=1. ! major radius of the cone Y0=1. ! initial fluid level ro=.1 ! minor radius of the cone re=.01 ! radius of the drain g=9.81 ! acceleration of gravity c read(5,*)R0,Y0,ro,re,g ! un-comment this line to change the c fixed values. dr=R0-ro! This is just a convenient definition to shorten the c Fortran statements below. c Read computational parameters. Here we need to choose a time c step size, dt. You will notice that the smaller dt, the more c nearly you can compute y=0, but the price for small dt is c a huge data file, since many time steps will be required to empty c the vessel. print *, 'Enter dt and itmax in that order.' read(5,*)dt,itmax c c Initialize. Here we set values at t=0. y(1)=Y0 yold=Y0 t(1)=0. r(1)=R0 i=1 c Begin time stepping. 1 i=i+1 ! This statement begins a time step. y(i)=y(i-1)! This gives a non-zero value of y to get the computation c started at each new time level. if(i .gt. imax)goto 2! Stop if we've exceeded the variable dimensions. t(i)=dt+t(i-1)! increment time do iter=1,itmax! begin the iteration loop. c Usually only 3 or 4 iterations are required at each time step. r(i)=ro+y(i)*dr/Y0 ! compute the radius of the fluid surface c The following statements just evaluate parts of Equation (21) term1=sqrt(2.*g*y(i))*re**2*Y0**2 ratio=re**4/r(i)**4 term2=sqrt(1.-ratio) term3=ro**2*Y0**2+2.*dr*ro*y(i)*Y0+ dr**2*y(i)**2 y(i)=y(i-1)-dt*term1/(term2*term3)! This is Equation (21) c Below is just a safety net to catch the solution should we c compute a non-physical negative fluid level as we work in the c neighborhood of y=0. if(y(i) .lt. 0.)then y(i)=0. r(i)=ro+y(i)*dr/Y0 goto 2 endif c The statement below just lets you monitor the progress. print *,'r(i),y(i),i ',r(i),y(i),i c Below tests to see if we are very close to y=0. if(y(i) .lt. sqrt(eps).and. y(i-1).lt.sqrt(eps))goto 2 c Below we test to see how closely two successive values of c y in the iteration agree. if(yold .gt. 0.)resid=abs((yold-y(i))/yold) yold=y(i)! yold is the value from the previous iteration if(resid .lt. sqrt(eps))goto 4! if we have converged, take another c time step enddo 4 goto 1 ! start another time step after the maximum number of c iterations regardless of whether the solution has converged 2 n=i c Below we compute the derivative of the fluid level. Clearly c dy/dt is the velocity of the fluid surface. Due to the c singular behavior of Equation (21) as y-->0, dy/dt becomes c infinite, so we don't try to compute the last couple of values. do i=2,n-1 dydt(i)=(y(i+1)-y(i-1))/dt dmax=max(abs(dydt(i)),dmax) enddo dydt(1)=(y(2)-y(1))/dt c Now we write out the solution for plotting. do i=1,n-1 dydt(i)=dydt(i)/dmax write(1,*)t(i),y(i),r(i),-dydt(i) enddo stop end</PRE> <P> A plot of the solution for the <B>cone</B> using this code is available <a href="cone_graph.html"><B>here</B></a>. Note that as the fluid level approaches zero, the rate at which it drops <nobr>(<I>dy</I>/<I>dt</I>)</nobr> increases rapidly, since the remain volume is decreasing very rapidly. <P> </BODY> </HTML> <HR><center> <table border=1> <TR> <td><A HREF="node11.shtml"><img height=36 width=36 alt="PREVIOUS" src="/UCESIcons/tools/previous.gif"></A></td> <td><A HREF="node11.shtml"><img height=36 width=36 alt="UP" src="/UCESIcons/tools/up.gif"></A></td> <td><A HREF="node13.shtml"><img height=36 width=36 alt="NEXT" src="/UCESIcons/tools/next.gif"></A></td> <td><a href="contents-node12.html"><img height=36 width=36 alt="CONTENTS" src="/UCESIcons/tools/toc.gif"></a></td> <td><a href="reset.hamlet"><img height=36 width=36 alt="RESET" src="/UCESIcons/tools/stop.gif"></a></td> <td><a href="copyright.html"><img height=36 width=36 alt="COPYRIGHT" src="/UCESIcons/tools/copyright.gif"></a></td> <td><a href="mailto:uces_info@krellinst.org"><img height=36 width=36 alt="FEEDBACK" src="/UCESIcons/tools/feedback.gif"></a></td> <td><a href="http://www.cs.utah.edu/%7Ehamlet/icons.html"><img height=36 width=36 alt="ABOUT" src="/UCESIcons/tools/about.gif"></a></td> </TR> </table> </center> <HR> <P> <ADDRESS> <I>Boyd Gatlin<BR> Associate Professor of Aerospace Engineering<BR> Mississippi State University</I>  </ADDRESS> </BODY>