Transcript Dimensioneringsuppgift

Ber¨akningsuppgift TMHL24 HT1 2014:
sv¨angnings- och utmattningsanalys
1
Inledning
Detta h¨afte beskriver den obligatoriska ber¨akningsuppgiften i kursen TMHL24
H˚
allfasthetsl¨ara-dimensioneringsmetoder h¨ostterminen 2014. Uppgiften genomf¨ors
individuellt med individuella data och rapporteras skriftligt. Rapporten ska
kunna l¨asas oberoende av detta h¨afte s˚
a det ska finnas beskrivningar och
f¨orklaringar till alla steg samt svar p˚
a de direkta fr˚
agor som st¨alls.
Ber¨akningarna ¨ar t¨ankta att genomf¨oras i Matlab. I slutet av h¨aftet finns
ett exempel p˚
a implementering av en sv¨angningsber¨akning i Matlab.
Matlab-koden ska bifogas rapporten och den ska vara l¨att att f¨olja och
utg˚
a direkt fr˚
an din dataupps¨attning, dvs om du inf¨or mellanresultat ska
dessa ber¨aknas direkt i Matlab-koden, inte p˚
a papper vid sidan av.
Rapporten l¨amnas in utskriven p˚
a papper till din lektionsledare. Om
rapporten l¨amnas in efter sista lektionstillf¨allet l¨amnas den ist¨allet i facket
utanf¨or Ulf Edlunds d¨orr. P˚
a kurshemsidan www.mechanics.iei.liu.se/edu ug/tmhl24
finns ett f¨ors¨attsblad som ska anv¨andas p˚
a rapporten. P˚
a detta f¨ors¨attsblad
finns ocks˚
a vissa instruktioner f¨or uppgiften och rapporten.
2
Ber¨
akningsuppgift
Figur 1: Strukturen som ska analyseras med avseende p˚
a deformation och
utmattning.
Strukturen i fig. 1 best˚
ar av en konsolbalk (fast insp¨and i ena ¨anden, fritt
r¨orlig i den andra), med tv˚
a punktmassor f¨astade. Sj¨alva balken antas approximativt vara massl¨os. Balken s¨atts i b¨ojsv¨angning av yttre krafter som
verkar p˚
a punktmassorna:
Q1 (t) = Q01 cos(ωp1t)
Q2 (t) = Q02 cos(ωp2t)
(1)
d¨ar Q01 och Q02 ¨ar givna konstanter medan ωp1 och ωp2 s¨atts till givna fraktioner av den l¨agsta od¨ampade egenvinkelfrekvensen, se nedan. I denna ber¨akningsuppgift
ska strukturen i fig. 1 analyseras, dels med avseende p˚
a utb¨ojningens storlek vid ytter¨anden, dels med avseende p˚
a risken f¨or utmattningsbrott vid
inf¨astningen.
Kraft-fo
¨rskjutningssamband
Betrakta f¨orst kraft-f¨orskjutningssambandet f¨or en statiskt belastad balk.
Vi anv¨ander klassisk balkteori (Euler-Bernoulli-balk), och eftersom differentialekvationen f¨or en s˚
adan ¨ar linj¨ar kan vi utnyttja superposition. D˚
a g¨aller
3
dels att f¨orskjutningen (utb¨ojningen) beror linj¨art av storleken p˚
a en p˚
alagd
kraft, dels att vi kan ber¨akna f¨orskjutningen f¨or en balk som belastas av flera
krafter som summan av f¨orskjutningarna som skulle f˚
as f¨or varje kraft f¨or sig.
Figur 2: En balk belastad med tv˚
a krafter.
F¨or balken i figuren kan vi allts˚
a skriva
w1 = α11 F1 + α12 F2
w2 = α21 F1 + α22 F2
P˚
a matrisform kan vi skriva ekv. (2)
!
!
α11 α12
w1
=
α21 α22
w2
F1
F2
(2)
!
(3)
Eller
w = AF
(4)
d¨ar A kallas flexibilitetsmatrisen.
Enligt Maxwell-Bettis teorem g¨aller α12 = α21 f¨or det linj¨art elastiska fallet.
Influenstalen αij kan vi ta fram med elementarfall f¨or balkb¨ojning. T.ex. ger
elementarfallet i fig. 3 influenstalet α22 som:
4
α22 =
L32
3EI
(5)
Figur 3: Elementarfall f¨or ber¨akning av influenstalen.
Ibland vill man r¨akna ut krafterna ur f¨orskjutningarna ist¨allet f¨or f¨orskjutningarna
ur krafterna. Om v˚
ar struktur ¨ar l˚
ast f¨or stelkroppsr¨orelser (vilket den ¨ar
eftersom ena ¨anden av balken ¨ar fast insp¨and) kan vi invertera ekv. (4) s˚
a
att vi f˚
ar
A−1 w = F ⇒ F = Kw
(6)
d¨ar K=A−1 ¨ar styvhetsmatrisen, som i m˚
anga fall ist¨allet ber¨aknas med
finita elementmetoden.
Uppgift 1
(1.a) Tag fram alla fyra influenstalen (tecknade med bokst¨aver, inte siffror).
Detta g˚
ar att g¨ora med elementarfallet i fig. 3 tillsammans med MaxwellBettis teorem och superposition, men det ¨ar ocks˚
a till˚
atet att anv¨anda ytterligare elementarfall. (1.b) Implementera influenstalen i Matlab f¨or att kunna ber¨akna flexibilitetsmatrisen A numeriskt med v¨arden fr˚
an din dataupps¨attning. Ber¨akna ¨aven styvhetsmatrisen K genom att invertera A. Matlabkoden ska vara l¨att att f¨olja och utg˚
a direkt fr˚
an din dataupps¨attning, dvs
om du inf¨or mellanresultat ska dessa ber¨aknas direkt i Matlab-koden, inte
p˚
a papper vid sidan av.
5
R¨
orelseekvationerna
Uppgift 2
(2.a) Fril¨agg punktmassorna i strukturen i fig. 1 och teckna deras r¨orelseekvationer
(p˚
a papper). Tag med tyngdkrafterna. Kalla reaktionskrafterna mellan balken
och punktmassorna f¨or R1 och R2 (allts˚
a inte F1 och F2 ; det framg˚
ar varf¨or
i uppgift 3 nedan). Uttryck de tv˚
a r¨orelseekvationerna tillsammans p˚
a matrisform. L˚
at t.ex. punktmassornas accelerationer bilda en kolumnmatris
!
w¨1
¨ =
w
(7)
w¨2
Skriv p˚
a en form s˚
a att det uppst˚
ar en matris, som vi kallar massmatris
M , som inneh˚
aller punktmassornas massor m1 och m2 , och ingenting mer
utom nollor; vi ska komma fram till den ekvation som motsvarar ekv. (15) i
exemplet i slutet av h¨aftet.
Kraft-fo
¨rskjutningssambandet, igen
I uppgift 1 ber¨aknades sambandet mellan krafter och f¨orskjutningar p˚
a balken
f¨or det statiska fallet, vilket gav sambandet F = Kw. Vi skall nu utvidga
detta genom att anta att d˚
a balken r¨or sig tillkommer ett bidrag s˚
a att
˙
R = F + Cw
(8)
d¨ar R ¨ar de totala reaktionskrafterna mellan balken och punktmassorna.1
F¨or det dynamiska fallet har vi allts˚
a:
˙
R = Kw + C w
Eller, p˚
a matrisform:
!
R1
=
R2
k11 k12
k21 k22
!
w1
w2
!
+
(9)
c11 c12
c21 c22
!
w˙ 1
w˙ 2
!
(10)
Det visar sig att antagandet att krafterna p˚
a balken har ett bidrag som a¨r
˙ ger en d¨ampning av strukturen.
proportionellt mot hastigheterna w
Vi antar vidare att d¨ampningsmatrisen kan skrivas
C = ck K + cm M
1
(11)
¨ i ekv. (8). Om balkens massa
Balken antas massl¨
os s˚
a det finns inget beroende av w
tas med m˚
aste vi a
orelseekvationer f¨or sj¨alva balken.
¨ven teckna r¨
6
Detta kallas proportionell d¨ampning och har sin huvudsakliga motivering i
att antagandet har vissa matematiska f¨ordelar, som vi inte utnyttjar h¨ar,
men ocks˚
a i att det ¨ar ett enkelt antagande som underl¨attar best¨ammande
av C. H¨ar m˚
aste vi rimligtvis s¨atta cm =0, ty d¨ampningen uppkommer vid
deformation av balkens material och i v˚
art fall representerar ju M punktmassor p˚
a en massl¨os balk, s˚
a M ¨ar inte relaterad till den materialm¨angd
som deformeras vid strukturens r¨orelse. K ¨ar d¨aremot ett m˚
att p˚
a balkens
storlek. Vi s¨atter allts˚
a
C = ck K
(12)
s˚
a att v˚
art slutresultat f¨or krafterna mellan balken och punktmassorna blir
˙
R = Kw + ck K w
(13)
Uppgift 3
(3.a) Avsluta h¨arledningen av r¨orelseekvationen i uppgift 2 genom att s¨atta
in ekv. (1) och (13) i r¨orelseekvationerna fr˚
an uppgift 2. V˚
ar analys a¨r nu
framme vid den punkt som motsvaras av ekv. (17) i exemplet. (3.b) Ber¨akna
¨
egenvinkelfrekvenserna f¨or det od¨ampade fallet (det ¨ar allts˚
a M w+Kw
=0
som ¨ar den underliggande ekvationen f¨or egenv¨ardesanalysen). Det ¨ar obligatoriskt att anv¨anda SI-enheter.
Uppgift 4
(4.a) Skriv om systemet av r¨orelseekvationer fr˚
an uppgift 3 som ett system av
fyra f¨orsta ordningens differentialekvationer. Vi ska allts˚
a genomf¨ora det som
motsvarar stegen fr˚
an ekv. (17) till (25) via ekv. (22) i exemplet. Vi g¨or detta
eftersom Matlabs l¨osare bara klarar f¨orsta ordningens differentialekvationer.
Uppgift 5
(5.a) Implementera ditt system av r¨orelseekvationer, som efter omskrivningen
i uppgift 4 ¨ar ett system av fyra f¨orsta ordningens ordin¨ara differentialekvationer. Anv¨and n˚
agon av Matlabs numeriska l¨osare t.ex. ode45 som i exemplet. Det ¨ar obligatoriskt att anv¨anda SI-enheter. (5.b) G¨or en ber¨akning d¨ar
Q01 och Q02 satts till noll och ber¨akningstiden till 5 sek. och plotta resultatet
(alla fyra funktionerna i din ber¨akning, utb¨ojningar och hastigheter). Balken
sl¨apps fr˚
an vila i horisontellt l¨age, s˚
a alla funktionerna s¨atts till noll som begynnelsevillkor. Bed¨om om ber¨akningstiden a¨r tillr¨acklig f¨or att ber¨akningen
7
ska n˚
a fortvarighet (steady-state), dvs om balkens sv¨angning hunnit d¨ampas
ut till v¨aldigt sm˚
a r¨orelser i f¨orh˚
allande till st¨orsta utb¨ojning. Om inte, ¨oka
ber¨akningstiden; 10 sek. ska vara tillr¨ackligt f¨or alla dataupps¨attningar.
Uppgift 6
(6.a) S¨att Q01 och Q02 till v¨arden enligt ditt datablad. Unders¨ok om f¨orskjutningen
w2 i balkens ¨ande vid fortvarighet ¨overskrider till˚
atet v¨arde enligt datablad
d˚
a ωp1 och ωp2 s¨atts lika och till 50% respektive 95% av l¨agsta od¨ampade
egenvinkelfrekvensen (denna ber¨aknades i uppgift 3).
Uppgift 7
(7.a) Ber¨akna hur krafterna R1 och R2 mellan balken och punktmassorna varierar i tiden genom att s¨atta in dina utr¨aknade v¨arden p˚
a w och
˙ i ekv. (13); detta g¨ors som en behandling av utdata fr˚
w
an differentialekvationsl¨osaren, se exemplet i slutet p˚
a h¨aftet. (7.b) H¨arled ekvationer f¨or
att ber¨akna snittmomentet i balken vid inf¨astningen i v¨aggen ur krafterna,
genom att teckna j¨amviktsekvationer f¨or balken frilagd fr˚
an punktmassorna
och v¨aggen (det blir j¨amviktsekvationer och inte r¨orelseekvationer eftersom
balken antas massl¨os). (7.c ) Anv¨and momentet f¨or att med element¨ar balkteori ber¨akna dragsp¨anningen p˚
a grund av b¨ojning vid inf¨astningen i balkens
¨overkant. (7.d ) Implementera sp¨anningsber¨akningen i din Matlab-kod och
ber¨akna och plotta sp¨anningen f¨or fallet d˚
a ωp1 och ωp2 s¨atts lika och till noll
respektive till 50% av l¨agsta od¨ampade egenvinkelfrekvensen.
Uppgift 8
Balkens inf¨astning i v¨aggen ¨ar utformad enligt fig. 4. (8.a) Tag fram reduktionsfaktorn Kf f¨or k¨alverkan. (8.b) Rita ett Haigh-diagram. Reduktion av
utmattningsdata g¨ors bara med avseende p˚
a k¨alverkan, ¨ovriga faktorer s¨atts
till ett. (8.c) Ber¨akna s¨akerheten mot utmattning avseende ¨andring av amplitudsp¨anningen d˚
a ωp1 och ωp2 s¨atts lika och till 50% respektive 95% av l¨agsta
od¨ampade egenvinkelfrekvensen. Sp¨anningarna vid fortvarighet anv¨ands vid
ber¨akning av s¨akerhetsfaktorerna, dvs ta inte med eventuella transienter i
b¨orjan av den numeriska l¨osningen. Det kan intr¨affa att du f˚
ar en s¨akerhetsfaktor
som ¨ar mindre ¨an ett f¨or fallet d˚
a ωp1 och ωp2 s¨atts till 95% av l¨agsta
od¨ampade egenvinkelfrekvensen. Observera att p˚
a grund av den konstanta lasten fr˚
an gravitationen kommer sp¨anningen variera kring en positiv
(drag) mittsp¨anning p˚
a ovansidan och en negativ (tryck) mittsp¨anning p˚
a
undersidan. Det a¨r h˚
alk¨alen med positiv mittsp¨anning som ska analyseras
8
Figur 4: Balkens inf¨astning vid v¨aggen.
eftersom en negativ mittsp¨anning (tryck) inte anses f¨ors¨amra utmattningsh˚
allfastheten. (8.d) Ber¨akna s¨akerheten mot att till beloppet st¨orsta nominella sp¨anning ¨overskrider str¨ackgr¨ansen, dvs σs /σmax .
9
Appendix: exempel
Figur 5: System vars r¨orelse analyseras i exemplet.
Vi ska st¨alla upp r¨orelseekvationerna f¨or systemet i fig. 5 och l¨osa dessa
numeriskt med hj¨alp av Matlab.
R¨
orelseekvationer
Fril¨agg massorna enligt fig. 6 och teckna r¨orelseekvationer i horisontell led:
−R1 + R2 = m1 w
¨1
−R3 + P = m2 w
¨2
(14)
eller, p˚
a matrisform:
m1 0
0 m2
w¨1
w¨2
!
10
=
−R1 + R2
−R3 + P
(15)
Figur 6: Fril¨aggning av massorna f¨or systemet i exemplet.
Vi antar att anordningarna som kopplar ihop massorna med varandra och
med v¨aggen ger krafter enligt
R1 = k1 w1 + c1 w˙ 1
(16)
R2 = R3 = k2 (w2 − w1 ) + c2 (w˙ 2 − w˙ 1 )
dvs respektive kraft f˚
ar ett bidrag beroende p˚
a fj¨aderns deformation och
ett bidrag beroende p˚
a d¨amparens deformationshastighet. Vi antar att kopplingsanordningarna ¨ar massl¨osa s˚
a att R2 =R3 .
Ekv. (15) och (16) ger
m1 0
0 m2
! w¨1
c1 + c2 −c2
+
−c2
c2
w¨2
! w˙ 1
k1 + k2 −k2
+
−k2
k2
w˙ 2
! w1
0
=
P (t)
w2
(17)
Eller, p˚
a symbolisk form
¨ + Cw
˙ + Kw = P
Mw
(18)
R¨
orelseekvationerna som ett f¨
orsta ordningens system
R¨orelseekvationerna enligt ekv. (17) ¨ar ett system av andra ordningens ordin¨ara differentialekvationer. F¨or att kunna l¨osa systemet numeriskt i Matlab
11
m˚
aste vi skriva om det som ett system av f¨orsta ordningens differentialekvationer. Detta ˚
astadkommer vi genom att inf¨ora massornas hastigheter som
egna storheter enligt definitionerna
w˙ 1 = v1
w˙ 2 = v2
(19)
w¨1 = v˙ 1
w¨2 = v˙ 2
(20)
Deriverar vi ekv. (19) f˚
ar vi
Ekv. (19) och (20) i (17) ger
m1 0
0 m2
! v˙ 1
c1 + c2 −c2
+
−c2
c2
v˙ 2
! v1
k1 + k2 −k2
+
−k2
k2
v2
w1
w2
!
=
0
P (t)
(21)
Nu har vi inte l¨angre n˚
agra andraderivator i systemet av differentialekvationer, men ekv. (21) g˚
ar inte att l¨osa eftersom det bara ¨ar tv˚
a ekvationer f¨or
fyra obekanta funktioner w1 (t), w2 (t), v1 (t) och v2 (t). Vi l¨oser detta problem
genom att notera att ekv. (19) ¨ar tv˚
a differentialekvationer f¨or w˙ 1 och w˙ 1 .
Genom att kombinera ekv. (21) och ekv. (19) kan vi skriva:

1
 0

 0
0
  w˙  
1
0 0
0
0
0
−1
0

 

w
˙
1 0
0   2 
0
0
0
−1
 +
0 m1 0   v˙ 1   k1 + k2 −k2 c1 + c2 −c2
0 0 m2
−k2
k2
−c2
c2
v˙ 2
w  
1
0

  w2 

0
 
=
  v1 
  0
P (t)
v2
(22)
Implementering i Matlab
F¨or att kunna implementera v˚
art system av ordin¨ara differentialekvationer,
ekv. (22), i Matlab m˚
aste vi ¨andra notationen n˚
agot. Vi f˚
ar inte ha s¨okta
funktioner med olika namn, utan vi m˚
aste ge dem samma namn och skilja dem ˚
at enbart med ett index som anger vilken position i en gemensam
kolumnmatris som respektive funktion har. Vi inf¨or d¨arf¨or:
   
y1
w1
y  w 
 2  2
(23)
 = 
 y3   v1 
y4
v2
12




  
w˙ 1
ydot1
 ydot   w˙ 
2
 2

= 

 ydot3   v˙ 1 
v˙ 2
ydot4

s˚
a att ekv. (22) kan skrivas

(24)
y  
1
0
0
y  

1
0
2

  +
c2 /m1   y3  
0
−c2 /m2
P (t)/m2
y4
(25)
d¨ar vi ¨aven multiplicerat med inversen av matrisen l¨angst till v¨anster i ekv.
(22) i v¨ansterledet, vilket i det h¨ar fallet blir samma sak som att dividera de
tv˚
a sista ekvationerna med m1 respektive m2 . (Vissa tecken har ¨andrats d˚
a
termer flyttats till andra sidan likhetstecknet.)
F¨or att l¨osa elv. (25) numeriskt i Matlab skapar vi f¨orst filen myode.m:
 
