numerical analysis chapra 5ESM30.docx

资源描述

numerical analysis chapra 5ESM30.docx

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

numerical analysis chapra 5ESM30.docx

numericalanalysischapra5ESM30

CHAPTER30

30.1ThekeytoapproachingthisproblemistorecastthePDEasasystemofODEs.Thus,bysubstitutingthefinite-differenceapproximationforthespatialderivative,wearriveatthefollowinggeneralequationforeachnode

Bywritingthisequationforeachnode,thesolutionreducestosolving4simultaneousODEswithHeun’smethod.Theresultsforthefirsttwostepsalongwithsomelaterselectedvaluesaretabulatedbelow.Inaddition,aplotsimilartoFig.30.4,isalsoshown

x=0

x=2

x=4

x=6

x=8

x=10

100

0.1

100

2.043923

0.021788

0.010894

1.021962

0.2

100

4.005178

0.084022

0.042672

2.002593

100

37.54054

10.27449

6.442319

18.95732

100

53.24294

24.66052

17.4603

27.92251

100

62.39032

36.64937

27.84901

34.34692

100

68.71331

46.03498

36.54213

39.5355

30.2Becausewenowhavederivativeboundaryconditions,theboundarynodesmustbesimulated.Fornode0,

（i）

Thisintroducesanexteriornodeintothesolutionati=1.Thederivativeboundaryconditioncanbeusedtoeliminatethisnode,

whichcanbesolvedfor

whichcanbesubstitutedintoEq.（i）togive

Forourcase,dT0/dx=1andx=2,andthereforeT1=T1–4.ThiscanbesubstitutedintoEq.（i）togive,

Asimilaranalysiscanbeusedtoembedthezeroderivativeintheequationforthenthnode,

（ii）

Thisintroducesanexteriornodeintothesolutionatn+1.Thederivativeboundaryconditioncanbeusedtoeliminatethisnode,

whichcanbesolvedfor

whichcanbesubstitutedintoEq.（ii）togive

Forourcase,n=5anddTn/dx=0,andtherefore

Togetherwiththeequationsfortheinteriornodes,theentiresystemcanbesolvedwithastepof0.1s.Theresultsforsomeoftheearlystepsalongwithsomelaterselectedvaluesaretabulatedbelow.Inaddition,aplotofthelaterresultsisalsoshown

x=0

x=2

x=4

x=6

x=8

x=10

25.0000

0.1

24.9165

25.0000

0.2

24.8365

24.9983

25.0000

0.3

24.7597

24.9949

25.0000

0.4

24.6861

24.9901

24.9999

25.0000

0.5

24.6153

24.9840

24.9997

25.0000

200

5.000081

6.800074

8.200059

9.200048

9.800042

10.00004

400

-11.6988

-9.89883

-8.49883

-7.49882

-6.89881

-6.69881

600

-28.4008

-26.6008

-25.2008

-24.2008

-23.6007

-23.4007

800

-45.1056

-43.3056

-41.9056

-40.9056

-40.3056

-40.1056

1000

-61.8104

-60.0104

-58.6104

-57.6104

-57.0104

-56.8104

Noticewhat’shappening.Therodneverreachesasteadystate,becauseoftheheatlossattheleftend（unitgradient）andtheinsulatedcondition（zerogradient）attheright.

30.3Thesolutionfort=0.1is（ascomputedinExample30.1）,

x=0

x=2

x=4

x=6

x=8

x=10

100

0.1

100

2.0875

1.04375

0.2

100

4.087847

0.043577

0.021788

2.043923

Fort=0.05,itis

x=0

x=2

x=4

x=6

x=8

x=10

100

0.000000

0.05

100

1.043750

0.000000

0.521875

0.1

100

2.065712

0.010894

0.005447

1.032856

0.15

100

3.066454

0.032284

0.016228

1.533227

0.2

100

4.046528

0.063786

0.032229

2.023265

Toassessthedifferencesbetweentheresults,weperformedthesimulationathirdtimeusingamoreaccurateapproach（theHeunmethod）withamuchsmallerstepsize（t=0.001）.Itwasassumedthatthismorerefinedapproachwouldyieldapredictionclosetotruesolution.ThesevaluescouldthenbeusedtoassesstherelativeerrorsofthetwoEulersolutions.Theresultsaresummarizedas

x=0

x=2

x=4

x=6

x=8

x=10

Heun（h=0.001）

100

4.006588

0.083044

0.042377

2.003302

Euler（h=0.1）

100

4.087847

0.043577

0.021788

2.043923

ErrorrelativetoHeun

2.0%

47.5%

48.6%

2.0%

Euler（h=0.05）

100

4.046528

0.063786

0.032229

2.023265

ErrorrelativetoHeun

1.0%

23.2%

23.9%

1.0%

Notice,thataswouldbeexpectedforEuler’smethod,halvingthestepsizeapproximatelyhalvestheglobalrelativeerror.

30.4TheapproachdescribedinExample30.2mustbemodifiedtoaccountforthezeroderivativeattherighthandnode（i=5）.Todothis,Eq.（30.8）isfirstwrittenforthatnodeas

（i）

Thevalueoutsidethesystem（i=6）canbeeliminatedbywritingthefinitedifferencerelationshipforthederivativeatnode5as

whichcanbesolvedfor

Forourcase,dT/dx=0,soT6=T4andEq.（i）becomes

Thus,thesimultaneousequationstobesolvedatthefirststepare

whichcanbesolvedfor

Forthesecondstep,theright-handsideismodifiedtoreflectthesecomputedvaluesofTatt=0.1,

whichcanbesolvedfor

