numerical analysis chapra 5ESM31.docx

numerical analysis chapra 5ESM31.docx

《numerical analysis chapra 5ESM31.docx》由会员分享，可在线阅读，更多相关《numerical analysis chapra 5ESM31.docx（19页珍藏版）》请在冰豆网上搜索。

numericalanalysischapra5ESM31

CHAPTER31

31.1Theequationtobesolvedis

AssumeasolutionoftheformT=ax2+bx+cwhichcanbedifferentiatedtwicetogiveT”=2a.Substitutingthisresultintothedifferentialequationgivesa=7.5.Theboundaryconditionscanbeusedtoevaluatetheremainingcoefficients.Forthefirstconditionatx=0,

orc=75.Similarly,forthesecondcondition

whichcanbesolvedforb=82.5.Therefore,thefinalsolutionis

Theresultsareplottedas

31.2TheheatsourceterminthefirstrowofEq.（31.26）canbeevaluatedbysubstitutingEq.（31.3）andintegratingtogive

Similarly,Eq.（31.4）canbesubstitutedintotheheatsourcetermofthesecondrowofEq.（31.26）,whichcanalsobeintegratedtoyield

TheseresultsalongwiththeotherparametervaluescanbesubstitutedintoEq.（31.26）togive

and

31.3InamannersimilartoFig.31.7,theequationscanbeassembledforthetotalsystem,

Theknownendtemperaturescanbesubstitutedtogive

Theseequationscanbesolvedfor

Thesolution,alongwiththeanalyticalsolution（dashedline）isshownbelow:

31.4

（1）

（2）

（3）

Term

（1）:

Term

（2）:

Term（3）:

Totalelementequation[

（1）+

（2）+（3）]

where

31.5Firstwecandevelopananalyticalfunctionforcomparison.Substitutingparametersgives

Assumeasolutionoftheform

Thiscanbedifferentiatedtwicetoyieldd2u/dx2=2a.ThiscanbesubstitutedintotheODE,whichcanthenbesolvedfora=2.5108.Theboundaryconditionscanthenbeusedtoevaluatetheremainingcoefficients.Attheleftside,u（0）=0and

andtherefore,c=0.Attherightsideofthebar,u（10）=0and

andtherefore,b=2.5107,andthesolutionis

whichcanbedisplayedas

Theelementequationcanbewrittenas

Thedistributedloadcanbeevaluatedas

Thus,thefinalelementequationis

Assemblyyields

whichcanbesolvedfor

Theseresults,alongwiththeanalyticalsolution（dashedline）aredisplayedbelow:

31.6HereisaVBAprogramtosolveforthesteadystatetemperaturedistributionforaheatedrodwithaconstantheatsource.ItissetuptosolvethesystemdescribedinExamples31.1and31.2.

OptionExplicit

SubFErod（）

DimnsAsInteger,iiAsInteger,iAsInteger,jAsInteger

DimkAsInteger,mAsInteger

Dimx（100）AsDouble,st（2,2）AsDouble,cAsDouble

Dims（2,2）AsDouble,a（100,100）AsDouble,b（100）AsDouble,d（100）AsDouble

DimTe（100）AsDouble,ffAsDouble

Dime（100）AsDouble,f（100）AsDouble,g（100）AsDouble,r（100）AsDouble

Dimdum1AsDouble,dum2AsDouble

DimdTeLeftAsDouble,dTeRightAsDouble

'setparameters

ns=4

x

（1）=0

x

（2）=2.5

x（3）=5

x（4）=7.5

x（5）=10

Te

（1）=40

Te（5）=200

ff=10

'constructsystemmatrix

st（1,1）=1:

st（1,2）=-1:

st（2,1）=-1:

st（2,2）=1

Forii=1Tons

c=1/（x（ii+1）-x（ii））

Fori=1To2

Forj=1To2

s（i,j）=c*st（i,j）

Nextj

Nexti

Fori=1To2

k=ii-1+i

Forj=1To2

m=ii-1+j

a（k,m）=a（k,m）+s（i,j）

Nextj

b（k）=b（k）+ff*（（x（ii+1）-x（ii））-（x（ii+1）-x（ii））/2）

Nexti

Nextii

'determineimpactofuniformsourceandboundaryconditions

CallMmult（a,Te,d,ns+1,ns+1,1）