ydot1
0
 ydot  
0
2


=

−(k
+
k2 )/m1
 ydot3 
1
k2 /m2
ydot4
0
1
0
0
k2 /m1
−(c1 + c2 )/m1
−k2 /m2
c2 /m2
%This is the file myode.m which, by calculating
%the column matrix of derivatives of the functions
%to be solved for by Matlab solver ode45, defines
%the governing system of ordinary differential equations.
function ydot=myode(t,y,m1,m2,k1,k2,c1,c2,P0,wp)
P=P0*sin(wp*t);
ydot(1)=y(3);
ydot(2)=y(4);
ydot(3)=-(k1+k2)*y(1)/m1+k2*y(2)/m1-(c1+c2)*y(3)/m1+c2*y(4)/m1;
ydot(4)=k2*y(1)/m2-k2*y(2)/m2+c2*y(3)/m2-c2*y(4)/m2+P/m2;
ydot=transpose(ydot);
I denna fil ber¨aknas en kolumnmatris som ger tidsderivatorna av alla s¨okta
funktioner i v˚
art system av ordin¨ara differentialekvationer, ekv. (25). Transponeringen p˚
a sista raden a¨r f¨or att man f˚
ar en radmatris n¨ar matrisen definieras
element f¨or element som i filen, medan Matlabs l¨osare kr¨aver att det ska vara
en kolumnmatris.
Filen myode.m anropas inte direkt; vi anropar Matlabfunktionen ode45
som i sin tur kommer utnyttja myode.m. F¨or detta ¨andam˚
al skapar vi filen
solveode.m nedan:
13