30.5ThesolutionisidenticaltoExample30.3,butwith9interiornodes.Thus,thesimultaneousequationstobesolvedatthefirststepare

whichcanbesolvedfor

Forthesecondstep,theright-handsideismodifiedtoreflectthesecomputedvaluesofTatt=0.1,

whichcanbesolvedfor

30.6UsingtheapproachfollowedinExample30.5,Eq.（30.20）isappliedtonodes（1,1）,（1,2）,and（1,3）toyieldthefollowingtridiagonalequations

whichcanbesolvedfor

Inasimilarfashion,tridiagonalequationscanbedevelopedandsolvedfor

and

Forthesecondsteptot=10,Eq.（30.22）isappliedtonodes（1,1）,（2,1）,and（3,1）toyield

whichcanbesolvedfor

Tridiagonalequationsfortheotherrowscanbedevelopedandsolvedfor

and

Thus,theresultattheendofthefirststepcanbesummarizedas

i=0

i=1

i=2

i=3

i=4

j=4

120

j=3

13.46293

9.820288

12.86944

j=2

5.050756

0.815685

4.29281

j=1

4.47395

0.416205

3.737833

j=0

Thecomputationcanberepeatedandtheresultsfort=2000sareshownbelow:

i=0

i=1

i=2

i=3

i=4

j=4

120

j=3

80.71429

83.83929

77.14286

j=2

59.01786

57.5

54.73214

j=1

37.85714

32.41071

34.28571

j=0

30.7Althoughthisproblemcanbemodeledwiththefinite-differenceapproach（seeSec.32.1）,thecontrol-volumemethodprovidesamorestraightforwardwaytohandletheboundaryconditions.Thus,thetankisidealizedasaseriesofcontrolvolumes:

Theboundaryfluxesandthereactiontermcanbeusedtodevelopthediscreteformoftheadvection-diffusionequationfortheinteriorvolumesas

ordividingbothsidesbyx,

whichispreciselytheformthatwouldhaveresultedbysubstitutingcenteredfinitedifferenceapproximationsintotheadvection-diffusionequation.

Forthefirstboundarynode,nodiffusionisalloweduptheentrancepipeandadvectionishandledwithabackwarddifference,

ordividingbothsidesbyx,

Forthelastboundarynode,nodiffusionisallowedthroughtheexitpipeandadvectionoutofthetankisagainhandledwithabackwarddifference,

ordividingbothsidesbyx,

Bywritingtheseequationsforeachequally-spacedvolume,thePDEistransformedintoasystemofODEs.ExplicitmethodslikeEuler’smethodorotherhigher-orderRKmethodscanthenbeusedtosolvethesystem.

Theresultswithaninitialconditionthatthereactorhaszeroconcentrationwithaninflowconcentrationof100（usingEulerwithastepsizeof0.005）fort=100are

0.5

1.5

2.5

3.5

4.5

5.5

6.5

7.5

8.5

9.5

42.0320

41.5128

41.0509

40.6463

40.2989

40.0087

39.7760

39.6008

39.4836

39.4248

Aplotoftheresultsisshownbelow:

30.8HereisaVBAprogramthatimplementstheexplicitmethod.ItissetuptoduplicateExample30.1.

OptionExplicit

SubExplicit（）

DimiAsInteger,jAsInteger,npAsInteger,nsAsInteger

DimTe（20）AsDouble,dTe（20）AsDouble,tpr（20）AsDouble,Tepr（20,20）AsDouble

DimkAsDouble,dxAsDouble,LAsDouble,tcAsDouble,tfAsDouble

DimtpAsDouble,tAsDouble,tendAsDouble,hAsDouble,tolAsDouble

DimxAsDouble

tol=0.000001

L=10

ns=5

dx=L/ns

k=0.835

Te（0）=100

Te（5）=50

tc=0.1

tf=12

tp=3

np=0

tpr（np）=t

Fori=0Tons

Tepr（i,np）=Te（i）

Nexti

tend=t+tp

Iftend>tfThentend=tf

h=tc

Ift+h>tendThenh=tend-t

CallDerivs（Te,dTe,ns,dx,k）

Forj=1Tons-1

Te（j）=Te（j）+dTe（j）*h

Nextj

t=t+h

Ift>=tendThenExitDo

Loop

np=np+1

tpr（np）=t

Forj=0Tons

Tepr（j,np）=Te（j）

Nextj

Ift+tol>=tfThenExitDo

Loop

Sheets（"sheet1"）.Select

Range（"a4:

bb5000"）.ClearContents

Range（"a4"）.Select

ActiveCell.Value="time"

x=0

Forj=0Tons

ActiveCell.Offset（0,1）.Select

ActiveCell.Value="x="&x

x=x+dx

Nextj

Range（"a5"）.Select

Fori=0Tonp

ActiveCell.Value=tpr（i）

Forj=0Tons

ActiveCell.Offset（0,1）.Select

ActiveCell.Value=Tepr（j,i）

Nextj

ActiveCell.Offset（1,-ns-1）.Select

Nexti

Range（"a5"）.Select

EndSub

SubDerivs（Te,dTe,ns,dx,k）

DimjAsInteger

Forj=1Tons-1

dTe（j）=k*（Te（j-1）-2*Te（j）+Te（j+1））/dx^2

Nextj

EndSub

Whentheprogramisrun,theresultis

30.9ThisVBAprogramissetuptoeitheruseDirichletorgradientboundaryconditionsdependingonthevaluesoftheparametersistrtandiend.ItissetuptosolvethefirstfewstepsofProb.30.2.

OptionExplicit

SubEulerPDE（）

DimiAsInteger,jAsInteger,npAsInteger,nsAsInteger

Dimis

展开阅读全文