Fori=1Tons+1

b（i）=b（i）-d（i）

Nexti

a（1,1）=1

a（2,1）=0

a（ns+1,ns+1）=-1

a（ns,ns+1）=0

'Transformsquarematrixintotridiagonalform

f

（1）=a（1,1）

g

（1）=a（1,2）

r

（1）=b

（1）

Fori=2Tons

e（i）=a（i,i-1）

f（i）=a（i,i）

g（i）=a（i,i+1）

r（i）=b（i）

Nexti

e（ns+1）=a（ns+1,ns）

f（ns+1）=a（ns+1,ns+1）

r（ns+1）=b（ns+1）

'Tridiagonalsolver

dum1=Te

（1）

dum2=Te（ns+1）

CallTridiag（e,f,g,r,ns+1,Te）

dTeLeft=Te

（1）

dTeRight=Te（ns+1）

Te

（1）=dum1

Te（ns+1）=dum2

'outputresults

Sheets（"sheet1"）.Select

Range（"a7:

b5000"）.ClearContents

Range（"a3"）.Select

ActiveCell.Value="dTe（x="&x（0）&"）/dx="

ActiveCell.Offset（0,1）.Select

ActiveCell.Value=dTeLeft

ActiveCell.Offset（1,-1）.Select

ActiveCell.Value="dTe（x="&x（ns+1）&"）/dx="

ActiveCell.Offset（0,1）.Select

ActiveCell.Value=dTeRight

ActiveCell.Offset（3,-1）.Select

Fori=1Tons+1

ActiveCell.Value=x（i）

ActiveCell.Offset（0,1）.Select

ActiveCell.Value=Te（i）

ActiveCell.Offset（1,-1）.Select

Nexti

Range（"b3"）.Select

EndSub

SubMmult（a,b,c,m,n,l）

DimiAsInteger,jAsInteger,kAsInteger

DimsumAsDouble

Fori=1Ton

sum=0

Fork=1Tom

sum=sum+a（i,k）*b（k）

Nextk

c（i）=sum

Nexti

EndSub

SubTridiag（e,f,g,r,n,x）

DimkAsInteger

'decompose

Fork=2Ton

e（k）=e（k）/f（k-1）

f（k）=f（k）-e（k）*g（k-1）

Nextk

'substitute

Fork=2Ton

r（k）=r（k）-e（k）*r（k-1）

Nextk

x（n）=r（n）/f（n）

Fork=n-1To1Step-1

x（k）=（r（k）-g（k）*x（k+1））/f（k）

Nextk

EndSub

Theoutputis

31.7Aftersettinguptheoriginalspreadsheet,thefollowingmodificationswouldbemadetoinsulatetherightedgeandaddthesink:

CellI1:

Setto100

CellI2:

=（I1+2*H2+I3）/4;ThisformulawouldthenbecopiedtocellsI3:

I8.

CellI9:

=（I8+H9）/2

CellC7:

=（C6+D7+C8+B7-150）/4

Theresultingspreadsheetisdisplayedbelow:

Correspondingcontourplotscanbegeneratedas

31.8Theresultsoftheprecedingproblemcanbesavedasatab-delimitedtextfile（inourcase,wecalledthefileprob3108.txt）.ThefollowingcommandscanthenbeusedtoloadthisfileintoMATLAB,aswellastogeneratethecontourplotalongwithheatflowvectors.

>>loadprob3108.txt>>[px,py]=gradient（prob3108）;

>>cs=contour（prob3108）;clabel（cs）;holdon>>quiver（-px,-py）;holdoff

31.9TheschemeforimplementingthiscalculationonExcelisshownbelow:

ThesimpleLaplaceequationappliesdirectlytotheblankwhitecells（recallFig.31.14）.However,fortheshadedcellsimpactedbytheheatsink,aheatbalanceequationmustbewritten.Forexample,forcellE3,thecontrolvolumeapproachcanbedevelopedas

Collectingandcancelingtermsyields

Substitutingthelengthdimensionsandthecoefficientofthermalconductivitygives,

Thisformulacanbecopiedtotheothershadedcells.Theresultisdepictedbelow,alongwiththecorrespondingcontourplots.

31.10Theresultsoftheprecedingproblemcanbesavedasatab-delimitedtextfile（inourcase,wecalledthefileprob3110.txt）.ThefollowingcommandscanthenbeusedtoloadthisfileintoMATLAB,aswellastogeneratethecontourplotalongwithheatflowvectors.

