numerical analysis chapra 5ESM31Word格式.docx
- 文档编号:21082089
- 上传时间:2023-01-27
- 格式:DOCX
- 页数:19
- 大小:418.01KB
numerical analysis chapra 5ESM31Word格式.docx
《numerical analysis chapra 5ESM31Word格式.docx》由会员分享,可在线阅读,更多相关《numerical analysis chapra 5ESM31Word格式.docx(19页珍藏版)》请在冰豆网上搜索。
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
k=ii-1+i
m=ii-1+j
a(k,m)=a(k,m)+s(i,j)
b(k)=b(k)+ff*((x(ii+1)-x(ii))-(x(ii+1)-x(ii))/2)
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)
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
a3"
ActiveCell.Value="
dTe(x="
&
x(0)&
"
)/dx="
ActiveCell.Offset(0,1).Select
ActiveCell.Value=dTeLeft
ActiveCell.Offset(1,-1).Select
x(ns+1)&
ActiveCell.Value=dTeRight
ActiveCell.Offset(3,-1).Select
ActiveCell.Value=x(i)
ActiveCell.Offset(0,1).Select
ActiveCell.Value=Te(i)
ActiveCell.Offset(1,-1).Select
b3"
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
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
r(k)=r(k)-e(k)*r(k-1)
x(n)=r(n)/f(n)
Fork=n-1To1Step-1
x(k)=(r(k)-g(k)*x(k+1))/f(k)
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);
holdon
31.11
ProgramPlate
UseIMSL
ImplicitNone
Integer:
:
ncval,nx,nxtabl,ny,nytabl
Parameter(ncval=11,nx=33,nxtabl=5,ny=33,nytabl=5)
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)
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
EndProgram
Functionprhs(x,y)
prhs,x,y
prhs=0
EndFunction
RealFunctionbrhs(iside,x,y)
iside
x,y
If(iside==1)Then
brhs=50
ElseIf(iside==2)Then
brhs=0
ElseIf(iside==3)Then
brhs=75
Else
brhs=100
EndIf
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:
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
Element3:
Node4:
Element4:
Node5:
Element5:
Node6:
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]'
7.5982
164.1913
200.6404
201.9494
157.3059
-14.8037
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- numerical analysis chapra 5ESM31 ESM31