%This is the file solveode.m where the call to ode45
%is made. The Matlab solver ode45 will calculate a
%numerical solution to the system defined in myode.m
%and return this solution to the matrices times and yout.
%*****data provided*****************************
m1=2.;
m2=2.;
k1=6.;
k2=4.;
c1=1.;
c2=1.;
P0=1.;
%*****egigenvalue calculation****************
%*****details omitted here*******************
we1=1.4142;
%******calculation of excitation frequency***
wp=0.9*we1;
%*****initial conditions*********************
yinit(1)=0.;
yinit(2)=0.;
yinit(3)=0.;
yinit(4)=0.;
%*****timespan*******************************
tstart=0;
tend=50.;
%******call to numerical solver**************
options=odeset(’RelTol’,1e-6,’AbsTol’,1e-6);
[times,yout]=ode45(@myode,[tstart,tend],yinit,options,...
m1,m2,k1,k2,c1,c2,P0,wp);
%*****plotting of positions and velocities***
for i=1:4
subplot(2,2,i);
plot(times,yout(:,i));
end
14
Kommentarer till Matlab-koden
Studera speciellt sj¨alva anropet till differentialekvationsl¨osaren:
options=odeset(’RelTol’,1e-6,’AbsTol’,1e-6);
[times,yout]=ode45(@myode,[tstart,tend],yinit,options,...
m1,m2,k1,k2,c1,c2,P0,wp);
H¨ar ¨ar ode45 namnet p˚
a en av flera l¨osare man kan v¨alja mellan i Matlab,
@myode anger att systemet som definieras av f¨orstaderivatorna i myode.m ska
l¨osas, [tstart,tend] ¨ar tidsintervallet d¨ar l¨osningen ska ber¨aknas, yinit
anger initialvillkor (startvillkor), options skapas av odeset p˚
a raden innan
och a¨ndrar vissa parametrar i l¨osarens beteende, ... anger att n¨asta rad i
filen ¨ar en forts¨attning p˚
a raden som p˚
ab¨orjats, och m1,m2,k1,k2,c1,c2,P0,wp
¨ar parametrar vi vill skicka vidare till myode.m via ode45. I exemplet anv¨ands
options f¨or att ¨oka noggrannheten i l¨osningen. Notera att det inte g˚
ar att
hoppa ¨over argumentet options om man vill forts¨atta med att ange parametrar som ska vidarebefordras till myode.m p˚
a det s¨att som gjorts i exemplet.
V¨ansterledet i anropet anger att utdata ska placeras i times och yout
d¨ar times ¨ar en kolumnmatris med de tidpunkter d¨ar ode45 valt att ber¨akna
den numeriska l¨osningen och yout ¨ar en matris med en kolumn f¨or varje s¨okt
funktion, i exemplet fyra kolumner.
Sist i filen finns instruktioner f¨or att plotta de ber¨aknade funktionerna:
for i=1:4
subplot(2,2,i);
plot(times,yout(:,i));
end
H¨ar betyder subplot(2,2,i) att f¨onstret d¨ar funktionerna plottas ska delas
upp i tv˚
a rader och tv˚
a kolumner till fyra delar, och vi plottar f¨or tillf¨allet i
f¨onsterdel nummer i. plot(times,yout(:,i)) plottar kolumn nummer i av
yout mot de tidpunkter som anges i kolumnmatrisen times.
Skriver vi solveode vid Matplabprompten utf¨ors instruktionerna, och
vi f˚
ar ber¨akningsresultatet enligt fig. (7). Skriv ocks˚
a times och yout vid
Matlabprompten (efter ber¨akningen) f¨or att f¨orst˚
a vilka data vi f˚
ar ut av
ber¨akningen.
Ber¨
akning av krafterna
Reaktionskrafterna mellan massorna och kopplingsanordningarna ber¨aknas
ur f¨orskjutningar och r¨orelser. Ekv. (19) och (23) i (16) ger
15
0.5
1
0.5
0
0
−0.5
−0.5
0
10
20
30
40
−1
50
1
0
10
20
30
40
50
0
10
20
30
40
50
1
0.5
0.5
0
0
−0.5
−0.5
−1
−1
0
10
20
30
40
50
−1.5
Figur 7: Resultat av ber¨akningen.
R1 = k1 y1 + c1 y3
R2 = R3 = k2 (y2 − y1 ) + c2 (y4 − y3 )
(26)
Vi kan r¨akna ut krafterna genom att l¨agga till en ber¨akning enligt ekv. (26)
i filen solveode.m, givetvis efter anropet till ode45:
n=size(times);
for i=1:n
R1(i)=k1*yout(i,1)+c1*yout(i,3);
R2(i)=k2*(yout(i,2)-yout(i,1))+c2*(yout(i,4)-yout(i,3));
end
Med instruktionen n=size(times) f˚
ar vi veta hur m˚
anga element kolumnmatrisen times har, dvs. vid hur m˚
anga tidpunkter ode45 ber¨aknat l¨osningen
numeriskt. I yout ges sedan motsvarande funktionsv¨arden s˚
a att varje kolumn ger en av de s¨okta funktionerna. Exempelvis ¨ar yout(i,4) den h¨ogra
massans hastighet vid den tidpunkt som ges av times(i).
16