>>loadprob3110.txt

>>[px,py]=gradient（prob3110）;

>>cs=contour（prob3110）;clabel（cs）;holdon

>>quiver（-px,-py）;holdoff

31.11

ProgramPlate

UseIMSL

ImplicitNone

Integer:

:

ncval,nx,nxtabl,ny,nytabl

Parameter（ncval=11,nx=33,nxtabl=5,ny=33,nytabl=5）

Integer:

:

i,ibcty（4）,iorder,j,nout

Real:

:

ax,ay,brhs,bx,by,coefu,prhs,u（nx,ny）,utabl,x,xdata（nx）,y,ydata（ny）

Externalbrhs,prhs

ax=0

bx=40

ay=0

by=40

ibcty

（1）=1

ibcty

（2）=2

ibcty（3）=1

ibcty（4）=1

coefu=0

iorder=4

CallFPS2H（prhs,brhs,coefu,nx,ny,ax,bx,ay,by,ibcty,iorder,u,nx）

Doi=1,nx

xdata（i）=ax+（bx-ax）*Float（i-1）/Float（nx-1）

EndDo

Doj=1,ny

ydata（j）=ay+（by-ay）*Float（j-1）/Float（ny-1）

EndDo

CallUMACH（2,nout）

Write（nout,'（8X,A,11X,A,11X,A）'）'X','Y','U'

Doj=1,nytabl

Doi=1,nxtabl

x=ax+（bx-ax）*Float（i-1）/Float（nxtabl-1）

y=ay+（by-ay）*Float（j-1）/Float（nytabl-1）

utabl=QD2VL（x,y,nx,xdata,ny,ydata,u,nx,.FALSE.）

Write（nout,'（4F12.4）'）x,y,utabl

EndDo

EndDo

EndProgram

Functionprhs（x,y）

ImplicitNone

Real:

:

prhs,x,y

prhs=0

EndFunction

RealFunctionbrhs（iside,x,y）

ImplicitNone

Integer:

:

iside

Real:

:

x,y

If（iside==1）Then

brhs=50

ElseIf（iside==2）Then

brhs=0

ElseIf（iside==3）Then

brhs=75

Else

brhs=100

EndIf

EndFunction

Output:

0.00000.000075.0000

10.00000.000071.6339

20.00000.000066.6152

30.00000.000059.1933

40.00000.000050.0000

0.000010.000075.0000

10.000010.000072.5423

20.000010.000067.9412

30.000010.000060.1914

40.000010.000050.0000

0.000020.000075.0000

10.000020.000075.8115

20.000020.000072.6947

30.000020.000064.0001

40.000020.000050.0000

0.000030.000075.0000

10.000030.000083.5385

20.000030.000083.0789

30.000030.000074.3008

40.000030.000050.0000

0.000040.000087.5000

10.000040.0000100.0000

20.000040.0000100.0000

30.000040.0000100.0000

40.000040.000075.0000

Pressanykeytocontinue

31.12

Element1:

Node1:

Node2:

Theothernodeequationscanbederivedsimilarly,

Element2:

Node2:

Node3:

Theotherelementequationsaresimilar.

Equationassembly:

Insertingboundaryconditions:

ThesolutioncanthenbeobtainedwithMATLAB:

>>C=[10-100000;

-1020-10000;

0-1020-1000;

00-1020-100;

000-10200;

0000-10-100];

>>b=[1253003003001300-850]';

>>x=C\b

x=

462.5000

450.0000

407.5000

335.0000

232.5000

-14.7500

31.13

Element1:

Node1:

Node2:

Theothernodeequationscanbederivedsimilarly,

Element2:

Node2:

Node3:

Element3:

Node3:

Node4:

Element4:

Node4:

Node5:

Element5:

Node5:

Node6:

Equationassembly:

Insertingboundaryconditions:

ThesolutioncanthenbeobtainedwithMATLAB:

>>C=[100-9.50000;

018-8.5000;

0-8.516-7.500;

00-7.514-6.50;

000-6.5120;

0000-5.5-50];

>>b=[-8001250300300575-125]';

>>x=C\b

x=

7.5982

164.1913

200.6404

201.9494

157.3059

-14.8037