HH/ITN/BN Mekanik-Dynamik och Mathematica 1 Något om Mekanik-Dynamik och Mathematica, läsanvisningar till Christer Nybergs bok i Mekanik Bertil Nilsson 2015-01-01 Jump Bungy Solve the initial value problem using the following rubber band model. Frb yOft k H y t L c y' t H y t L 0 H y t L NDSolve m y'' t m g Frb , y 0 . m 75, g 9.81, H 50, y t , t, 0, 15 ; 0 ; 0 H, y' 0 L 30, k 0 125, c A picture illuminates the situation Plot Evaluate PlotStyle AxesLabel y t . yOft, D y t . yOft, t , t, 0, 15 , Red, Blue , AxesStyle Directive Bold, 12 , "t s ", "y t m ,y t m s " yt m,yt ms 40 20 2 20 4 6 8 10 12 14 t s 50 , 2 Mekanik-Dynamik och Mathematica HH/ITN/BN ť Dynamik Till skillnad från statik som behandlar kroppar i vila, handlar dynamik om kroppar i rörelse. Dynamiken delas dessutom in i - kinematik som behandlar kroppars rörelse utan påverkan av krafter. Här handlar det om att den bana som en kropp rör sig i är känd och man söker ofta de krafter som detta ger upphov till, exempelvis tåg eller karusell på Liseberg. Enligt Newtons andra lag nedan så upplever en resenär i en karusell tröghetskraften massa acceleration. - kinetik handlar om den rörelse eller bana som en kropp tvingas till utgående från givna yttre krafter, exempelvis en fotboll. Krafter uppkommer på grund av påverkan från omgivningen, exempelvis tyngdkraft och kontaktkrafter. Centralt i dynamiken är Newtons tre lagar. Dessa talar om hur en kraft påverkar en kropp: N1: Tröghetslagen. En kropp förblir i sitt tillstånd av vila eller likformig rörelse om den inte av verkande krafter tvingas att ändra detta tillstånd. N2: Accelerationslagen. En kropp som påverkas av kraften F får en acceleration a sådan att F ma, där konstanten m är kroppens (tröga) massa. Egentligen sa Newton något annat, se vidare under avsnittet kinetik. N3: Lagen om verkan och motverkan. Två kroppars ömsesidiga verkningar på varandra är alltid lika stora och riktade åt motsatt håll. Kinematik I kinematik läggs den matematiska grunden för dynamiken, alltså även kinetiken. Notera att hela dynamiken vilar på att man kan frilägga, det vill säga statiken måste vara väl inhämtad. Samband mellan läge, hastighet och acceleration utreds och är mycket centralt. Arbeta som vanligt i ett universellt fastspikat koordinatsystem och håll fast vid det! Vanligtvis brukar vi göra studiet i vårt vanliga koordinatsystem med beteckningarna x, y och z. Arbeta hela tiden i SI-enheter! Enheten radianer måste användas i dynamik ty deriveringsregler gäller inte annars. Grader kan möjligtvis tillåtas så länge vinkeln inte ska deriveras. Läget för tyngdpunkten av en kropp eller punkt som funktion av tiden beskrivs sedan av sin ortsvektor t x t , y t , z t . Koordinaten, "läget", "vägen" skall betraktas som en vektor med både längd och riktning! Derivator och integraler av en vektor eller matris faller naturligt ut på komponenterna, t x t , y t , z t respektive t t x t t, y t t, z t t . Dimensionsanalys av ekvationer och svar är viktigt. Tillskott av läget x kallas x. Tillskottet av läge per tidsenhet på derivata har vi att gränsvärdet lim x t t 0 fysik, därför har man infört beteckningen x hastighet per tidsenhet x . t x t x , t x t kallas för medelhastighet under tidsintervallet t. Enligt definition som kallas (momentan)hastighet vid tiden t. Tidsderivator är mycket vanliga i "x-prick" för tidsderivata. På samma sätt har vi accelerationen som ett tillskott av .. För denna inför vi beteckningen x x , t "x-prickprick". Så x . t Hastighet vektor är tidsderivatan av läget vektor x .. Acceleration vektor är tidsderivatan av hastigheten vektor x x t 2 x t2 . I många fall är man inte primärt intresserad av x t och x t utan hastigheten som funktion av läget x x . Då är det mycket vanligt att .. göra omskrivningen x x t KR x x x t x x . x Detta brukar ofta leda direkt till en separabel (ODE). Då vi betraktar x t säger vi att vi arbetar i tidsplanet och i fasplanet när x x . Lös inte problem i dynamik med färdiga formler av den typ som härleds i läroböcker och som förekommer i alla formelsamlingar för gymnasiet. Det är nämligen mer arbete att utreda om formeln gäller med hänsyn till begynnelsevärden och övriga förutsättningar än att lösa rörelseekvationerna själv. Många formler förutsätter accelerationen konstant. Använd inte beteckningarna s, v och a för varierande läge, hastighet och acceleration - även om - eller i synnerhet inte - om du från gymnasiet är inarbetad på dessa beteckningar! Du får bara allt svårare för att frigöra dig från "färdiga-formler-tänkandet" och kommer inte vidare!! Använd .. istället x, x och x tillsammans med dina kunskaper från kurs i ODE!! Läs i "Något om (ODE) och Mathematica". De olika arbetsmoment som kinematik och kinetik erbjuder kan sammanfattas i följande schema. För att förflytta sig mellan de önskade storheterna läge, hastighet och acceleration är det derivation som gäller i ena rikningen och integration i den andra. Notera att lösningen till ett begynnelsevärdesproblem (BVP) tar oss i allmänhet hela vägen "upp" till x t . Hastighet fås sedan genom .. derivering. På grund av Newtons accelerationslag, m , är modeller i kinetik naturligt förmulerade som en andra ordningens (ODE), men kinematik är lika ofta formulerat som en första ordningens (ODE). För alla dessa moment i schemat finns Mathematica till mycket god hjälp. HH/ITN/BN Mekanik-Dynamik och Mathematica Läge, x x t Begynnnelsevärdesproblem .. x g x, x ODE BVP x 0 x0 BV x 0 x0 Mathematica DSolve NDSolve x x t Hastighet, x x t 3 .. x t .. Acceleration, x f t Fig 1. Lösningsmetoder Även i dynamik behöver vi våra nyttiga funktioner från statiken. Först storleken av vektorn . Tyvärr kan vi inte använda vanliga | som i linjär algebra när vi definierar funktionen, eftersom den är upptagen till annat i Mathematica, utan väljer några snarlika som skrivs l| till vänster (left) och r| till höger (right) om argumentet. Det väsentliga i funktionskroppen är "Matematicapornografi" för att förenkla symboliska uttryck där alla variabler är positiva. : PowerExpandSimplify . , Thread Cases , Symbol, Infinity . , resten är 0 Ständigt behöver man bestämma en enhetsvektor i samma riktning som en given vektor : liksom beräkning av enhetsvektor i riktning från punkt A till punkt B A ,B : B A samt slutligen projektion av en vektor på en annan vektor . proj : . Om man inte uppskattar svåra ackord på tangentbordet, kommer här samma funktioner i mer klassisk form. Använd vilka versioner du vill eller andra namn som du trivs med. Speciellt kan det som alternativ vara läge att använda de i Mathematica inbyggda . funktionerna Norm[ ] som motsvarar och Normalize[ ] som motsvarar . längd : ea : eAB A , B : B A proj , : proj I utdata använder Mathematica gärna funktionerna sec Α 1 cos Α , Sec Α och cosec Α 1 sin Α , Csc Α . Dessa är en kvarleva från navigering på de gamla segelfartygens tid och anses numera, åtminstone i Europe, som lite exotiska , och ingår därför inte i den svenska skolundervisningen. Det samma får nog gälla för cot Α 1 tan Α , Cot Α . Vill man inte se dessa är det lämpligt att också aktivera följande. $PrePrint . Csc z Sec z Cot z 1 1 1 Defer Sin z , Defer Cos z , Defer Tan z &; sid 123-129: Läs igenom. CN bryter naturligtvis mot den gula rutan ovan, även om han punktvis skådar ljuset! Att alltid arbeta i ett koordinatsystem är viktigt. Dimensionsanalys av ekvationer och svar är viktigt. I många fall är man inte primärt intresserad av x t och x t utan hastigheten som funktion av läget x x . Då är det mycket vanligt att ta hjälp av kedjeregeln och göra omskrivningen .. x x t KR x x x t x x . x Detta brukar ofta leda direkt till en separabel (ODE). sid 130-131: Lite om självklara geometriska insikter, typ att "arean under v–t-kurvan", det vill säga integralen där area räknas med tecken (CN skiljer sig inte från de flesta andra mekanikböcker när det gäller att uttrycka sig korrekt), är den tillryggalagda vägen. sid 132-136: Att använda styckvis definierade funktioner är en tidig IQ-test av CN och lägger en rejäl teknisk dimridå över det som ska förmedlas! Vi börjar med ett enkelt exempel med en kropp som rör sig längs x-axeln (sk rätlinjig rörelse) som funktion av tiden x t sin t . Först definerar vi funktionen i Mathematica sedan ritar vi en graf med kroppens läge, hastighet och acceleration som funktion av tiden. Notera att Mathematica förstår den vanliga sparv -notationen för derivata men inte prickar för tidsderivata. 4 Mekanik-Dynamik och Mathematica x t : Sin t Plot x t , x ' t , x '' t PlotStyle HH/ITN/BN , t, 0, 2 Π , Red, Blue, Green , AxesLabel .. "t", "x,x, x " .. x,x,x 1.0 0.5 t 1 2 3 4 5 6 0.5 1.0 Enligt diagramläsning ska "arean under v–t-kurvan", dvs integralen, vara den tillryggalagda vägen. Vi gör en inspektion då t Plotx ' t , t, 0, 2 Π , PlotStyle Blue, Filling Axis, AxesLabel 2Π. "t", "x" Enligt den röda kurvan är vi tillbaka där vi startade vid t 0, så läget är detsamma. Däremot har vi varit ute på en liten semestertur som kan avläsas på trippmätaren vid hemkomsten. Vi kollar båda 2Π x ' t , Abs x ' t t 0 0, 4 Exempel 6.1 sid 132: Typexempel på då det är bra att göra en egendefinierad funktion. Om funktionen är styckvis definierad finns Piecewise pw , eller de lite mer programmeringsmässiga If och Which. x t : t2 t 20 t 100 t 10 10 Plot x t , t, 0, 30 , PlotStyle Red, AxesLabel "t", "x t " xt 500 400 300 200 100 5 10 15 20 25 30 t Plot x ' t , t, 0, 30 , PlotStyle Blue, PlotRange All, AxesLabel vt 20 15 10 5 t 5 10 15 20 25 30 Arean under v t -kurvan. Stämmer bra med värdet som avläses i första figuren. 30 0 500 x' t t "t", "v t " HH/ITN/BN Mekanik-Dynamik och Mathematica Plot x '' t , t, 0, 30 , PlotStyle Brown, AxesLabel 5 "t", "a t " at 2.0 1.5 1.0 0.5 t 5 10 15 20 25 30 Exempel 6.2 sid 133: Väsentligen samma som Exempel 6.1, fast tvärtom ;-). Samtidigt tar vi chansen att visa på hur enkelt det är att spara papper. Först acceleration, sedan hastighet med start från vila, slutligen läget med start i origo. Allt enligt integrationsvägen i fig 1. Efter bilderna introducerar vi DSolve och NDSolve för att kunna gå (BVP)-vägen i fig 1. Därefter använder vi oss av detta i en repris av detta exempel. a t 10 t 2 t t 10 ;v t 10 t a Τ Τ; x t v Ω 0 Ω; 0 Plot a t , t, 0, 20 , PlotStyle Plot v t , t, 0, 60 , PlotStyle Plot x t , t, 0, 60 , PlotStyle Brown, AxesLabel "t", "a t " , Blue, AxesLabel "t", "v t " , Red, PlotRange All, AxesLabel "t", "x t " xt vt at 100 10 8 6 4 2 1000 80 5 2 10 15 20 800 , 60 , 600 40 400 t 20 200 10 20 30 40 50 60 t 10 20 30 40 50 60 t I Mathematica används funktionen DSolve för att lösa en stor klass av differentialekvationer, allt från enkla separabla och linjära av godtycklig ordning till mycket komplicerade olinjära. Vi börjar med en enkel linjär första ordningens differentialekvation y' x y x 0 och dess lösning. DSolve y ' x yx y x 0, y x , x x c1 DSolve löser även system av differentialekvationer och är mycket lättanvänd. Strängt taget handlar det om att skriva av rätt! Lägg märke till att Mathematica förstår den vanliga nomenklaturen när vi menar derivata. Notera de nödvändiga []-parenteserna eftersom y x ska vara en funktion av x! För övrigt kan man naturligtvis använda vilka namn man vill. Observera dubbla likhetstecken eftersom det är en ekvation! Det är inte bara namnet som antyder släktskap med Solve, utan även hantering av indata och resultat. Som vanligt gäller att när man väl förstått filosofin bakom Mathematica är det mesta självklart! DSolve y ' x y x sin x cos x 2 DSolve1 y x Sin x , y x , x 1 x c1 y x x2 y ' x 1 c1 x2 1 x2 2xy x 1 tan 1 1 x2 ArcTan x , y x , x x 2 2 Man kan givetvis ta med begynnelsevärden för att få konstanterna bestämda. Dessa paketeras då tillsammans med differentialekvationen i en lista. Tänk på att även begynnelsevärdena ska anges som ekvationer, det vill säga med två likhetstecken. Så begynnelsevärdesproblemet BVP y' x y0 y x 2 sin x ODE BV är bara att muppa rakt in i Mathematica, sedan piggar vi upp oss med en bild. yAvx DSolve 1 y x 5 2 x y' x sin x y x cos x Sin x , y 0 2 ,y x ,x Simplify 6 Mekanik-Dynamik och Mathematica Plot y x . yAvx, x, 0, 5 , PlotStyle Red, AxesLabel HH/ITN/BN "x", "y x " yx 2.0 1.5 1.0 0.5 1 2 3 4 5 x 0.5 2t I nästa exempel noterar Mathematica naturligtvis att både partikulärlösningen därefter. DSolvex '' t 1 2t x t 2 4 x' t 2 c2 t 2t 4x t och t , x t , t 2t ingår i homogena lösningen och korrigerar ansatsen till Simplify t2 2 c1 Då vi inte kan finna lösningen analytiskt, vilket naturligtvis är långt ifrån ovanligt ute i verkliga livet, finns funktionen NDSolve till vår hjälp för att göra en ren numerisk lösning. Utdata från denna är en InterpolatingFunction som kan verka lite märkvärdig innan man blivit vän med den. Den fungerar dock som vilken annan funktion som helst. För övrigt är den ett kraftfullt redskap om man vill göra interpolation i t.ex. mätdata. Som exempel kör vi en repris på begynnelsevärdesproblemet ovan. Det enda som skiljer i menyn jämfört med DSolve är att man, likt Plot, måste ange i vilket intervall man vill att spektaklet ska utspela sig. NyAvx y x NDSolve y' x y x Sin x , y 0 Domain: 0. 5. InterpolatingFunction Plot y x 2 , y x , x, 0, 5 Output: scalar . NyAvx, x, 0, 5 , PlotStyle x Red, AxesLabel "x", "y x " yx 2.0 1.5 1.0 0.5 x 1 2 3 4 5 0.5 Exempel 6.2 sid 133 igen: Nu tar vi exemplet en gång till. Förra gången integrerade vi allt själva för hand enligt vänstra kolumnen i fig 1, nu tar vi högra kolumnen istället och formulerar det som ett begynnelsevärdesproblem (BVP) och låter DSolve, eller NDSolve vid svårare (ODE), göra hela jobbet. Smidigt! Accelerationen under resans gång varierar enligt en styckvis definierad funktion, därför passar Piecewise pw bra, eller den lite mer programmeringsmässiga If. xAvt DSolvex '' t x 0 x t 5 t2 t2 120 t 10 t 2 t 0, x ' 0 10 , 10 0, x t , t t 10 600 True Restid tillbaka till stillastående? T t FindRoot x ' t 0 . D xAvt, t , t, 50 60. Lite grafer kanske? Först läget, sedan hastighet och slutligen accelerationen, i detta fall reproduktion av indata. HH/ITN/BN Mekanik-Dynamik och Mathematica Plot x t . xAvt, t, 0, t . T , PlotStyle 7 Darker Green , AxesLabel "t s ", "x m " x m 3000 2500 2000 1500 1000 500 t s 10 20 30 40 PlotEvaluate x ' t Filling 50 60 . D xAvt, t Axis, FillingStyle , t, 0, t . T , PlotStyle LightBlue, AxesLabel "t Blue, s ", "x m s " x ms 100 80 60 40 20 10 20 30 40 PlotEvaluate x '' t Filling 50 t s 60 . D xAvt, t, 2 Axis, FillingStyle , t, 0, t . T , PlotRange LightRed, AxesLabel "t s ", "x 5, 15 , PlotStyle Red, m s2 " x m s2 15 10 5 10 20 30 40 50 t s 60 5 Läge och hastighet då den precis ska börja bromsa? Jämför fig ovan! xAvt, D xAvt, t .t 9.99 x 9.99 499.001 99.9 x 9.99 Total körsträcka. Jämför fig ovan! xAvt . T x 60. 3000. Exempel 6.3 sid 135: Ännu en höghöjdsträning på (ODE)-lösaren i Mathematica. Accelerationen är styckvis definierad. Under accelerationsfasen är den konstant, så v a0 t v t a0 30 . 4 Placera geparden i origo vid tidens början och lös begynnelsevärde- sproblemet .. x BVP 30 4 t 4 0 t 4 x 0 x0 0 0 ODE BV : Start i origo då t 0. BV : Start från stillastående, då t 0. 30 xAvt DSolvex '' t Ift , 0, x 0 4, 4 15 t2 x t t 4 30 t 2 4 True Denna gång ritar vi i samma diagram för omväxlings skull 0, x ' 0 0, x t , t First 8 Mekanik-Dynamik och Mathematica HH/ITN/BN Plot Evaluate x t . xAvt, x ' t . D xAvt, t , t, 0, 10 , PlotStyle Red, Blue , PlotRange All, AxesLabel "t", "x t ,v t " x t ,v t 250 200 150 100 50 t 2 4 6 8 10 Slutligen svaret på frågan xAvt . t x 10 10 240 Notera fördelen med att behålla "Rule" formatet! Vi får ett självdokumenterande svar x T svar. Exempel 6.4 sid 136: Vi tar polisjakten helt numeriskt. Först accelerationsfunktionen för polisen. Låt polisen vara i origo vid tidens början. xPAvt DSolve xP '' t 2.5, xP 0 0, xP ' 0 0 , xP t , t First 1.25 t2 xP t Polisens restid till 50 m/s. Vi bekänner oss bara till SI-enheter!! t1 t . Solve xP ' t 50 . D xPAvt, t , t First 20. Så accelerationsfunktionen och bilarnas lägen aP t 2.5 t 0 t t1 ; t1 Nu bilarnas läge som funktion av tiden. xPAvt xP t xAAvt xA t DSolve xP '' t aP t , xP 0 0, xP ' 0 0 , xP t , t First 1.25 t2 t 20 50. t 500. True DSolve xA '' t 0, xA 0 0, xA ' 0 40 , xA t , t First 40 t Plot Evaluate xP t . xPAvt, xA t . xAAvt , t, 0, 60 2500 2000 1500 1000 500 10 20 30 40 50 60 Så restid och sträcka till ingripande. Solve klarar inte att hantera styckvis definierade funktioner. xPAvt . FindRoot xP t xP 50. xA t . xPAvt . xAAvt, t, 60 2000. sid 137-142: Läs igenom. Är olika typfall som slentrianmässigt brukar behandlas i gamla mekanikböcker. Vi tar I på sid 137 med Mathematica. HH/ITN/BN Mekanik-Dynamik och Mathematica DSolve x '' t 1 k t2 x t k, x 0 2 t x0 v0 , x ' 0 9 x0 , x t , t 2 v0 2 Vi har löst ett begynnelsevärdeproblem. Utan begynnelsevärden kommer det som vanligt ut lika många konstanter ci som vi har ordning på differentialekvationen. DSolve x '' t k, x t , t k t2 x t c2 t c1 2 Exempel 6.5 sid 138: Placera origo vid inbromsningens början då vi även startar klockan! xAvt DSolve 1 x t 2 x '' t r, x ' 0 v0 , x 0 0 ,x t ,t First r t2 2 t v0 Hur mycket är klockan då tåget stannar? T Solve x ' t v0 t 0 . D xAvt, t , t First r Slutligen bromssträckan, det vill säga resvägen till stillastående. xAvt . T x v0 v20 r 2r Exempel: Inte sällan ingår en av konstanterna i differentialekvationen som okänd. Det är då tänkt att bestämma denna ur ett (extra) randvillkor (RV). Som bekant kan vi inte sätta fler (BV)+(RV) än vi har "sparvar" i (ODE):n. Detta kan vi lösa med ett litet trick: Låt DSolve betrakta den okända konstanten som en funktion vars derivata då är noll! Vi får då ett system med två (ODE) och kan nu klämma in ytterligare ett (RV) eller (BV). Som typexempel väljer vi att söka den konstanta accelerationen a under inbromsning från 30 m/s till 10 m/s på 50 m. .. Lösningsförslag: Vi börjar med den populära omskrivningen x DSolve v x v' x a x , a' x v 0 30, v 50 10 , v x ,a x ,x a x 8, v x 2 225 x x x a och låter hastigheten x v x . Då får vi 0, 4 x Normalt skulle vi gjort separation av variablerna, integration samt efterföljande ekvationslösning för bestämning av a. 10 ekv 50 x x 30 400 a x 0 50 a Solve ekv a 8 Med samma svar som ovan. Man kan notera att varianten med DSolve är lite rikare eftersom vi som biprodukt även får v x . Exempel 6.7 sid 140: Typiskt begynnelsevärdesproblem igen. xAvt DSolve x '' t k x' t , x' 0 v0 , x 0 0 ,x t ,t Simplify k t v0 1 x t k Exempel 6.9 sid 142: Begynnelsevärdesproblem verkar vara vårt ständigt återkommande problem. First 10 Mekanik-Dynamik och Mathematica xAvt Ω2 x t , x ' 0 DSolvex '' t v0 sin t Ω x t v0 , x 0 HH/ITN/BN 0, x t , t Simplify First Ω sid 143: Viktigt! Exempel 6.10 sid 143: Deriveringsövning. Man kan naturligtvis också skriva respektive acceleration, se nästa exempel. t : , och sedan ' t och '' t för hastighet b t, c t2 , d t3 ; D ,t b, 2 c t, 3 d t2 D ,t 0, 2 c, 6 d t eller D , t, 2 0, 2 c, 6 d t Exempel 6.11 sid 143: Integrationsövning med begynnelsevärden som gränser eller med DSolve som måste skedmatas komponentvis DSolve x '' y '' z '' x t t 1 x t b t3 t b t, x 0 x0, x ' 0 v0x, t c, y 0 y0, y ' 0 v0y, t d Sin Ω t , z 0 z0, z ' 0 ,y t ,z t . Simplify 1 6 t v0x 6 x0, y t 6 c t2 t v0x x0, 6 d t v0y v0z Ω v0x c t bt t v0z Ω2 d sin t Ω Ω2 First Ω2 z0 Ω2 z0 '' t 2 2 2 y0, z t ,t d sin t Ω y0, t 2 ' t , bt dtΩ 2 t v0y 2 b t3 c t2 v0z , x t , y t , z t v0y cos t Ω d d Ω Ω c v0z d sin t Ω Exempel 6.12 sid 144: Återigen begynnelsevärdesproblem plus lite efterpyssel. Mathematica klarar lätt system av (ODE). Detta är ett klart gränsfall till kinetik xyAvt DSolve x '' t 0, y '' t g, x' 0 v0 Cos Β , x 0 0, y' 0 v0 Sin Β , y 0 0 , x t ,y t ,t Simplify First g t2 x t t v0 cos Β , y t t v0 sin Β 2 Speciellt kastparabeln. Lös ut y x . ekv xyAvt . Rule Equal . x t x, y t g t2 x t v0 cos Β , y yAvx t v0 sin Β Solve ekv, y, t g x2 y 2 x tan Β 2 v20 cos Β Simplify x 2 ,t v0 cos Β Kastvidden ur kastparabeln. Solve y 0 . yAvx, x Simplify First y HH/ITN/BN Mekanik-Dynamik och Mathematica x v20 sin 2 Β 0 , x 11 g Eller om vi också är intresserade av restiden. Lägg märke till hur Mathematica levererar självdokumenterande. Man behöver inte skriva Svar som i småskolan. xyAvt . Solve y t x0 0, y 0 0 , x 0 . xyAvt, t 2 v0 sin Β v20 sin 2 Β g g Stighöjd och stigtid ur villkoret y ,y 2 v0 sin Β 0 g 0 i högsta punkten. Läget i x-led får vi på köpet! xyAvt . Solve y ' t x Simplify 0 . D xyAvt, t , t v0 sin Β v20 sin Β cos Β g g ,y Simplify v0 sin Β v20 sin2 Β g 2g sid 145-164: Hårdsmält materia! Höghöjdsträning på derivering av vektorer i diverse olika exotiska koordinatsystem. Skumma igenom och fascineras av hur svårt man kan göra enkla saker. Grundidén håller, det vill säga stanna i det vanliga rätvinkliga koordix t , y t , z t alltid är natsystemet och projicera på andra riktningar om så önskas! Bra att veta är att hastighetsvektorn t tangent till bankurvan. Vi sammanfattar de vanligaste ON-systemen, det vill säga koordinatsystem där basvektorerna är enhetsvektorer (Normerade) och parvis vinkelräta (Ortogonala). y 1.0 y Θ 0.8 t r 0.6 x, y, r, t, Θ, n, z Kartesiskt, vårt vanliga xyz–system. z Polärt eller cirkulärt 2D , cylindriskt 3D . Naturligt tangent och normal , t tangentvektor som alltid pekar i parameterriktningen, n pekar alltid in mot krökningscentrum, b t n kallas binormal. b x 0.4 t 0.2 n Vi vet att läget 0.2 0.4 0.6 0.8 t xt,yt,zt 1.0 1.2 1.4 , hastigheten x t xt,yt,zt och accelerationen .. t .. .. .. x t , y t , z t alltid uttrycks i det kartesiska. Om man vill projicera på något av de andra systemet måste först deras ortonormerade basvektorer bestämmas. Vi tar 2D fallet men räknar i 3D som vanligt (vektorprodukt ;-) och följer kokboken. r Pekar i samma riktning som Θ 0, 0, 1 Hastighetsvektorn t så n t . I vårt språk t är alltid tangentvektor till banan . I vårt språk t . ' t t och pekar alltid i parameterrikningen, . t. 0, 0, 1 n t t r r. Θ t t t , så Exempel 6.14 sid 150: Dessa ständiga (BVP). Nu är det plan cirkelrörelse. s t 1 S t 3 c t2 . DSolve S '' t k t3 6 Accelerationsvektorn s' t s '' t , R t2 2 c c 4R och dess belopp a kt 2 k t, 2 Simplify c k t, S ' 0 0, S 0 0 ,S t ,t First 12 Mekanik-Dynamik och Mathematica t4 2 c HH/ITN/BN 4 kt c 16 R2 kt 2 samt slutligen en ögonblicksbild vid t 10s. 3 öb a . c ,R 2, k 400, t 10 10 8801 16 Alltså felräknat i boken. Det symboliska uttrycket (6) stämmer dock med det ovan för a med 4001 8 . Tycker inte om att CN svarar exakt , som dessutom är fel, då k är given på decimal form. Man bör svara på samma form som indata! Som vi vet är Mathemat- ica konsekvent. öb N 5.86335 a . c 2, k 0.3, R 400, t 10 5.86335 sid 158: Ett ofta använt specialfall är cirkulär rörelse i planet, översta tredjedelen av sidan. Lägg speciellt märke till begreppet centripetalacceleration. Det är detta man (förmodligen) menar när man pratar om en utåtriktad centrifugalkraft, vilken inte finns. Uttrycken för får vi lätt med hjälp av vår kokbok ovan. Cos Θ t r , Sin Θ t ,0 cos Θ t , sin Θ t , 0 0, 0, 1 Θ r sin Θ t , cos Θ t , 0 t r r r cos Θ t , r sin Θ t , 0 Först hastighet och acceleration i vårt vanliga kartesiska system. ' t r Θ t sin Θ t , r Θ t cos Θ t , 0 '' t r Θ t sin Θ t r Θ t 2 cos Θ t , r Θ t cos Θ t r Θ t 2 sin Θ t , 0 Sedan med projektion i det polära systemet, som sammanfaller med det naturliga vid cirkelrörelse. ' t . r, ' t . Simplify Θ 0, r Θ t '' t . r, '' t . Θ Simplify r Θ t 2, r Θ t Exempel 6.20 sid 159: På vårt sätt, utan (1)&(2)! Läget av P i det vanliga xyz-systemet ges av t a tΩ a Θ cos t Ω , a Cos Θ , Sin Θ , 0 tΩ .Θ Ωt sin t Ω , 0 Hastigheten är tidsderivatan av läget och är alltid tangentvektor till bankurvan och pekar alltid i parameterriktningen, se figur. ' t a Ω tΩ cos t Ω aΩ tΩ sin t Ω , a Ω tΩ sin t Ω aΩ tΩ cos t Ω , 0 Accelerationen är tidsderivatan av hastigheten, eller andraderiavatan av läget, se figur. HH/ITN/BN Mekanik-Dynamik och Mathematica 13 '' t 2 a Ω2 tΩ sin t Ω , 2 a Ω2 tΩ cos t Ω , 0 Färdig, hastighet och acceleration uttryckt i det vanliga kartesiska xyz-systemet! Om vi speciellt vill ha dessa uppdelade i komposanter längs några andra för problemet naturliga riktningar, CN:s (5)&(6), är det projektion som gäller. Vi vet att r cos Θ , sin Θ , men för sportens skull r t 0, 0, 1 Θ r cos t Ω , sin t Ω , 0 sin t Ω , cos t Ω , 0 Så äntligen . p r, . Θ Simplify . Θ Simplify r .a Ωt a Ω tΩ tΩ ,aΩ r Ω, r Ω . p r, r .a Ωt 2 tΩ 0, 2 a Ω 0, 2 r Ω2 Exempel 6.22 sid 161: Detta är ett typexempel på där derivata (eller implicit dito) fungerar alldeles utmärkt som lösningsmetod! Ställ upp geometrin statiskt vid godtycklig tidpunkt, sedan tar matematiken hand om att generera rörelseekvationen med tecken och allt!! Mycket smidigt! Först Pyttans sats ekv rt 2 2 r t h2 xt x t 2 h2 2 Sedan parar vi ihop med dess tidsderivator och lite trigonometri samt att vi vevar in (därav negativtecknet i r ' t stant fart (r'' t 0). v1 ) med kon- h Solveekv, D ekv, t , D ekv, t, 2 , r' t , Tan Θ v1 , r '' t 0, x t x t , x ' t , x '' t , r t , r ' t , r '' t h v1 ,x t x t h v1 ,x t PowerExpand v1 , r t 0, h ,r t ,r t h cos Θ tan Θ ,r t sin Θ v21 tan3 Θ ,x t Simplify h ,r t h cos Θ tan Θ x t v21 tan3 Θ ,x t sin Θ v1 , r t 0 Här gäller den andra lösningen som CN:s (5)&(8). Den första är en spöklösning under golvet! Lägg märke till CN:s slarv =( med tecken! Åt vilket håll är v p definierad? Mycket bättre att hålla sig till koordinatsytemet så kommer x ' t ut med rätt tecken, här negativ eftersom vagnen rör sig i negativ x-rikting mot origo. Exempel 6.24 sid 163: Även detta är ett typexempel på där derivata (eller implicit dito) fungerar alldeles utmärkt som lösningsmetod! Ställ upp geometrin statiskt sedan tar matematiken hand om rörelseproblemet med tecken och allt!! Mycket smidigt! Med uppenbar rätvinklig triangel får vi hylsans läge y i höjdled. y t : b Tan Θ t Derivera en och två gånger med avseende på tiden för att få hylsans hastighet respektive acceleration, CN:s (5) och (9). y ' t , y '' t bΘ t cos Θ t 2 b Θ t 2 tan Θ t bΘ t 2 , cos Θ t 2 cos Θ t 2 14 Mekanik-Dynamik och Mathematica HH/ITN/BN Vi avslutar med lite reklamfilm Manipulate Graphics Thickness 0.02 , Line 0, 0 , 1.5, 0 , Line 1, 0 , 1, 1.1 Green, Line 0, 0 , 1.5 Cos Θ , Sin Θ , Red, Rectangle 0.95, Tan Θ 0.1 , 1.05, Tan Θ 0.1 , PlotRange 0.1, 1.1 , Θ, 0, 45 , sid 164: Nej tack, verkligen inte! Vi räknar som vanligt! Problem 6.68. Typexempel på problem som löses genom att ställa upp geometrin för en (statisk) ögonblicksbild vid godtycklig tidpunkt och sedan derivera fram dynamiken! Först Θ med rätvinklig triangel tan Θ y x h . vt h Θ t : ArcTan vt Θ ' t , Θ '' t 2 h t v3 hv h2 Simplify t 2 v2 , h2 t 2 v2 2 Va? Enklare kan det inte bli! Så nu en repris för r och annan tillämpning på geometri i rätvinklig triangel, nämligen Pyttans sats r2 x2 y2 r r t : x2 vt y2 2 vt Simplify h2 v2 , h2 h2 . h2 r ' t , r '' t t v2 2 h2 t 2 v2 32 t 2 v2 Exempel: Målarens mardröm. En målare befinner sig på en L m lång stege då dess kontaktpunkt med marken plötsligt släpper och glider ut med konstant fart längs marken. För vår vän på stegen väntar en obehaglig nedfärd. Sök hastigheten för stegens kontaktpunkt mot huset. Lösningsförslag: Lägg in stegen i ett koordinatsystem. Geometrin bestäms av Pytagoras sats. Vi söker hastigheter, vilket är derivata med avseende på tiden, så derivera implicit map t. Om man använder Dt, som deriverar allt som kommer i dess väg, så ser det snyggt ut och precis som när man deriverar för hand med kedjeregeln, men man måste hålla reda på konstanter själv, och framför allt, mycket mera tekniskt när man ska sätta in numeriska värden. Så det rekommenderas D och det är ju inte fel att ha lite koll i sitt modellerande på vad som är beroende på den variabel som man deriverar med avseende på. ekv Dx t 2x t x t 2 2yt y t Lös ut den sökta hastigheten y 2 y t L2 , t 0 y t y längs väggen. Solve ekv, y ' t xt x t y t yt Stämmer ju bra att resan går i negativ riktning då x ökar. T.ex. Stegens längd 5 m ger y då x 3 m och x 2 m/s. Mata in med rätt tecken så kommer svaret ut med rätt tecken i förhållande till de koordinatriktingar vi valt. Så blir det alltid! Smidigt! y . x t 3, x ' t 3 y t 2 2, y t 52 32 HH/ITN/BN Mekanik-Dynamik och Mathematica 15 Exempel: Lille Kalle under lampans sken. Sök hur längden av Kalles skugga ändras då han promenerar mot lampan med konstant fart. Lösningsförslag: Figuren ovan åskådliggör modellen med intressanta geometriska storheter. Det räcker med likformiga trianglar för att koppla det som ändras med tiden, nämligen skuggans längd s och Kalles läge a i förhållande till lampan. Vi söker hastigheter, vilket är derivata med avseende på tiden, så derivera implicit map t och lös ut skuggans ändringshastighet a t ekv s t s t D , t H a t s t s t H s L L Solve ekv, s ' t La t s t H L T.ex. numeriskt: H 10 m, L 2 m och a 2 m/s, negativ eftersom han rör sig mot lampan, ger slutligen s. ALLA indata med tecken ger ALLTID resultatet med rätt tecken! Matematik med vektorer kan man lita på Skuggans längd minskar som sig bör! s . H 10, L 2, a ' t 2 1 s t 2 Exempel: I en ladugård finns en traktordriven höbalslyft enligt figur. Sök höbalens fart upp mot taket då traktorn kör iväg med konstant fart. Lösningsförslag: Med beteckningar enligt figur har vi dels Pytagoras sats och ett samband för repets längd L. Repet utgör ett så kallat kopplingsvillkor mellan x och y. Sådana är mycket vanliga i mekanik. Låt S vara den okända del av linans längd som är upplindad i taljans skivor. ekv h x t 2 xt 2 2 h2 2 lt ,L l t 2 h 2 ,L yt 2 h lt y t l t S S Lös ut y. Notera hur kraftfullt det är att ta med både geometrin och dess derivata i samma Solve. y Solve Join ekv, D ekv, t 1 y t h2 xt 2 , y t , y' t , l t , l' t xt x t 2h L S ,y t 2 1 y t 2 h2 2 h2 xt 2 xt L S ,y t ,l t 2 h2 xt 2 xt 2 xt x t ,l t , 2 xt x t 2h h2 ,l t h2 h2 xt 2 xt 2 xt x t ,l t h2 xt 2 Här duger bara den sista lösningen eftersom y ' t 0. Den andra lösningen är modellens spegelbild under markplanet! Även denna gruvvariant ryms i formuleringen. Återigen får vi ut allt , man ska inte vara snål mot sig själv balens läge och hastighet, i koordinatriktningen, som funktion av traktorns läge x och hastighet x. Var noga med tecken på t.ex. x om du vill exemplifiera med numeriska data. Gör det! 16 Mekanik-Dynamik och Mathematica HH/ITN/BN Exempel: Studera kolvrörelsen enligt figur. Sök läge, hastighet och acceleration för punkterna A, det vill säga kolven, och punkten G på vevstaken som funktion av vinkeln Θ på vevaxeln. Låt det hela veva på med konstant vinkelhastighet Ω 25 rad s. Lösningsförslag: Typisk implicit derivation. Ställ upp geometrin statiskt. Derivation ger sedan smidigt både hastigheter och accelerationer eftersom detta är just tidsderivator av läget med avseende på tiden. Arbeta med ortsvektorer för punkterna. Om vevstaken har längden L så får vi punkten A:s läge med hjälp av Pythagoras sats och två uppenbara rätvinkliga trianglar. L2 A t : B t : r Cos Θ t G t : r Cos Θ t , 0 , Sin Θ t 100 250 t A 250 2 r Sin Θ t B 100 250 t 100 Lite numeriska data. 250 data 100 125 ,r L , Θ' t 1000 25 2 Π, Θ '' t 0; 1000 Hastighet och fart. A' t , B' t , G' 25 Π cos Θ t sin Θ t 32 49 1 400 64 4 32 v 1 400 64 0 25 Π sin Θ t 25 Π cos Θ t sin Θ t 49 Π sin Θ t sin Θ t 4 7 25 . data 2 25 2 t 25 4 4 125 Π sin Θ t 28 125 Π sin Θ t 28 Π cos Θ t Π cos Θ t 2 sin Θ t & Simplify 2 25 Π 5 cos Θ t sin2 Θ t 4 196 25 25 Π 1 , , 4 25 sin2 Θ t 2 25 Π2 cos2 Θ t 2 7Π 196 25 sin2 Θ t sin Θ t 5 Π sin 2 Θ t 25 cos 2 Θ t 367 28 PlotEvaluate v , Θ t , 0, 2 Π , PlotStyle Brown, Blue, Red , AxesLabel "Θ rad ", "vA ,vB ,vG Π 3Π , 2 Π, Automatic All, Ticks , Π, 2 2 PlotRange v A ,vB ,vG m s 20 15 10 5 Π 3Π Π 2 2 Θ rad 2Π Acceleration och dess belopp. A '' t , B '' t , G '' t . data Simplify m s " , HH/ITN/BN Mekanik-Dynamik och Mathematica 625 Π2 759 50 cos 2 Θ t 734 cos Θ t 7340 cos 2 Θ t 4 2 25 cos 2 Θ t 625 2 625 Π2 5313 50 cos 2 Θ t 734 cos Θ t a & 1 4 5 cos 4 Θ t 15 367 0 625 2 25 7 2 25 cos 2 Θ t 50 cos 2 Θ t 367 734 cos 3 Θ t 10 cos 4 Θ t 30 3125 14 32 Π2 sin Θ t Π2 sin Θ t Simplify 625 Π2 759 734 cos 3 Θ t 32 Π2 cos Θ t 14 680 cos 2 Θ t 28 25 50 cos 2 Θ t 17 2 50 cos 2 Θ t 734 cos Θ t 7340 cos 2 Θ t 25 50 cos 2 Θ t 734 cos 3 Θ t 5 cos 4 Θ t 15 2 25 cos 2 Θ t 367 3 625 Π2 , 9 765 625 , 2 1 Π4 sin2 Θ t 196 1568 25 cos 2 Θ t 367 3 390 625 Π4 2 5313 50 cos 2 Θ t 734 cos Θ t 14 680 cos 2 Θ t PlotEvaluate a , Θ t , 0, 2 Π , PlotStyle AxesLabel "Θ rad ", "aA ,aB ,aG 25 7 50 cos 2 Θ t 734 cos 3 Θ t 10 cos 4 Θ t 30 Brown, Blue, Red , m s2 ", Ticks Π 3Π , 2 Π, Automatic , Π, 2 2 a A ,aB ,aG m s2 4000 3000 2000 1000 Π Θ rad 3Π Π 2Π 2 2 Exempel: På ett nöjesfält finns en karusell som roterar med konstant vinkelhastighet Θ, se fig. Banprofilen ges av z h 2 1 cos 2Θ . Bestäm banan på parameterform och rita ut fart och acceleration i varje punkt. Använd R 5 m, h 1 2 1 m och Θ rad s. Lösningsförslag: I tidigare kurs har vi stiftat bekantskap med en funktion av en reell variabel. Detta begrepp kan enkelt spilla över till en vektorvärd funktion av en reell variabel. Vi låter helt enkelt varje komponent i vektorn vara en funktion av samma variabel. Komponenterna kallas då ofta för koordinatfunktioner. En avsikt med konstruktionen är kanske att följa en punkts läge i rummet eller planet som funktion av tiden t. Om vi låter punktens läge beskrivas med ortsvektorn t x t , y t , z t får vi också naturliga definitioner av begreppen definitionsmängd D och värdemängd V . Vi säger att t är en parameterkurva med parametern t. Definition av derivation och integration gör knappast någon förvånad. t t t t xt t , yt t xt , t, zt t , kallas tangentvektor till kurvan yt t, z t t Med figurens beteckningar får vi banprofilen t : R Cos Φ Cos Θ t , R Cos Φ Sin Θ t , z z . Φ ArcSin . z R h 1 Cos 2 Θ t 2 Numeriska data. 1 data R 5, h 1, Θ ' t , Θ '' t 2 Hastighet och fart. 0; 18 Mekanik-Dynamik och Mathematica ' t . data 394 sin Θ t Simplify 7 sin 3 Θ t 3 sin 5 Θ t 394 cos Θ t 9 cos 3 Θ t 3 cos 5 Θ t , 16 v cos2 2 Θ t HH/ITN/BN 2 cos 2 Θ t 99 , sin Θ t cos Θ t cos2 2 Θ t 16 2 cos 2 Θ t 99 Simplify 3144 cos 2 Θ t 2372 cos 4 Θ t 8 cos 6 Θ t Plotv, Θ t , 0, 2 Π , PlotStyle AxesLabel "Θ rad ", "v 79 235 8 cos 8 Θ t 4 cos 2 Θ t cos 4 Θ t 197 Red, Π 3Π , 2 Π, Automatic , Π, 2 2 m s " , Ticks v ms 2.54 2.52 2.50 2.48 2.46 Π Θ rad 3Π Π 2Π 2 2 Beloppet av accelerationen. PlotEvaluate AxesLabel '' t . data , Θ t , 0, 2 Π , PlotStyle rad ", "a "Θ m s2 ", Ticks Blue, Π 3Π , Π, , 2 Π, Automatic 2 2 a m s2 1.34 1.32 1.30 1.28 1.26 1.24 Π 2 Π Θ rad 3Π 2Π 2 En liten rundtur längs en färgglad bana piggar alltid upp. Ju rödare det är ju fortare går det ParametricPlot3DEvaluate t . data , Θ t , 0, 2 Π , PlotStyle v ColorFunction 2.4 Function x, y, z, u , Hue .Θ t 2.54 Thickness 0.03 , 2.4 2 Π u HH/ITN/BN Mekanik-Dynamik och Mathematica 19 Övningar som ska redovisas .. Betyg 3: 6.1-4, 6.7, 6.11(Välj lite numeriska data och rita x, x, x, Θ t, n med projektion), 6.58(först och Betyg 4-5: 6.63(först och i x, y i x, sedan i y r, sedan i Θ r, Km2: Sök yB som funktion av xA . Rita en graf och i x, y sedan i med projektion). med projektion. Animera rörelsen med Manipulate), 6.64, 6.65, Animera rörelsen i 6.11 med Manipulate, samt i Km1 och Km2 nedan. Km1: Sök y som funktion av x. Rita en graf Θ 0, 2Π ), 6.23, 6.26, 6.37, 6.38(först 20 Mekanik-Dynamik och Mathematica HH/ITN/BN Kinetik Att lösa kinetikproblem, det vill säga Newtons rörelseekvationer, ska följa kokboken 1. Placera allt i ett jordfast koordinatsystem. Arbeta alltid i tre dimensioner. Identifiera koordinater för läget av tyngdpunkten för varje del i modellen och frilägg modellen vid godtycklig tidpunkt 2. Friläggning. Det vill säga frilägg modellens alla delar från varandra. Inför kända yttre krafter och moment. Inför okända kontaktkrafter och moment mellan delarna enligt tabell sid 61–62, dessa blir alltid på den ena delen och på den andra. Rita tydliga figurer över delarna med alla kraft– och momentpilar tydligt markerade, gärna i olika färger. 3. Uppställning av rörelseekvationer för läget av varje del vid godtycklig tidpunkt med avseende på xyz–riktningarna i m .. Tyngdpunktens rörelse. .. Θ Rotation kring tyngdpunkten. i i i Observera att momentekvationerna skall formuleras som rotation kring tyngdpunkten. Använd inte godtycklig momentpunkt, det ger fel resultat Valfri momentpunkt gäller bara vid statik då Vi får alltså 2 vektorekvationer för varje del, Fx , F y , Fz .. m och Mx , M y , Mz i 0. .. Θ. Om vi har n st kroppar får vi 6n skalära ekvationer att lösa. 4. Identifiering av begynnelsevillkor BV och eller randvillkor RV . 5. Lösning av rörelseekvationerna som ett BVP är inte mekanik utan matematik. Så använd enbart kunskaper från kurs i ODE och DSolve i Mathematica. I praktiken så är enda ändringen från statik till dynamik att man byter nollan i .. .. högerledet i jämviktsekvationerna till m respektive Θ och Solve till DSolve 6. Tolkning av lösning arna , dimensionsanalys, grafer, ställ frågor till modellen Egentligen var det inte m .. som Newton indikerade i sin andra lag, utan förändringshastigheten av den sk rörelsemängden m är lika med den pålagda kraften, dvs t m m m .. . Oftast betraktar man massan konstant under resan, exempelvis för en boll som kastas iväg eller en bil som accelererar, vilket leder till den välkända accelerationslagen. För en rymdraket däremot gäller inte detta eftersom den eldar upp det mesta av sig själv under de första minuterna efter start, så m kan inte försummas. Mathematica är en ovärdelig hjälpreda under punkt 2 och 4. Man kan helt koncentrera sig på punkterna 1, 3 och 5, det vill säga modellering och mekanik. Så precis som vid statik är det alltså friläggning som man absolut måste vara duktig på själv! Om man gör en så kallad Lagrangeformulering kan man enkelt få Mathematica att generera rörelseekvationerna. Kanske hinner vi med ett exempel längre fram. Innan du går vidare så repetera inledningen av kinematiken i detta häfte. sid 165-175: Läs igenom. Rätt bra för att vara CN. Inertialsystem kan ni läsa som vårt enda vanliga jordfasta xyz-system. sid 176: Näe tack, verkligen!! sid 177-180: Läs, inget nytt! Precis enligt kokboken ovan, dvs frilägg och sätt upp rörelseekvationerna för varje kropp. Som illustration på detta börjar vi med en kropp. Exempel: Inbromsning med hjälp av friktion, enligt Newton. Sök bromssträcka och tid. Ta för vana att först, om möjligt, alltid göra en symbolisk lösning. Detta lär dig mycket mekanik, dimensionskontroll med mera samt en djupare kunskap om just den problemtyp du löser Till slut tränar man storleksordning och rimlighet genom att sätta in lite numeriska data. .. Lösningsförslag: Med Newton mx xAvt DSolve 1 x t 2 2 t v0 F får vi direkt följande begynnelsevärdesproblem (BVP) och dess lösning. m x '' t Μk m g, x 0 0, x ' 0 g t2 Μk Först inbromsningstid till stillastående. tid Solve x ' t 0 . D xAvt, t , t First v0 , x t , t First HH/ITN/BN Mekanik-Dynamik och Mathematica v0 t 21 g Μk Sedan bromssträcka. xAvt . tid v0 v20 g Μk 2 g Μk x Lägg märke till det självdokumenterande svaret! Avslutningsvis lite numerik. xAvt . tid . v0 x 1.78389 7, g 9.81, Μk 0.4 6.24363 Exempel 7.1 sid 180: Ekv (1) är typexempel på så kallad solidifiering. Man ska inte behöva lägga ihop 3 2 1 6 i huvudet! Istället gör vi det vi sagt, nämligen frilägga och ställa upp rörelseekvationerna för varje kropps tyngdpunkt, det vill säga friläggning och Newton för varje kropp. Som vanligt placerar vi ett koordiantsystem med origo fast i en lyktstolpe till vänster om "P-pilen" och x-axeln pekandes åt höger. Med krafter enligt bild får vi med start från vila (BVP) och dess lösning. klossar DSolve P t2 x3 t N1 t2 2 m3 m3 x3 '' t P N1 , x3 ' 0 0, x3 0 0, m2 x2 '' t N1 N2 , x2 ' 0 0, x2 0 0, m1 x1 '' t N2 , x1 ' 0 0, x1 0 0 , x1 t , x2 t , x3 t , t First N1 t2 , x2 t N2 t2 2 m2 N2 t2 , x1 t 2 m1 Med hjälp av kopplingsvillkoren att klossarna rör sig tillsammans kan vi nu så lösa ut kontaktkrafterna, som mycket riktigt blir oberoende av tiden. ekv x3 t P t2 x2 t , x2 t N1 t2 N1 t2 N2 t2 N1 t2 N2 t2 , 2 m2 2 m2 2 m3 kontakt N1 x1 t . klossar N2 t2 2 m1 Solve ekv, N1 , N2 m1 m1 m2 P m2 m3 m1 P , N2 m1 m2 m3 Notera nyttan av att räkna symboliskt!! Slutligen CN:s uppsättning klossar. kontakt . m1 P N1 2 1, m2 2, m3 3 P , N2 6 sid 181-186: Naturligtvis räknar vi alltid i det vanliga, dvs det vänstra i (7.18) och (7.19). I övrigt handlar det om att sätta upp BVP ODE BV . Kopplade system är inget problem i Mathematica, det använde vi ju i ex 7.1 ovan. Nyttiga exempel. Vi tar två av dem, de andra gör du själv! Exempel 7.3 sid 182: Låt partikelns läge vara boll x t DSolve m x '' t m y '' t 6 m t v0 cos Α t xt,yt så med (BV) från text får direkt (BVP) och dess lösning. k y t , x' 0 v0 Cos Α , x 0 0, 0, y ' 0 v0 Sin Α , y 0 0 , x t ,y t k t3 v0 sin Α ,yt 6m ,t First t v0 sin Α Bankurvan Solve boll . Rule y x tan Α Equal . x t k y3 6 m v20 sin Α x, y t y , x, t Simplify y 2 ,t v0 sin Α Slutligen hittar vi på lite data och ritar. Till vänster läget, i mitten farten och till höger accelerationen under de T första sekunderna. 22 Mekanik-Dynamik och Mathematica HH/ITN/BN data m 1, v0 1, k 1, Α 20 ; T 5; t x t ,y t . boll . data; ParametricPlot Evaluate t , t, 0, T , PlotStyle Brown, AxesLabel "x", "y" Plot Evaluate ' t , t, 0, T , PlotStyle Blue, AxesLabel "t", "v" , Plot Evaluate '' t , t, 0, T , PlotStyle Red, AxesLabel "t", "a" v y a 3.5 3.0 2.5 , 2.0 1.5 1.0 0.5 1.5 1.0 0.5 1.5 , 1.0 1 1 0.5 x 2 , 0 1 2 3 4 t t 5 1 2 3 4 5 Exempel 7.4 sid 184: Låt den hängande massan ha koordinaten s t räknat från rullen. I övrigt det sneda xy-systemet. Så med (BV) från text får vi direkt (BVP) och dess lösning. Givetvis använder vi stora F isf f för friktionskraft CN! lådor DSolve t2 F 4 m x '' t 4 m y '' t 3 m s '' t 4 g m sin Β S 4 m g Sin Β F, x ' 0 0, x 0 0, N 4 m g Cos Β , y ' 0 0, y 0 0, 3 m g S, s ' 0 0, s 0 0 , x t ,y t ,s t t2 S x t 4 g m cos Β N ,yt 3g m t2 S First ,s t 8m ,t t2 6m 8m Med villkoret att snöret håller, övre lådan stannar på banan och fullt utbildad friktion (annars rör det ju inte på sig ;) kan vi lösa ut FNåS Solve x t s t ,y t 4 g m cos Β Μk , N 4 g m cos Β , S 0, F Μk N . lådor, F, N, S First 12 F g m sin Β g m cos Β Μk 7 gm Vi ser att alla dessa är oberoende av tiden. Så är det inte alltid, speciellt inte för snörkrafter. Nu kan vi meka ihop en "riktig" lådlösning som endast är funktion av tiden. lådor lådor . FNåS 1 t2 x t 8m yt 12 4 g m cos Β Μk 1 12 6m 7 0, s t t2 g m cos Β Μk 7 g m cos Β Μk g m sin Β gm gm 3 g m t2 g m sin Β 4 g m sin Β , CN:s lösning går via fasplanet för att snabbast nå målet gällande hastighet som funktion av läget. Detta kunde vi också gjort, men det funkar alltid att stanna i tidsplanet så här. På köpet får vi då även restider. Speciellt restiden T till läget l där CN undrar över farten. T Solve x t l . lådor, t 14 l 14 l , t t 4 g sin Β 4 g cos Β Μk 3g 4 g sin Β 4 g cos Β Μk 3g Eftersom vi inte befattar oss med negativa restider duger bara den andra lösningen. Så svaret på CN:s undran om lådornas fart då de glidit sträckan l. (Lägg märke till det självdokumenterande svaret.) D lådor, t .T 2 14 Simplify l 2 l x g 4 sin Β 4 cos Β Μk 14 3 g 4 sin Β 7 l 14 y 4 cos Β Μk l 3 , 2 0, s g 4 sin Β 4 cos Β Μk 3 l g 4 sin Β 4 cos Β Μk 3 7 g 4 sin Β 4 cos Β Μk 3 sid 187-198: Hoppa över cirkelrörelse med polära koordinater! Du har nog förstått vid det här laget att vårt vanliga xy-system fungerar här med. Lösenordet är projektion. Det räcker med det tunga kinematikpasset! HH/ITN/BN Mekanik-Dynamik och Mathematica 23 Vi kortsluter nu det mesta av resten i boken. Jag sammanfattar med lite nedslag här och var i den. Energilagar - Arbete, Potentiell och Kinetisk energi, Effekt sid 199-200: Läs igenom! Energilagar är ibland smidigt att använda vid problemlösning, men inget nytt! Vi kan lika väl köra på med våra (BVP) enligt ovan eftersom energilagarna är härledda från Newton! Om inte det fungerar så gör inte energilagar det heller! sid 202-203 "Om partikelns , 204: Arbete är kraft gånger väg som naturligt beräknas med skalärprodukt. Exempel: Sök det arbete som kraften om 100 N uträttar under skådespelet som återges i figuren. Lösningsförslag: Arbete är ju lika med kraftens storlek i vägens riktning gånger sträckan, det vill säga A F cos Α s vilket inte är något annat än skalärprodukt om vi betraktar både kraft och väg som vektorer A . Mathematica arbetar alltid med vinklar i Π radianer. Omräkning av grader till radianer görs som vanligt med hjälp av 180 eller vackrare med deg som då resulterar i ett . Först kraften och vägen på vektorform 100 Cos 60 50, 50 , Sin 60 3 50, 0 0, 0 50, 0 sedan arbetet med skalärprodukt . 2500 y Exempel: Genom att använda integral kan vi även ta hand om fallet då både väg x och kraft varierar under resan. Sök det arbete som kraften , 1 N uträttar då den släpar en grön boll uppför cosinusbacken y x 1 cos x m, x 0, Π . 3.0 x x 2.5 ,1 2.0 1.5 1.0 0.5 0.0 0.0 0.5 1.0 1.5 2.0 2.5 3.0 x Lösningsförslag: Här varierar både kraft och väg med parametern x. Genom att ta på oss integralglasögen inser vi att en studie av en liten del av resan med efterföljande sammanslagning leder till målet. Så under den lilla förflyttningsvektorn vid x uträttas det lilla arbetet A x, y Π 0 får vi A y 1, x A 0 A x ,1 x x. Lägg sedan samman alla små bidrag som vanligt A 1, sin x Π 0 x x sin x x x cos x Π 0 A Π 3 A 3 y x x. Här . Eller med Mathematica x , 1 . 1, D 1 Cos x , x x TrigToExp Π A x, y Fx s , F y s Exempel: Vilket arbete krävs för att dra ut en fjäder om man vet att kraften 400 N drar ut den 3 100 y x , s 1 20 s s A A A x s , y s . Vi får då sb Fx sa s , Fy s x , s y s s. m m? Lösningsförslag: Låt fjädern vara utdragen x m. Det lilla arbetet att dra ut den ytterligare ett litet stycke F x 1, 0 Naturligtvis går det bra att generalisera till allmänna kurvor på parameterform s A b a Π A x Π 0 y 1, x kx x 400 x 3 100 x och slutligen A A x blir då 24 Mekanik-Dynamik och Mathematica A 400 1 20 A x x 0 HH/ITN/BN 3 0 100 50 A 3 Vi har två typer av energi som mäts i Nm i mekanik. När det gäller värme används J eller Ws, eller ännu vanligare kWh. 1 2 Kinetisk energi T: T 1 2 IΘ , 2 mv2 . där som vanligt m är kroppens massa och v dess tyngdpunktsfart, I tröghetsmomentet och Θ dess rotationshastighet kring tyngdpunkten. Potentell energi V: Om de krafter som verkar på kroppen är oberoende av bankurvan kallas de konservativa krafter. Två sådana viktiga tillämpningar är Tyngdkraft (överst sid. 212, mitten sid. 213) V mgh. Man brukar tala om lägesenergi. 1 kx2 , 2 Fjäderkraft (överst sid. 214) V där k är fjäderkonstanten och x den aktuella förlängningen. Se föregående exempel. Om det inte försvinner någon energi under rörelsen, t.ex. friktion, gäller att den totala mekaniska energin är konstant T V T0 V0 Exempel: En sten släpps från ett fönster på höjden h ovan mark. Sök dess fart v då den når mark! Vi får direkt med marken som referensnivå Tf Vf Tm Vm 0 1 2 mgh mv2 0 v 2gh . Vi kontrollerar med vårt vanliga artilleri. Låt y-axeln vara riktad nedåt med origo vid fönstret. Vi får då restiden på köpet. yAvt DSolve m y '' t m g, y ' 0 0, y 0 0 ,y t ,t First g t2 y t 2 Restiden till marken T Solve y t 2 h . yAvt, t h t 2 h , t g g Samma visa med tiden hela tiden. En av fysikens förbannelser. All existerande teori, Newton, Maxwell eller kvantmekanik, kan bara leverera t2 . Fysiker kan alltså inte avgöra om tiden går framåt eller bakåt men vi vet D yAvt, t 2 .T 2 h 2 y g h g Effekt: P Fv, där som vanligt F är kraften som verkar på kroppen och v dess tyngdpunktsfart. Effekt mäts i W. Äldre enhet är hästkraft (hk eller hp). Momentekvationen - rotation kring en fix axel sid 297: Bra, läs igenom! sid 305 tom (14.18): Vi ska återvända till kokboken och studera den till m .. .. , det vill säga tyngdpunktens rörelse längs koordi- nataxlarna, analoga Θ som är kroppens vridning kring sin tyngdpunkt runt koordinataxlarna. Här är och tröghetsmomentet (egentligen masströghetsmomentet). sid 311-316: Beräkning av tröghetsmoment 2 m det pålagda momentet m. Ungefär samma resa som beräkning av tyngdpunkt. Finns även i samma tabell. Tröghetsmomentet brukar ibland anges genom den så kallade tröghetsradien k, som I mk 2 . För att inte förväxla tröghetsmomentet med den komplexa enheten (på gammal form I) i Mathematica brukar man använda beteckningen J i stället, som för övrigt får anses vara synonymt eftersom det används i många mekanikböcker. HH/ITN/BN Mekanik-Dynamik och Mathematica Exempel: Bestäm tröghetsmomentet m r2 m för en 25 y homogen cylinder med radien R och massan m, med avseende på rotationsaxeln. x R Lösningsförslag: Låt cylindern ha längden L. Först har vi volymdensiteten Ρ tunna lökringar som vid radien r har arean 2Πr tröghetsmomentet från ett sådant är J r2 m J 0 2 r 0 . Dela sedan upp den cirkulära tvärsnittsytan i r. En sådan tunn lökring växer till ett tunnväggigt rör med längden L. Bidraget till r2 ΡL 2Πr r. Nu är det bara att lägga samman. m R J m ΠR2 L L2Πr r Π R2 L m R2 J 2 Exempel: Bestäm tröghetsmomentet m r2 m för en tunn rektangel med bredden a, höjden b och massan m, med avseende på en axel längs kanten b. y b x a Lösningsförslag: Först har vi ytdensiteten Ρ tröghetsmomentet från en sådan är J J a J 0 2 x 0 2 x m m . Klipp sedan upp rektangeln i smala rektangulära ab 2 x Ρb x. Nu är det bara att lägga samman. strimlor b x. Bidraget till m b x ab a2 m J 3 y Exempel: En tunn pappskiva i form av en rätvinklig triangel med massan m är uppriggad enligt figur. Sök tröghetsmomentet m r2 m då den roterar kring y–axel. b x a Lösningsförslag: Först har vi ytdensiteten Ρ m 1 2 rektangulära strimlor y x samman. J 0 2 x 0 x. Bidraget till tröghetsmomentet från en sådan är J m a J 2 ab b1 x2 m x . a Klipp sedan upp triangeln i smala x2 Ρy x x. Nu är det bara att lägga x b 1 1 och hypotenusans ekvation y x ab x a a2 m J 6 Exempel: En tunn tråd med densiteten Ρ böjs till en spiral kΘ , 0 Θ 4Π. Bestäm spiralens med radien r Θ tröghetsmoment kring origo. y 1.5 1.0 k 0.05 0.5 1.5 1.0 0.5 0.5 0.5 1.0 1.5 x 1.0 1.5 Lösningsförslag: Klipp upp spiralen i små bitar s r Θ. Det lilla tröghetsmomentet ges sedan av J Sedan är det bara att lägga samman alla de små bidragen r2 m r2 Ρ s r2 Ρr Θ. 26 Mekanik-Dynamik och Mathematica J 4Π J Ρ 0 kΘ 3 HH/ITN/BN Θ 0 12 Π k 1 Ρ J 3k Exempel: En homogen cylinder med massan M och radien R är lagrad i sin rotationsaxel. Med hjälp av ett upplindat snöre accelereras den igång från stillastående med hjälp av en i snöränden hängande sten med massan m. Vi har alltså situationen återgiven i figuren till höger. Studera rörelsen som funktion av tiden. Θ S y mg Lösningsförslag: I problemtextens figur har vi passande nog frilagt de två kropparna och infört koordinater för deras lägen, snörkraften S och den verkande tyngdkraften mg. Med nyligen beräknat tröghetsmoment för en cylinder får vi så (BVP) baserat på rörelseekvationerna. 1 ΘÅy DSolve M R2 Θ '' t R S, Θ 0 0, Θ ' 0 0, 2 m y '' t S t2 g m t2 Θ t mg S, y 0 0, y ' 0 0, Θ t , y t , t First S t2 ,yt MR 2m Kopplingsvillkoret att snöret håller bestämmer S och knyter samman Θ och y. snörkraft Solvey t RΘ t . ΘÅy, S First gmM S 2m M Så till slut de två kropparnas rörelse som funktion av tiden ΘÅy ΘÅy . snörkraft g m t2 Θ t Simplify g m t2 ,yt 2mR MR 2m M Svängningar sid 260-281: Vår omvärld är full av företeelser som karakteriseras av svängningar, exempelvis vibrationer i en maskin, stötar i ett knä vid löpning eller brus i en mobiltelefon. Många av dessa kan med fördel lineariseras till en diskret modell med en massa, en fjäder och en dämpare. Inte sällan ser man flera sådana "atomer" sammankopplade i större strukturer. Med hjälp av Newtons rörelselagar får vi så en modell som bör finnas i verktygslådan hos varje praktiserande ingenjör. För fjädern gäller att kraften som krävs för att förlänga den en sträcka x är proportionell mot denna, det vill säga F f kx, där k kallas fjäderkonstanten med enheten N/m. Observera att uttrycket gäller med tecken på x, det vill säga även då den trycks ihop. Ff k x Ff x För dämparen gäller däremot att kraften är proportionell mot förlängningshastigheten (farten), det vill säga Fd cx , där c kallas dämpningskonstanten med enheten Ns/m. En dämpare består i princip av en cylinder fylld med olja. En rörlig bricka med hål i delar in cylindern i två kammare. Oljans viskositet, det vill säga hur trögflytande den är, antal hål i brickan och deras storlek avgör hur trögt det är att flytta den, eftersom olja från ena kammaren då skall flyttas till den andra. Detta ger dämparen dess karakteristiska funktion som också brukar kallas viskös dämpning. Kraften som krävs för att flytta på brickan är alltså inte beroende av dess läge HH/ITN/BN Mekanik-Dynamik och Mathematica 27 utan endast på hastigheten (farten). Jämför potatispress! Fd c Fd x x Massa och fjäder utan yttre last (fri odämpad svängning) Vi börjar med en massa kopplad till en fjäder. Vi frilägger och noterar att den motriktade kraften i fjädern är proportionell mot läget med fjäderkonstanten k. Formulera nu rörelseekvationen med hjälp av Newton .. .. kx mx mx kx .. 0 k x m x 0 k m Vi känner igen en linjär andra ordningens (ODE) med konstanta koefficienter. Här kallas Ωe enheten [rad/s]. Man pratar också om egenfrekvensen fe Ωe 2Π och egensvängningstiden Te 1 . fe för egenvinkelfrekvens med Prefixet "egen" kommer av att Ωe endast beror på systemets "egna" inneboende parametrar. Lösningen är en så kallad harmonisk svängning (sinusformad) som fortgår i all evighet om systemet störs från sitt jämviktsläge. Notera speciellt var Ωe dyker upp! DSolve m x '' t kx t , x t , t k t x t k t c2 sin c1 cos m m Som väntat får vi två konstanter som fixeras av (BV). Uttrycket ovan brukar ofta skrivas om som en ren sinus- eller cosinusvåg. a sin Ωt a2 b cos Ωt a b2 a2 och identifierar amplitud R och fasvinkeln a2 a2 b2 cos a2 b2 sin Ωt b sin Ωt b2 a2 sin Ωt b2 och fasvinkel sin cos Ωt 1 b2 a a2 cos Ωt b 1, 1 b2 a2 1 b2 Summaformel för sinus a som naturligtvis mäts i radianer. Punkten R , bestäms likt argumentet för ett komplext tal. Var och en av de oändligt många vinklar b R ligger på enhetscirkeln som löser ekvationerna a R cos b R sin kan göra anspråk på att kallas för fasvinkel. På grund av periodiciteten hos cosinus och sinus skiljer de sig åt med en multipel av 2Π så alla ger "samma effekt". Ofta nöjer man sig med den så kallade principalvinkeln som ligger i intervallet Π, Π . När man räknar för hand gäller det att se till så man hamnar i rätt kvadrant!! Om a 0 eller b 0 är det ju enkelt annars går det bra att beräkna b arctan a fasvinkeln som Π eftersom arctan() levererar vinklar i första och fjärde kvadranten. Den avslutande korrektionen om a 0 kommer sig naturligtvis av att vi kan ha dividerat bort "negativ" information, ty b a b a och b a b . a I Mathematica är det inga problem, speciellt inte om vi använder versionen med två argument ArcTan[a,b]som alltid levererar rätt vinkel i intervallet Π, Π och givetvis i radianer. Vi tar ett exempelvis med m 2 kg, k 18 N/m och begynnelsevärdena x 0 1 och x 0 1. xAvt DSolve 2 x '' t 18 x t , x 0 1, x ' 0 1 ,x t ,t First 1 x t sin 3 t 3 cos 3 t 3 Plot x t . xAvt, t, 0, 10 , PlotStyle Red, AxesLabel "t", "x t " 28 Mekanik-Dynamik och Mathematica HH/ITN/BN xt 1.0 0.5 t 2 4 6 8 10 0.5 1.0 a, b Coefficient x t . xAvt, Sin 3 t , Cos 3 t 1 , 1 3 a2 R b2 10 3 ArcTan a, b tan 1 3 En sista ängslig koll R Sin 3 t x t . xAvt Simplify 0 Massa, fjäder och dämpare utan yttre last (fri odämpad svängning) Massa kopplad till en fjäder och dämpare. Vi frilägger och noterar att den motriktade kraften i fjädern är proportionell mot läget och i dämparen mot hastigheten med dämpningskonstanten c. Formulera rörelseekvationen med hjälp av Newton .. mx kx cx .. mx kx cx .. 0 c x m x k x m 0 Vi känner igen en linjär andra ordningens (ODE) med konstanta koefficienter. Beroende på om rötterna till den karakteristiska ekvationen är reella och olika, reella och lika respektive komplexkonjugerade använder vi beteckningen överkritisk, kritiskt respektive underkritisk dämpning. Detta avgörs av den nya aktören, dämparen, som i praktiken bidrar med att "äta energi", vilket ger sig tillkänna som en exponentiellt avtagande amplitud hos lösningen. Även den tidigare nämnda egenvinkelfrekvensen Ωe c kommer med i lösningen på "samma ställe" som tidigare. Man brukar också definiera dämpningsfaktorn Ζ genom 2Ωe Ζ m . För de tre fallen ovan gäller Ζ figurer. 1, Ζ 0 respektive Ζ 1. De tre möjliga fallen av homogena lösningar illustreras kvalitativt i följande xt xt xt t t t Överkritisk dämpning. Inte så vanligt förekommande i ingenjörssammanhang. Exempelvis underkritisk dämpning med m Kritisk dämpning. Lösningskurvan skär t axeln högst en gång. Typiskt utseende för en ny hjulupphängning på en bil. 3 kg, c 1 Ns m, k Underkritisk dämpning. Mycket vanligt insvängningsförlopp. Typiskt utseende för en gammal hjulupphängning på en bil. 5 N m med (BV): x 0 1 och x 0 1. HH/ITN/BN Mekanik-Dynamik och Mathematica xAvt DSolve 1 x t 3 x '' t 5x t 1 x' t , x 0 59 t t6 7 1 ,x t ,t 59 cos 6 6 59 t R ab.ab . ab Coefficientx t 59 t , Cos . xAvt, Sin 6 3 First 59 t 59 sin 59 1, x ' 0 29 6 t3 6 59 Plot R, x t . xAvt, R , t, 0, 20 , PlotStyle PlotRange All, AxesLabel "t", "x t , R t " Blue, Red , xt, Rt 1.5 1.0 0.5 5 10 15 t 20 0.5 1.0 1.5 Massa, fjäder och dämpare med yttre last (påtvingad svängning) Avslutningsvis kan man naturligtvis utsätta de båda modellerna ovan för en yttre, ofta kallad störande eller påtvingad, kraft. I nästan alla intressanta applikationer är denna av periodisk natur F t F0 sin Ωt . Vi väljer att exemplifiera med en full massa–fjäder–dämpare–modell. Formulera nu rörelseekvationen med hjälp av Newton .. mx kx cx .. Ft mx kx cx Ft .. x c x m k x m Ft m Från kurs i (ODE) känner vi till att lösningen delas upp i "systemets inneboende" homogena lösning xh samt en partikulärlösning x p beroende på högerledet, ofta kallad fortvarighetslösning (eng. steady-state solution). Allmänna lösningen blir sedan x xh x p. Homogena lösningen brukar kallas för transienta lösningen (eng. transient solution) eftersom den klingar av ganska snabbt med tiden. Kvar blir x p , därav namnet. När vi löser en differentialekvation med Mathematica får vi inte denna, ibland önskade, uppdelningen. Ett trick är att lösa systemet utan begynnelsevärden. De termer som då har "konstanterna" Ci på sig tillhör den homogena lösningen och öppnar därmed för en separation av lösningen. Exempelvis med m (BVP). xBVP 1 kg, c DSolve 1 x t t 20 2 Ns m, k 5 N m, F0 x '' t 2 x' t x 0 0, x ' 0 t 37 sin 2 t 80 10 N, Ω 2 rad s med (BV): x 0 5x t 10 Sin 2 t , 1 ,x t ,t First t 0 och x 0 Simplify 1 cos 2 t 34 Plot x t . xBVP, t, 0, 15 , PlotStyle Red, AxesLabel xt 2 1 2 1 2 4 6 8 10 12 14 t "t", "x t " 1. Först hela lösningen till 30 Mekanik-Dynamik och Mathematica Nu över till separationen x xODE 1 t 17 0, c2 xp x p . Samma som ovan men utan begynnelsevärden (BV). DSolve x '' t x t Så c1 xh 17 c1 HH/ITN/BN 2 x' t 10 t sin 2 t 5x t 10 Sin 2 t , x t , t First Simplify 40 t cos 2 t 17 c2 0 ger oss partikulärlösningen x p och sedan xh x t . xODE . c1 0, c2 0 Simplify 10 4 cos 2 t sin 2 t 17 xh x t 1 t . xBVP 37 sin 2 t xp Simplify 80 cos 2 t 34 Nu kan vi inspektera de olika delarna. Allmänna lösningen anpassar sig raskt till partikulärlösningen eftersom den homogena lösningen snabbt klingar av med tiden beroende på exponentialfunktioner med negativa exponenter. Plot x t . xBVP, xh, xp , t, 0, 15 , PlotStyle AxesLabel "t", "x t ,xh t ,xp t " Red, Blue, Green , x t ,xh t ,x p t 2 1 2 4 6 8 10 12 14 t 1 2 Även här brukar fortvarighetslösningen skrivas om till ren sinusvåg x p t vinkelfrekvens Ω, här lika med 2, och fasvinkel a, b 10 i radianer naturligtvis. Coefficient xp, Sin 2 t , Cos 2 t 40 , 17 17 a2 R b2 10 17 ArcTan a, b tan 1 4 En sista ängslig koll R Sin 2 t xp Simplify 0 Amplitud över alla gränser Exempel: En enkel massa med fjäder och dämpare är uppriggat enligt figur. Sök fortvarighetslösningens amplitud orsakad av den störande kraften. asin Ωt bcos Ωt Rsin Ωt med amplitud R, HH/ITN/BN Mekanik-Dynamik och Mathematica 31 k m Lösningsförslag: Om den störande kraftens vinkelfrekvens Ω närmar sig egenvinkelfrekvensen Ωe kommer lösningens amplitud att öka dramatiskt. Man talar om resonans (eng. resonance). Vid konstruktion med liten dämpning ska man alltså undvika detta tillstånd eftersom ett haveri är att vänta. Vi illustrerar med en bukett dämpare samt Ωe 1 och ritar upp partikulärlösningens amplitud R. amp ; Doxt x t r . DSolve x '' t c x' t Coefficient Simplify First xt AppendToamp, 3 Sin Ω t , x t , t , Cos Ω t , Sin Ω t . c1 0, c2 0 ; Chop; r.r 5 , 1, , c, x t 10 ; 100 Plot Evaluate amp , Ω, 0, 2 , PlotStyle "Rainbow", PlotRange 0, 4 , AxesLabel "Ω Ωe ", "R" R 4 3 2 1 0.0 0.5 1.0 1.5 2.0 Ω Ωe Vi poängterar ytterligare de två fundamentala och viktiga fallen då den störande vinkelfrekvensen Ω blir extra tydlig om vi illustrerar utan dämpare. Fallet Ω Ωe och Ω Ωe . Effekten Ωe för modell utan dämpare Speciella tillämpningar har vi då Ωe xAvt DSolve Ω är mycket mindre än Ωe x '' t 36 x t x 0 0, x ' 0 Sin 5 t , 1 ,x t ,t Ω, så exempelvis med Ωe First 6 och Ω 5. Simplify 1 x t sin 5 t sin 6 t 11 Denna kan, med ett litet trick, skrivas om till xt TrigFactor x t 2 11 t . 6 t cos sin 11 . xAvt . 6 2 2 Här kan vi nu betrakta 2 t cos 2 11 som en med tiden "sakta" varierande amplitud till den "fortare" svängande sinusvågen sin 11t . 2 Bland musiker brukar detta kallas för svävning (eng. beat) och kan höras t.ex. strax innan en gitarrsträng är riktigt stämd. Elektroingenjörer identifierar fenomenet som amplitudmodulering (eng. amplitude modulation) och hittas på radions AM-band. Där känner vi också igen FM-bandet baserat på frekvensmodulering (eng. frequency modulation). 2 Plotxt, t Cos , 11 2 AxesLabel 2 t Cos , t, 0, 6 Π , PlotStyle 11 2 2 t cos " "t", "x t , 11 2 Red, Blue, Blue , 32 Mekanik-Dynamik och Mathematica 2 HH/ITN/BN t xt, cos 2 11 0.2 0.1 t 5 10 15 0.1 0.2 Fallet Ω Ωe för modell utan dämpare Detta resonansfall måste ovillkorligen undvikas vid all konstruktion eftersom det leder till haveri. Samma som föregående med Ωe Ω 6. xAvt DSolve 13 x t x '' t 36 x t x 0 0, x ' 0 Sin 6 t , 1 ,x t ,t First Simplify Apart 1 sin 6 t t cos 6 t 72 12 Här kommer nu partikulärlösningens tidsberoende amplitud allt mera kraftiga svängningar eftersom t Plotx t . xAvt, t 12 då t t 12 att snabbt dominera förloppet. Om inget hejdar utvecklingen får vi , med ett oundvikligt haveri som följd. t , 12 , t, 0, 6 Π , PlotStyle Red, Blue, Blue , 12 t AxesLabel " "t", "x t , 12 t xt, 12 1.5 1.0 0.5 5 10 15 t 0.5 1.0 1.5 Ett gammalt välkänt exempel på denna typ av haveri är kollapsen av Tacoma Narrows Bridge 1940, där den måttliga vindens frekvens matchade längden på bron och satte en torsionsmod i resonanssvängning. På www hittar man lätt filmsnuttar från skådespelet! Ett exempel från senare tid är rymdfärjan Challenger där en turbopump oavsiktligt konstruerades med resonans. Lyckligtvis upptäcktes det före start, men felet kostade NASA flera miljoner dollar. Vid olyckan 1985 kom problemställningen att diskuteras på nytt. HH/ITN/BN Mekanik-Dynamik och Mathematica 33 Övningar som ska redovisas Betyg 3: 7.11, 7.14, 7.15, 7.20, 14.23, 14.25, 12.1. Betyg 4-5: 7.19, 14.26, 12.44, samt Kt1 och Kt2 nedan. Kt1 : En solid pappersrulle lastas osurrad på ett lastbilsflak enligt figur. Så kör lastbilen iväg med accelerationen a. Hur lång tid tar det och hur långt har då lastbilen hunnit när papersrullen rullar av flaket? d Kt2 : Studera en enkel robot som bara kan röra sig xy–planet. Låt armarnas längder vara L1 1 m, L2 0.8 m och Θ Θ1 , Θ2 , där Θ1 , Θ2 0, Π rad. Roboten ska måla längs en bana som beskrivs av den räta linjen från 0.9, 0.1 m till 1.2, 1.1 m. Det är viktigt att farten 0.1 m s hålles, annars blir färglagret för tunt eller för tjockt. a Bestäm framåtekvationerna, det vill säga P Θ xΘ,yΘ . b Bestäm restiden T och dela in denna i n st små tidssteg. c Stega fram för i 0, , n och beräkna först position b och y L2 P Θ2 L1 Θ1 hastighet b från banan sedan robotens fyra styrvariabler Θ, Θ ur x det olinjära ekvationssystemet b P, b P med FindRoot. d Rita ut robotens fyra styrvariabler Θ, Θ som funktioner av i 0, , n med ListPlot. Animera med Manipulate ť Blandade tillämpningar för resten Exempel: En sprinter startar med konstant acceleration och når sin maxfart vmax efter 2.5 s, sedan behåller han denna fram till målgång efter 100 m med sluttiden 10.4 s. Sök sprinterns maxfart vmax! Lösningsförslag: Typisk drillövning av vt-diagram! Konstant acceleration innebär att farten ökar linjärt upp till maxfart. Vi har alltså situationen vt ms vmax 2.5 10.4 t v 0 Den tillryggalagda sträckan är lika med arean under kurvan eftersom s t s t. Detta möblerar tillräckligt med villkor för att bestämma vmax 2.5 ekv sacc svmax 100, sacc 0 sacc svmax 100, sacc 1.25 vmax , svmax vmax 10.4 t t, svmax 2.5 7.9 vmax Nu är det bara att lösa ut vmax, de två delsträckorna får vi på köpet. 2.5 vmax t 34 Mekanik-Dynamik och Mathematica HH/ITN/BN Solve ekv sacc 13.6612, svmax 86.3388, vmax 10.929 Exempel: Under tiden efter andra världskriget gjorde sig överste John P. Stapp 1910–1999 berömd för att utsätta sin egen kropp för flera hundra makabra experiment i syfte att utreda vilka påkänningar människokroppen tål. Mest känt är då han 1954 på en släde bromsades in från över 1000 km h till stillastående på 1.4 s Efter detta fick han bland annat problem med synen eftersom blodkärlen i näthinnan blev skadade. Förloppet övervakades naturligtvis och hastigheten under inbromsningen kunde med god noggrannhet beskrivas med uttrycket v t 170 1.4 t 3 2 . Bestäm hans maximala retardation samt bromssträckan. Lösningsförslag: Typexempel på förståelse av vt-kurvor. v t : 170 1.4 3 2 t Plot v t , t, 0, 1.4 , PlotStyle Red, AxesLabel "t s ", "v t m s " vt ms 250 200 150 100 50 0.2 0.6 0.8 1.0 t s 1.2 1.4 En direkt kalkyl med hjälp av definitionerna a t v t a t 0.4 och v t x . t D v t ,t 255 1.4 t Med uppenbart störst retardation i ms2 då t 0. a 0 301.72 Inbromsningssträckan i meter fås efter separation och integration av v t x 0 x 1.4 v 0 x 0 x t t x 170 32 1 1.4 t 32 1 1.4 0 340 5 1.45 2 x , t ty x = arean under vt-kurvan. Vi får 157.7 m. 1.4 x v t t 0 157.699 Exempel: Två vagnar med massorna m1 2 kg respektive m2 3 kg är i vila på rak räls. En fjäder med fjäderkonstanten k 10 N m förbinder dem med varandra. Plötsligt ges vagnarna hastigheterna 1 m s. Sök deras läge v1 2 m s respektive v2 och hastighet som funktion av tiden. k m1 x1(t) m2 x2(t) HH/ITN/BN Mekanik-Dynamik och Mathematica 35 Lösningsförslag: Låt x1 t och x2 t vara vagnarnas läge som funktion av tiden. Efter friläggning får vi direkt med Newton .. 2 x1 t 10 x2 t x1 t ODE .. BVP 10 x2 t x1 t 3 x2 t x1 0 0, x1 0 2 x2 0 0, x2 0 1 BV Det är bara att muppa rakt in i Mathematica. x12 DSolve 2 x1 '' t 10 x2 t x1 t , 3 x2 '' t 10 x2 t x1 t , x1 0 0, x1 ' 0 2, x2 0 0, x2 ' 0 1 , x1 t , x2 t , t First FullSimplify 1 x1 t 5t 5t 9 1 3 sin 25 , x2 t 5t 5t 6 3 sin 25 3 3 Nu kan vi inspektera skådespelet! Läge, hastighet och acceleration. Första vagnen röd. Plot Evaluate PlotStyle x1 t , x2 t . x12 , t, 0, 10 , Red, Blue , AxesLabel "t s ", "x1 t ,x2 t m " x1 t ,x2 t m 2.5 2.0 1.5 1.0 0.5 2 PlotEvaluate D PlotStyle 4 6 8 x1 t , x2 t 10 t s . x12, t Red, Blue , AxesLabel , t, 0, 10 , "t s ", "x1 t ,x2 t m s " x1 t ,x2 t m s 2.0 1.5 1.0 0.5 t s 2 0.5 4 6 8 10 1.0 1.5 PlotEvaluate D PlotStyle .. x1 t , x2 t . x12, t, 2 Red, Blue , AxesLabel "t .. x1 t ,x2 t m s2 4 2 2 2 4 4 6 8 10 t s , t, 0, 10 , .. .. s ", " x 1 t , x 2 t m s2 " 36 Mekanik-Dynamik och Mathematica HH/ITN/BN Exempel: En basebollspelare kastar iväg bollen enligt figur. Sök bollens högsta höjd, den så kallade stighöjden, kastlängd samt tidpunkterna för dessa två tillfällen Rita upp bollbanan Lösningsförslag: Typexempel! Ta alltid tillfället i akt att göra en symbolisk lösning eftersom den är utbildande jämfört med ett själslöst numeriskt exempel hämtat "ur luften". Ta hjälp av Newton och ställ upp rörelseekvationerna i x- och y-riktningarna. Applicera begynnelsevärden och lös begynnelsevärdesproblemet. xy DSolve m x '' t m y '' t 1 x t t v0 cos Α , y t 0, x 0 0, x ' 0 v0 Cos Α , m g, y 0 h, y ' 0 v0 Sin Α g t2 t Solve y ' t v0 sin Α ,t First 2 t v0 sin Α 2 Vid högsta punkten på banan är y tH 2h , x t ,y t 0. Detta bestämmer önskad information vad gäller stighöjd. 0 . D xy, t , t First g xy . tH x v0 sin Α v20 sin Α cos Α g g ,y v0 sin Α g 1 v20 sin2 Α 2h g 2 Med numeriska data data v0 30, Α 30 , h 2, g 9.81 ; xy . tH . data x 1.52905 39.7259, y 1.52905 13.4679 På ungefär samma sätt vaskar vi fram tillståndet vid nedslag. tL Solve y t 0 . xy, t v0 sin Α 2gh t v20 sin2 Α Simplify 2gh , t v20 sin2 Α g v0 sin Α g tL . data t 0.127978 , t 3.18608 Första lösningen motsvarar kastparabelns "historiska" startpunkt på marken bakom kastaren. xy . tL 2 2gh x Simplify v20 sin2 Α v0 cos Α v0 sin Α 2gh g Med numeriska data har vi exemplets kastlängd och restid. x 3.18608 . data 82.7768, y 3.18608 v0 sin Α 2gh ,y g xy . tL 2 v20 sin2 Α 0. Slutligen två små bilder över bollens resa genom vädret. v20 sin2 Α g v0 sin Α 0 HH/ITN/BN Mekanik-Dynamik och Mathematica ParametricPlot 37 x t ,y t . xy . data, t, 0, t . tL 2 . data , PlotStyle Orange, AxesLabel "x m ", "y m " y m 14 12 10 8 6 4 20 40 Plot Evaluate PlotStyle 60 80 x m x t ,y t . xy . data , t, 0, t . tL 2 . data , Red, Blue , AxesLabel "t s ", "x t ,y t m " x t ,y t m 80 60 40 20 0.5 1.0 1.5 2.0 2.5 3.0 t s Exempel: En basketsituation är uppriggad enligt figur. Med given elevationsvinkeln Θ 60 måste bollen ges en speciell utgångsfart v0 för att träffa korgen. Sök denna samt restiden fram till korgen. Försumma luftmotståndet. Lösningsförslag: Placera ett naturligt xy-koordinatsystem enligt firgur. Vid godtycklig tidpunkt är det endast tyngdkraften som verkar på bollen. Formulera rörelseekvationerna, applicera begynnelsevärden och lös begynnelsevärdesproblemet xÅy x t DSolve m x '' t 0, m y '' t m g, x 0 x0 , x ' 0 v0 Cos Θ , y 0 y0 , y ' 0 v0 Sin Θ . x0 0, y0 2.1, Θ 60 , g 9.81 , x t , y t , t 0.5 t v0 , y t 4.905 t 2 0.866025 t v0 First 2.1 Svaret på de båda delfrågorna ges av nedslagsplatsen. Vi har två ekvationer och två obekanta v0 och restiden t. vÅt t Solve x t 1.1086, v0 4, y t 7.21632 , t 3 . xÅy 1.1086, v0 7.21632 Sista lösningen är den enda vi befattar oss med, så det får bli svaret på uppgiften. Avslutningsvis en liten reseberättelse 38 Mekanik-Dynamik och Mathematica HH/ITN/BN bana Last xÅy . vÅt 2, 2 ; ParametricPlot bana, t, 0, 1.1 , PlotRange 0, 4.5 , AxesLabel x, y , PlotStyle Blue, Dashing 0.025 , AspectRatio Automatic, Epilog Green, Thickness 0.03 , Line 3.9, 3 , 4.1, 3 , Red, PointSize 0.05 , Point bana . t 0.8 y 4 3 2 1 x 0 1 2 3 4 Exempel: En ögonblicksbild av en höjdhoppare precis vid upphoppet kan beskådas i figuren. Sök v0 och Θ så att hopparens tyngdpunkt G precis klarar höjden. Lösningsförslag: Dessa eviga kastparabler. Nu är det dax igen. Placera ett koordinatsystem med origo där G befinner sig enligt figur. Ta sedan hjälp av Newton och ställ upp rörelseekvationerna i x- och y-riktningarna. Applicera begynnelsevärden och lös begynnelsevärdesproblemet. xy DSolve m x '' t m y '' t 0, x 0 0, x ' 0 v0 Cos Θ , m g, y 0 0, y ' 0 v0 Sin Θ 1 x t t v0 cos Θ , y t 2 2 t v0 sin Θ , x t ,y t ,t First g t2 När klockan är t har vi nått högsta punkten på banan och då är y 0. Detta möblerar ett ekvationssystem vars lösning är den önskade informationen. Vi har att uppfylla tre villkor (ekvationer) med tre obekanta, v0 , Θ och stigtiden t som bonus. Solve v0 v0 x t 5.04228, Θ 5.04228, Θ 1, y t 1.13005, t 2.01155, t 1.06, y ' t 0.464872 , v0 0.464872 , v0 0 . xy . D xy, t 5.04228, Θ 5.04228, Θ .g 9.81, v0 , Θ, t 2.01155, t 0.464872 , 1.13005, t 0.464872 Eftersom vi varken befattar oss med negativa farter, tider eller vinklar är det bara den sista lösningen som duger. I höjdhopp behöver inte tyngdpunkten passera ovanför ribban. Detta var svårt att uppnå på "Hedenhös" tid då man använde den så kallade saxstilen. Lite bättre blev det med dykstilen och dess utveckling med så kallat "hängande ben", men det verkliga genombrottet kom under OS 1968 i Mexico City då amerikanen Dick Fosbury (1947-), på bild nedan, chockade friidrottsvärlden med att visa upp en helt ny hoppstil, "the Fosbury Flop", som han vann guld med på 2.24. För den riktigt vige erbjuder denna stil just en rejäl sänkning av tyngdpunkten under det att kroppen krånglar sig över ribban. Vi som känner till mekanikens energilagar vet att detta betyder ökad hopphöjd för given utgångsfart. HH/ITN/BN Mekanik-Dynamik och Mathematica 39 Exempel: En fisk med massan m slutar simma vid hastigheten v0 och glider horisontellt genom vattnet. Bestäm fiskens hastighet som funktion av läget samt hastighet och läge som funktion av tiden om friktionskraften mot vattnet är proportionell mot hastigheten med proportionalitetskonstanten k. Hur långt rör den sig ? När stannar den? .. Lösningsförslag: Vi känner igen Newton, mx F. Vi börjar med att rulla ut måttbandet, det vill säga x-axeln, precis då den slutar simma. Den enda kraft som verkar i färdriktningen är då motståndet från vattnet kx. Vi får då .. BVP mx x0 kx 0 x0 v0 ODE BV .. I första fallet är vi primärt intresserade av hastighet som funktion av läget x x . Vi gör därför omskrivningen x x t KR x x x t x . x x Därmed är vi över i en separabel (ODE) som vi integrerar direkt med (BV) som gränser x x mx x v0 kx x 0 x k m x x v0 x x k x m 0 x k x. m v0 Så hastigheten avtar linjärt med läget eller sträckan den glidit. För att få hastighet och läge som funktion av tiden kan vi antingen ge oss på den nyss erhållna x .. k x m v0 eller den ursprungliga mx kx. Notera att i båda fallen är x x t . Vi provar båda vägarna. Den förstnämda är både separabel och linjär. Eftersom vi har haft separationsångest ett tag så väljer vi att öva lite på en linjär av första ordningen k x k x m v0 k x m x v0 k m Slutligen BV x 0 t k mv0 k x k mv0 k 0: 0 IF C1 m m m k t t m t C1 0 k x m mv0 k t k x m k mv0 k x C1 C1 m t v0 Separabel t lösningen till BVP är x t mv0 1 k k m t Hastigheten får vi genom att derivera x t x t k mv0 1 k m t mv0 0 k k k m m t k v0 m t Så hastigheten avtar exponentiellt med tiden den glidit. Man ska alltid hålla vad man lovat varför det nu är dax för den andra vägen med ursprungliga (ODE) som är en linjär andra ordningens homogen med konstanta koefficienter. .. mx karakteristiska ekvationen och dess två rötter r1,2 .. kx k x m x k 2m k 0 2 2m r1 0 r2 0 k k m Fall 1 C2 m k Eftersom ODE är homogen så är x t xh t xp t C1 xh t C1 C2 m t t 0 k Slutligen BV x0 C1 0 x0 v0 C2 k m C2 mv0 k Glidsträckan blir ett gränsvärde x t m k då t m 0 0 0 v0 C1 C2 mv0 k mv0 k återigen lösningen till BVP x t mv0 1 k DSolve m x '' t k x' t , x 0 0, x ' 0 v0 , x t , t kt m v0 1 m x t k Om vi som exempel väljer m k v0 t . Så resan är begränsad i rummet men inte i tiden. Den stackarn glider för gott! Avslutningsvis kollar vi om Mathematica delar vår åsikt. xAvt k m 1 får vi en kvalitativ överblick av resan. First Simplify 40 Mekanik-Dynamik och Mathematica PlotEvaluate xAvt 1, 2 , D xAvt 1, 2 , t PlotStyle Red, Blue , PlotRange . m 1, k All, AxesLabel HH/ITN/BN 1, v0 1 , t, 0, 5 , "t", "x t , x t " xt,xt 1.0 0.8 0.6 0.4 0.2 t 1 2 3 4 5 Exempel: Ett så kallat bungy jump är väl bekant för de flesta. Från en hoppställning högt ovan mark kastar man sig handlöst ut. Livlinan utgörs av en gummisnodd vars ena ände är fäst vid uthoppsplatsen och den andra surrad runt vristerna på offret. Välj lite numeriska data och gör en simulering av ett hopp .. Lösningsförslag: Om vi placerar en y-axel med origo vid marken och pekande uppåt har vi återigen Newton, m y F. Krafterna som verkar är tyngdkraften och den från gummisnodden. Eftersom en gummisnodd, för enkelhets skull antar vi, fungerar ungefär som en fjäder frånsett att den inte kan ta trycklaster, är den lite kinkig att simulera och lösa analytiskt. Men med storsläggan NDSolve går det både snabbt och enkelt. För att offret inte ska svänga upp och ned i all oändlighet, vilket är fallet med en fjäder, lägger vi in lite dämpning i snodden. Med fjäderkonstanten k, dämpning c och den naturliga längden L gäller om hoppställningen har höjden H att kraften i snodden blir .. kH 0 Fgs y L c y om H om H y y L L 0 0 och BVP my mg Fgs y0 H y0 0 ODE BV Nu gäller det bara att översätta den styckvis definierade funktionen Fgs till Mathematica, vilket lätt görs med funktionen Piecewise eller If. Vidare kräver NDSolve att resan görs helt numeriskt, så välj t.ex. m c 50 Ns/m, så har vi direkt en studie över de T första sekunderna T 75 kg, H 50 m, L 15; Fgs yAvt k H y t L 0 c y' t H H y t y t L L 0 ; 0 m y '' t m g Fgs , y 0 H, y ' 0 0 . m 75, g 9.81, H 50, L y t , t, 0, T First; NDSolve 30, k 125, c 50 , ODE BV Data En bild över spektaklet piggar alltid upp PlotEvaluate y t PlotStyle . yAvt, D y t . yAvt, t Red, Blue , AxesLabel "t , t, 0, T , s ", "y t m , y t yt m,yt ms 40 20 2 4 6 8 10 12 14 t s 20 Man ser tydligt på hastighetsgrafen när gummisnodden börjar bromsa in det fria fallet. m s " 30 m, k 125 N/m och HH/ITN/BN Mekanik-Dynamik och Mathematica 41 Exempel: En fallskärmshoppare, vars vikt är 75 kg, hoppar ut från en helikopter 4 km ovan markytan. Antag att luftmotståndet är proportionellt mot hopparens fart. Proportionalitetsfaktorn är 15 kg s när fallskärmen är outlöst och 105 kg s när fallskärmen är utvecklad. Antag att fallskärmen utvecklas 1 min efter uthoppet från helikoptern. Sök läge och hastighet under resan samt bestäm hur lång tid hoppet tar. .. Lösningsförslag: Om vi placerar en y-axel med origo vid marken och pekande uppåt har vi återigen Newton, my F. Krafterna som verkar är tyngdkraften och den från fallskärmen. Vi har att ta ställning till två stycken (BVP). Ett under den första minuten med kända (BV), sedan ett annat under resten av nedfärden. Denna måste matchas med (BV) = (RV) vad gäller tid, läge och hastighet så att den sammanfogas med lösningen under den första minuten. Eller lite mera precist .. BVP1 my1 y1 0 mg k1 y1 H y1 0 0 .. ODE BV för t 0, 60 och BVP2 my2 y2 60 y2 60 mg k2 y2 y1 60 y1 60 ODE BV för t 60, När man så äntligen trött kommer i mål gäller det att mata Plot med rätt lösning i rätt intervall . Räddaren heter naturligtvis NDSolve. Vi löser (BVP1) för hela resan t 0, och byter helt enkelt proportionalitetskonstant k1 k2 vid 60 s! Smidigt . Nu är det bara att sätta igång! Kör ordentligt i backen så tar vi reda på efteråt hur lång tid hoppet tog! yAvt y t NDSolve m y '' t m g If t Τ, k1 , k2 y ' t , y 0 H, y ' 0 0 . m 75, g 9.81, H 4000, Τ 60, k1 y t , t, 0, 500 First Domain: 0. 500. InterpolatingFunction Output: scalar 15, k2 105 , ODE BV Data t Hopptid Thopp t . FindRoot y t 0 . yAvt, t, 200 241.56 Några bilder över spektaklet piggar alltid upp Plot y t . yAvt, t, 0, Thopp , PlotRange AxesLabel "t PlotEvaluate D y t PlotStyle s ", "y t . yAvt, t Blue, AxesLabel All, PlotStyle Red, m " , , t, 0, Thopp , PlotRange All, "t s ", "y t yt m m s " yt ms 4000 50 150 200 t s 10 3000 100 , 20 2000 30 1000 40 50 100 150 200 t s 50 Man ser tydligt på hastighetsgrafen när fallskärmen börjar bromsa in det fria fallet. Även de två gränshastigheterna är lätta att identifiera, det vill säga då vi inte har någon acceleration längre. Dessa blir utan respektive med fallskärm D y t . yAvt, t 49.0492, 7.00714 .t 55, 200 42 Mekanik-Dynamik och Mathematica HH/ITN/BN Exempel: En vattenraket skjuts iväg i geografin enligt figur. Bestäm sträckan R längs backen upp till nedslagsplatsen samt restiden. Sök även nedslagshastigheten. Lösningsförslag: Placera ett naturligt xy-koordinatsystem vid A. Vi får då rörelseekvationerna med begynnelsevärden. Lös ut x t och y t . xÅy DSolve m x '' t 0, m y '' t m g, x 0 0, x ' 0 v0 Cos Θ , y 0 , x t ,y t ,t First 1 x t t v0 cos Θ , y t 2 2 t v0 sin Θ 0, y ' 0 v0 Sin Θ g t2 Svaret på de båda första delfrågorna ges av nedslagsplatsen. RÅt R Solve x t 0, t R Cos Α , y t R Sin Α 2 v20 sec2 Α cos Θ sin Α 0 , R Θ . xÅy, R, t 2 v0 sec Α sin Α ,t g Θ g Första lösningen är startpunkten A och den andra den sökta nedslagsplatsen B. Nu är Θ upp oss med ett numeriskt exempel. data v0 50.0, Θ 60 , Α Simplify 30 , g Α så minustecknen efter 9.81 ; xÅy . RÅt, RÅt . data x0 R 0, y 0 0 0, t 0 x 5.88532 R 147.133, y 5.88532 84.9473 169.895, t 5.88532 Nedslagshastighet och fart. Som biprodukt får vi även en kontroll av att utgångshastighet och fart är de föreskrivna. DxÅy, t . RÅt . data x 0 x' t 25., y 0 2 43.3013 , x 5.88532 y' t 2 25., y 5.88532 14.4338 . DxÅy, t . RÅt . data 50., 28.8675 ParametricPlot Evaluate x t , y t . xy . data , t, 0, t . Last Rocht PlotStyle Red, AspectRatio Automatic, PlotRange All, AxesLabel "x", "y" , Epilog Green, Thickness 0.02 , Line 0, 0 , 200 Cos Α , 200 Sin Α . data y 100 80 60 40 20 50 100 150 x . data , är ok! Vi piggar HH/ITN/BN Mekanik-Dynamik och Mathematica 43 Exempel: För ett bilbälte är den inbromsande kraften proportionell mot kvadratroten ur bröstkorgens läge under krockförloppet. För att det inte ska bli sista vilan för en passagerare med massan m krävs att denne kommer till vila efter sträckan a. Bestäm proportionalitetskonstanten k så att detta uppfylls vid en krock med hastigheten u .. Lösningsförslag: Standard Newton mx .. x x x . x F. Eftersom vi är intresserade av hastighet som funktion av väg gör vi omskrivningen Då k är okänd kan vi göra ett litet trick och låta k x vara en okänd konstant funktion. Vi har då möjlighet att få med randvillkoret x a vAvx 0 också, med resultat att DSolve gör hela jobbet direkt. Dessutom får vi ut x x redo för uppritning. DSolvem v x v ' x v 0 u, v a v x ,k x 3 m u2 k x 4 a3 2 k x 0, 0, , x u2 a3 2 ,v x x , k' x x3 2 3 m u2 , k x a3 2 4 a3 2 u2 a3 2 ,v x a3 2 x3 2 Numersikt exempel. Plot Evaluate v x . vAvx 2 . m 70, u 20, a 0.5 , x, 0, 0.5 , PlotStyle Red, PlotLabel "Hastighet som funktion av x", AxesLabel "x", "v x " Hastighet som funktion av x vx 20 15 10 5 0.1 0.2 0.3 0.4 x 0.5 Annars går det naturligtvis lika bra att låta k vara en konstant svar DSolvem v x v ' x 3 m u2 4 k x3 2 k x ,v 0 3 m u2 4 k x3 2 m m , v x v x 3 3 och sedan bestämma den ur randvillkoret x a Solve v x 0 . svar 2 .x 0 a, k 3 m u2 k 4 a3 2 Exempel: Kraftöverföringen i en enkel ångmaskin kan studeras i vidstående figur. Låt nu vevaxeln gå med konstant varvtal 2Π rad s. Bestäm sedan hastigheter .. och accelerationer r, r, Θ, Θ. Rita ut dem u, v x , x 44 Mekanik-Dynamik och Mathematica HH/ITN/BN Lösningsförslag: Här passar det alldeles utmärkt med en liten vektorbetraktelse. Inför ett koordinatsystem i O med x-axeln åt höger och y-axeln rakt upp som vanligt. Eftersom vi har en typisk 2D modell nöjer vi oss med en sådan studie och får om L R CA och t data L OA . Vinkeln Θ får vi lätt ur sinussatsen. 300 90 ,R ; 1000 1000 L, 0 R Cos Π Β t , Sin Π R Sin Β t ArcSin . data; L t Θ t Nu är det bara att rita. Vi börjar med PlotEvaluate t , Β t . data; t och dess komponenter. t , Β t , 0, 2 Π , PlotStyle Green, Blue, Red , AxesLabel 3Π , 2 Π, Automatic , Π, 2 2 "Β rad ", "rx ,ry , m " , Π Ticks rx ,r y , m 0.4 0.3 0.2 0.1 Β rad Π Π 2 0.1 PlotEvaluate PlotStyle 2Π 2 ' t . Β' t 2 Π , Β t , 0, 2 Π , Green, Blue, Red , AxesLabel "Β rad ", "rx ,ry , m s ", Π 3Π , 2 Π, Automatic , Π, 2 2 Ticks rx ,r y , ' t , 3Π ms 0.6 0.4 0.2 Β rad Π 0.2 Π 2 3Π 2Π 2 0.4 0.6 PlotEvaluate PlotStyle .. .. .. '' t . Β' t Green, Blue, Red , AxesLabel 2 Π, Β '' t "Β 0 rad ", , Β t , 0, 2 Π , .. .. .. " r x, r y, m s2 ", Π 3Π , Π, , 2 Π, Automatic 2 2 Ticks rx ,r y , '' t , m s2 3 2 1 Β rad Π 1 Π 2 3Π 2 2Π 2 3 PlotEvaluate PlotStyle Ticks Θ t , Θ' t . Β' t Blue, Red , AxesLabel 2 Π, Β '' t "Β Π 3Π , 2 Π, Automatic , Π, 2 2 0 rad ", "Θ , Β t , 0, 2 Π , rad ,Θ rad s ", OC , HH/ITN/BN Mekanik-Dynamik och Mathematica 45 Θ rad ,Θ rad s 2 1 Β rad Π 3Π Π 2Π 2 2 1 2 PlotΘ '' t . Β' t PlotStyle 2 Π, Β '' t Orange, AxesLabel 0 , Β t , 0, 2 Π , "Β rad ", "Θ rad s2 ", Π 3Π , 2 Π, Automatic , Π, 2 2 Ticks Θ rad s2 10 5 Β rad Π 3Π Π 2 2 5 2Π 10 Exempel: En karusell har banprofilen t 5 cos 3t , 5 cos t 3 sin 2t , 5 sin t förhållande till övriga mått. a) Bestäm banans längd. b) Bestäm hastighet som funtion av tiden samt största och minsta fart. c) Bestäm acceleration som funtion av tiden samt beloppet av denna. m, t 0, 2Π s. Låt vagnarna vara små i Lösningsförslag: Banprofilen, det vill säga läget av en vagn som funktion av tiden. t : 5 Cos 3 t , 3 Plot Evaluate AxesLabel Cos t Sin 2 t , Sin t ; t , t, 0, 2 Π , PlotStyle Red, Blue, Green , "t s ", "x t ,y t ,z t m " x t ,y t ,z t m 20 10 t s 1 2 3 4 5 6 10 20 a. Först den lilla båglängden vid tidpunkten t. ds ' t 5 36 sin2 3 t 4 cos2 t cos t 3 4 cos 2 t cos 3 t 2 2 Sedan hela banlängden. 2Π N ds t 0 143.958 b. Hastighet är tidsderivatan av läget, ' t , det vill säga tangentvektorn. Farten är beloppet av hastigheten, v PlotEvaluate AxesLabel ' t , "t ' t , t, 0, 2 Π , PlotStyle s ", "x t ,y t ,z t ,v t m s " Red, Blue, Green, Orange , . 46 Mekanik-Dynamik och Mathematica HH/ITN/BN x t ,y t ,z t ,v t m s 40 30 20 10 t s 1 10 2 3 4 5 6 20 30 Högsta och lägsta fart samt hur mycket klockan är då det inträffar. FindMinimum 6.5423 40.3113 t t ' t , t, 2 , FindMaximum ' t , t, 0 2.19013 0. En liten rundtur längs en färgglad bana piggar alltid upp. Ju rödare det är ju fortare går det ParametricPlot3DEvaluate ColorFunction t . data , t, 0, 2 Π , PlotStyle .t Function x, y, z, u , Hue Thickness 0.03 , 2 Π u 40 c. Accelerationen är andraderivatan av läget med avseende på tiden, '' t , samt dess belopp. PlotEvaluate '' t , AxesLabel .. .. '' t .. , t, 0, 2 Π , PlotStyle .. .. s ", " x t , y t , z t ,a t "t Red, Blue, Green, Orange , m s2 " .. x t ,y t ,z t ,a t m s2 50 1 2 3 4 5 6 t s 50 Högsta och lägsta värde på beloppet av accelerationen samt hur mycket klockan är då det inträffar. Eftersom en indikation på vilken påkänning en resenär upplever under resan. FindMinimum 15.335 84.0751 t t '' t 1.49337 0.766509 , t, 1 , FindMaximum '' t , t, 1 m ger detta HH/ITN/BN Mekanik-Dynamik och Mathematica 47 Exempel: En bil kör i berg och dalbana. Farten i punkt A är 100 km h och bromsas med konstant acceleration ned till 50 km h i punkt C. Bilens tyndpunktshöjd är försumbar jämfört med vägens krökningsradie. Bestäm nu a vägens krökningsradie i punkten A om passagerarna upplever att storleken av den totala accelerationen är 3 ms2 . b storleken av totala accelerationen i punkt C. Lösningsförslag: Vi börjar med att bestämma tangentiella accelerationen at under inbromsningen från A till C. Med den populära .. omskrivningen x x x x at får vi 50 60 60 3.6 acc Solve 100 x x at x First 0 3.6 2.41127 at Nu är vi rustade för att besvara frågor! Först a med hjälp av "cirkelrörelseformler". NSolvea2tot v a2t 27.7778, atot v2 a2n , an 3., an 100 , atot ,v 1.78488, r 3 . acc 3.6 r 432.301 , v 27.7778, atot 3., an 1.78488, r 432.301 varav den negativa förkastas. Med samma formelsida uppslagen får vi även svaret på delfråga b. Solvea2tot a2t v2 a2n , an 50 150., v 1.28601, atot 13.8889, an 150 . acc 3.6 r r ,r ,v 2.73277 , r 150., v 13.8889, an 1.28601, atot 2.73277 Här duger bara den sista lösningen eftersom vi inte befattar oss med vektorer som har negativ längd. Exempel: Tarzans lian har en säkerhetsmarginal på 5 då han hänger rakt ned i den, det vill säga den håller för 5 st Tarzan. Antag nu att Tarzan svingar sig från en klippa på samma höjd som den gren lianen sitter fast i. Hur många LillTarzan kan han ta med sig utan att lianen brister om en LillTarzan väger 1 4 av vad Tarzan väger? Lösningsförslag: Vi antar att lianen har längden r och att han tar n st LillTarzan med sig. Energibetraktelse T1 farten v i lägsta punkten där ju påkänningen är som störst. 1 ekv n mLT 02 mT 1 mT n mLT g r 2 1 v n mLT v2 mT n mLT g 0 2 g r n mLT vmax mT V1 mT 2 v2 n mLT Solve ekv 2 g mT Simplify, v r , v Med centripetalaccelerationen an 2 v2 r g r och spännkraften S i lianen har vi i lägsta punkten rörelseekvationen i radiell led T2 V2 ger 48 Mekanik-Dynamik och Mathematica HH/ITN/BN v2 ekv mT S n mLT mT n mLT g . vmax 2 r 2 g n mLT mT S g n mLT mT Varav spännkraften i lianen. Slian S Solve ekv, S Simplify First mT 3 g n mLT Med given vikt på en LillTarzan bestäms antalet n av kravet på säker resa. 1 ekv 5 mT g mT S . Slian . mLT 4 5 g mT n mT 3g mT 4 Solve ekv, n n N 2.66667 Han kan ta med sig högst 2 st LillTarzan. Exempel: På ett nöjesfält finns en märklig bollbana. Antag att en pojke står på avståndet x 3R. Vilken fart ska han ge bollen för att den ska komma tillbaka ? Betrakta bollen som en partikel Vilket är det minsta avstånd x som pojken kan stå på för att bollen ska kunna komma tillbaka på detta sätt? Lösningsförslag: Inför ett koordinatsystem med origo i B där x-axeln är riktad mot pojken och y-axeln uppåt. Låt bollens hastighet vid A ha farten u och riktad i negativ x-riktning och dess hastighet vid C ha farten v och riktad i positiv x-riktning. Lämna sedan över till Solve att reda ut resan med energibetraktelse B C samt kastparabel C A med restid T. m u2 uvT m v2 mg2R Solve 2 3 g 2 g g R R ,v 2 5 g R g g R 3 g , g R 2 R ,T 2 R R ,T 2 ,v 2 2 g 2 , u 2 3 ,v 5 R ,T 2 v T, u, v, T , u g 3 g 2 R, 3 R R R ,T 2 5 u 2 R ,v g T2 2 2 5 u 1 ,0 g Eftersom vi varken befattar oss med negativ restid T eller negativa farter u, v (belopp av hastighet!) är det bara den sista lösningen som duger. Det är svaret till första delfrågan. Även sista delfrågan passar Solve som handsken. xMin Solve m v2min m g, xmin vmin T, xmin , vmin R xmin g Last xMin xmin 2R R T, xmin . Last uvT g R T HH/ITN/BN Mekanik-Dynamik och Mathematica Exempel: Lilla Lisa vill åka loopen på Liseberg. Men innan hon vågar sig upp vill hon gärna kontrollera hur högt startpunkten är i förhållande till själva loopen. Hon uppskattar både loopens radie r och startpunktens höjd. Nu vill hon ha hjälp med en formel så att hon kan jämföra den uppskattade höjden med den erforderliga höjden h för att hon ska komma oskadd ur äventyret 49 h r r Lösningsförslag: Newton i normalriktning och energiekvationen. Det blir kritiskt i övre läge om kontaktkraften N ekv m v2min 1 mg m v2min m g 2 r, N 0 2 r m v2min N, m g h 0. gm N, g h m 2gmr r m v2min ,N 0 2 Med den sökta lösningen på h, och hastigheten i övre läget. Solve ekv, h, vmin , N 5r h 2 5r g , vmin r ,N 0, h 2 g , vmin r ,N 0 Exempel: En friktionsfri rak glidbana fungerar som utskjutningsramp för en låda, se fig. Lådan släpps vid banans övre ände och lämnar den vid den nedre.Vilken vinkel Α skall banan luta för att hastighetens horisontella komposant vh ska blir så stor som möjligt när lådan lämnar banan ? α vh v Lösningsförslag: Anta att banans längd är L. Energibetraktelse, med nollnivå vid banans lägsta punkt, ger nu direkt lådans fart v längs banan precis då den lämnar banan. 1 Solve m 02 1 m g L Sin Α 2 v 2 m v2 m g 0, v 2 g L sin Α , v 2 g L sin Α Den eftersökta horisontella komposanten. vh v Cos Α 2 g L . 2 sin Α cos Α Π Plotvh . g 1, L 1 , Α, 0, , PlotStyle Orange, AxesLabel 3 vh 0.8 0.6 0.4 0.2 Α 0.2 0.4 0.6 0.8 1.0 Numeriskt extremvärde ges enklast av FindMaximum eller av NMaximize. Π FindMaximumvh . g 1, L 1 , Α, 0, 6 "Α", "vh " 50 Mekanik-Dynamik och Mathematica 0.877383, Α HH/ITN/BN 0.61548 Π NMaximizevh . g 1, L 1 ,0 , Α Α 4 0.877383, Α 0.61548 Annars kan man ju söka nollställe till derivatan. ekv D vh , Α g L cos2 Α 0 3 2 2 g L sin 2 Α alfa Solve ekv, Α Α cos 2 1 , Α 2 1 cos , Α 3 alfa Α 0 sin Α cos 1 3 2 , Α 2 1 cos 3 3 N 2.52611 , Α 2.52611 , Α 0.61548 , Α 0.61548 Här är det bara sista lösningen som duger. 180.0 . Last alfa Α Π 35.2644 För hand (varför då?) är det lättare att studera v2h . Den har samma extremvärde på Α. Α sin Α cos2 Α 0 cos Α cos2 Α sin Α 2 cos Α cos Α cos Α cos2 Α 1 ArcSin 2 sin2 Α 0 cos2 Α 0 Α sin Α Π 2 2 sin2 Α 0 vh 0 0, ointressant 1 3 sin2 Α 0 sin Α 1 3 180.0 3 Π 35.2644 Exempel: En friktionsfritt lagrad trumma är uppriggad enligt figur. Jämför utvecklingen över tiden för de två olika metoderna att sätta den i rotation. Bestäm även kraften i snöret Tröghetsradien för trumman är 375 mm och dess vikt 100 kg. Lösningsförslag: Inledningsvis tar vi och samlar ihop lite numeriska data. data M 20, m 100, r 0.25, k 0.375, g 9.81 ; Först metoden med vikt, det vill säga a. Klipp av snöret och kalla snittkraften S. Rörelseekvationer för trumma och vikt samt kopplingsvillkor och begynnelsevärden riggar sedan upp begynnelsevärdesproblemet. Tröghetsmomentet för trumman är J mk 2 där k är tröghetsradien. HH/ITN/BN Mekanik-Dynamik och Mathematica bvpa SolveJ '' t r S, M x '' t '' t , x '' t , S, J t k2 m M r2 ,x t S, x '' t r '' t , J m k2 , First g M r2 gMr Mg 51 k2 m g k2 m M M r2 ,S k2 m M r2 k 2 m ,J Här kan vi direkt avläsa de konstanta accelerationerna samt den konstanta kraften i snöret. Nu över till utvecklingen över tiden. x a DSolve bvpa 1, 2 . Rule Equal, x 0 t ,x t ,t First g M r t2 t 2 k2 m 0, x ' 0 0, 0 0, ' 0 0 , g M r2 t 2 ,xt M r2 2 k2 m M r2 Sedan metoden med konstant kraft, det vill säga b. bvpb SolveJ '' t r S, S 196.2, J k 2 m 196.2 r t k2 m x b ,S DSolve bvpb 1 . Rule 20 9.81, J Equal, 0 m k2 , 0, ' 0 '' t , S, J 0 , t ,t 3.488, S 196.2, J First First 98.1 r t2 t k2 m Inspektion av .. avgör tävlingen till fördel för metod b. bvpa, bvpb t . data 3.20327, x t 0.800816, S 180.184, J 14.0625 , t 14.0625 Avslutningsvis en bild över skådespelet. Metod a med tunn och b med fet linje. c Red, Green, Blue, Orange ; PlotEvaluate Join x t , t , x' t , PlotStyle AxesLabel xa t , a ' t . D x a, t . x a, r t , t ,r ' t , ' t . D x b, t .x b . data , t, 0, 5 , Join c, Thick, & c , "t", None , PlotLabel "xa t , a t ,xa t , a t ,xb t , b t ,xb t , t ,xa t , a t ,xb t , b t ,xb t , b b t " t 30 25 20 15 10 5 1 2 3 4 5 t Exempel: När skridskoprinsessan Miss Inertia gör en piruett håller hon inledningsvis armarna enligt den övre figuren och snurrar med vinkelhastigheten 5 rad s. Bestäm vinkelhastigheten efter det att hon har fällt in armarna enligt den undre figuren. Lite data, Huvud: Klot med radie 0.1 m och vikten 4 kg, Bål: Cylinder med radie 0.2 m, höjd 0.6 m och vikten 15 kg, Arm: Cylinder med radie 0.05 m, höjd 0.7 m och vikten 5 kg, Ben: Kon med basradie 0.2 m, höjd 1 m och vikten 15 kg. Lösningsförslag: Vi behöver veta Miss Inertias masströghetsmoment kring rotationsaxeln i de två tillstånden. Vi tar fram lärobok i dynamik, slår upp formelsamling för masströghetsmoment och beräknar bidraget från huvud, bål och ben. 52 Mekanik-Dynamik och Mathematica 2 1 0.12 4 hbb 5 3 0.22 15 15 2 0.22 HH/ITN/BN Huvud Bål Ben 10 0.496 Sedan masströghetsmoment för en utfälld arm, använd Steiner's sats. 1 1 0.052 5 ua 4 2 0.7 0.72 5 5 0.2 12 2 1.71979 samt masströghetsmoment för en infälld arm, använd Steiner's sats igen. 1 0.052 5 ia 5 0.05 2 0.2 2 0.31875 Slutligen vid frånvaro av yttre moment bevaras rörelsemängdsmomentet, det vill säga Iföre Ωföre vinkelhastigheten sedan hon fällt in armarna Solve hbb 2 Ωefter 17.3603 5 ua 2 hbb ia Iefter Ωefter . Detta ger oss direkt Ωefter Exempel: En trådrulle ligger på bordet enligt figur när sömmerskan drar i tråden vid A åt höger. Åt vilket håll rullar då trådrullen? Lösningsförslag: Frilägg trådrullen och benämn trådkraften åt höger för P. Vi får ekvationerna om x räknas positiv åt höger och positiv medurs. Vi får ekvationerna ekv mg 0, m x '' t P J mg '' t 0, m x t .. Ingen glidning x F, rP F P, J t FR Pr .. R . Lägg denna ekvation till de övriga. AppendTo ekv, x '' t 0, m x t mg O RF F R '' t P, J t FR P r, x t R t .. Svaret på vår brännande fråga ges av x. svar Solve ekv, x '' t , PR R x t r m R2 J P R , t J '' t , , F r m R2 Simplify P J , mg, F J mrR m R2 Här är både P, r, R, m och J positiva, så resan går åt höger eftersom R erforderlig friktionskoefficient för att upprätthålla rullning. Μ F Abs . svar P J Μ mrR mg m R2 J First r. Trådrullen rullar alltså upp sig på snörstumpen! Slutligen HH/ITN/BN Mekanik-Dynamik och Mathematica 53 Exempel: En tavla är uppriggad i ett bilsäte enligt figur. Plötsligt bromsar bilen in med en konstant g horisontell retardation a 3 . Friktionen kan försummas i punkten A men ej i punkten B där friktionskoefficienten är 0.5. Undersök om tavlan kommer att börja glida vid B. Avgör frågan genom att bestämma erforderlig Μe för att undvika glidning. Sök även maximal retardation a som bilen kan ha för att tavlan fortfarande ska stå kvar. Lösningsförslag: Antag att tavlan har höjden h och massan m. Efter friläggning och lite geometriska överläggningar (rita!) får vi rörelseekvationen i x-led, jämvikt i y-led samt momentjämvikt kring tyngdpunkten. Lös ut alla krafter. NF Solvem x '' t NA Cos 7 NA Sin 7 h Cos 30 2 NB Sin 7 FB Cos 7 NB Cos 7 FB Sin 7 h 7 NA Sin 30 7 2 N A , NB , FB First mg , 0, h NB Cos 30 7 FB 0, Tp 2 Simplify 1 N A m cos 7 g tan 7 tan 23 sin 7 x t , FB 1 tan 7 x t , tan 23 2 1 m g cos 7 NB m cos 7 g tan 7 tan 23 tan 7 tan 23 1 x t 2 Tavlan ska ha samma retardation som bilen. g NF . x '' t Simplify 3 1 N A g m cos 7 3 tan 23 tan 7 6 1 NB tan 23 3 1, 1 g m sin 7 3 cos 7 3 , FB g m cos 7 tan 7 tan 23 3 3 tan 23 1 6 Erforderlig Μe . FB g . NF . x '' t Μe Simplify 3 NB cos 7 tan 7 tan 23 2 sin 7 Μe 3 3 tan 23 1 3 cos 7 N 0.313373 Mindre än given möjlig 0.5, så svaret på första delfrågan är att tavlan inte glider vid B. Vad gäller maximal retardation har vi två fall att ta ställning till. Glidning vid B eller att tavlan lättar vid A, det vill säga "roterar" framåt i sätet. Låt oss undersöka det sista fallet som ges av villkoret N A 0. ekv NA 0 . NF 1 m cos 7 g tan 7 tan 23 1 tan 7 tan 23 x t 2 Bestäm motsvarande maximala retardation på bilen i detta fall. amax Solve ekv, x '' t g tan 7 First tan 23 x t tan 7 tan 23 1 Kontrollera om friktionen vid B klarar av detta detta. 0 54 Mekanik-Dynamik och Mathematica HH/ITN/BN FB ΜNA . NF . amax 0 Simplify NB tan 23 ΜNA N 0 0.424475 Mindre än given möjlig 0.5 så detta fall sätter gränsen för bilens retardation till Exempel: Två lådor ligger ovanpå varandra på ett glatt bord, se övre fig. En kraft anbringas på den övre lådan. Bestäm hur lång tid det tar innan vi har situationen i den undre figuren. Hur lång sträcka har den undre lådan då hunnit glida ? g för att tavlan inte ska "rotera" framåt i sätet. 3 0.5 m 100 N μ=0.8 5 kg 8 kg μ=0 100 N 5 kg 8 kg .. Lösningsförslag: Efter friläggning och införande av koordinater och krafter enligt figur nedan har vi standard Newton mx F för de två kropparna. Som sagt, separation är vårt ständigt återkommande problem. Nu är det dax igen. Här är det ganska lätt eftersom alla integrationskonstanter blir noll. Genomför gärna kalkylen för hand! Nu låter vi DSolve göra hela jobbet direkt. P x1 F=μm1g x2 x1å2 DSolve P t2 x1 t m2 m1 x1 '' t P Μ m1 g, m2 x2 '' t Μ m1 g, x1 0 0, x1 ' 0 0, x2 0 0, x2 ' 0 0 , x1 t , x2 t , t First g Μ m1 t 2 2 m1 g Μ m1 t 2 , x2 t m1 ODE för m1 ODE för m2 BV för m1 BV för m2 2 m2 Vi inleder med en liten reseberättelse. Plot Evaluate x1 t , x2 t , x1 t x2 t , L . x1å2 . m1 5, m2 8, P 100, L 0.5, Μ 0.8, g 9.81 PlotStyle Orange, Green, Blue, Red , AxesLabel "t s ", "x1 t ,x2 t ,x1 t x2 t ,L" , t, 0, 0.4 , x1 t ,x2 t ,x1 t x2 t ,L 1.0 0.8 0.6 0.4 0.2 0.1 0.2 0.3 0.4 t s När klockan är t har den lilla lådan hunnit till framkanten på den stora och således flyttat sig sträckan x1 t x2 t L, där L är given i uppgiften. I grafen ovan inträffar detta då den blå kurvan (x1 t x2 t ) skär den röda (L). Detta bestämmer restiden till det att vi har situationen i den undre figuren. ekv x1 t x2 t L . x1å2 HH/ITN/BN Mekanik-Dynamik och Mathematica P t2 g Μ m1 t 2 g Μ m1 t 2 2 m1 55 L 2 m2 Nu är det bara att lösa ut den. T Solve ekv, t 2 L 2 t L , t gΜ g Μ m1 P m2 m1 gΜ g Μ m1 P m2 m1 Eftersom vi inte befattar oss med negativa restider får vi resvägarna x1å2 . T 2 Simplify 2 L m2 g Μ L m1 x1 gΜ g Μ m21 g Μ m1 P m2 m1 LP g Μ m2 m1 2 m2 P g Μ L m21 L , x2 gΜ g Μ m1 P m2 m1 g Μ m21 g Μ m2 m1 m2 P Så med lite numerik. x1å2 . T 2 x1 0.371468 . m1 5, m2 8, P 0.838416, x2 0.371468 100, L 0.5, Μ 0.8, g 9.81 0.338416 Det kanske är av intresse att veta att om det finns en minsta kraft P för att den lilla lådan överhuvudtaget ska nå fram. Detta kritiska värde ges av att lådorna glider parallellt i all evighet Pkrit Solve x1 t g Μ m1 m1 P m2 x2 t . x1å2, P Simplify m2 Som med våra övriga värden blir Pkrit P . m1 5, m2 8, Μ 0.8, g 9.81 63.765 Exempel: Ett löphjul med mått enligt figur har massan 15 kg och tröghetsradien 160 mm. Det släpps med en medurs rotation på 20 rad s i ett spår som lutar 10 . Mellan axeltapparna och kanten gäller att friktionskoefficienten Μ 0.3. För övrigt är hjulet ej i kontakt med spåret. a Bestäm tidpunkt och läge då glidning upphör. b Bestäm tidpunkt och läge när det vänder. c Sök slutligen tidpunkt, hastighet och vinkelhastighet då det återkommer till startpunkten. Lösningsförslag: Inför x-axeln uppåt längs banan, y-axeln vinkelrät däremot samt krafterna. Teckna sedan rörelseekvationerna data m ekv m x '' t m y '' t J mx t 15, r '' t F 0.15, J N 15 0.1602 , Ω 20, Μ 0.3, Α medurs. Frilägg och för in de "vanliga" 10 , g 9.81; m g Sin Α F, m g Cos Α , rF g m sin Α , m y t .. Tp N g m cos Α , J t Fr Här är nu y 0 det vill säga y 0 under hela resan. Från början har vi dessutom glidning, det vill säga fullt utbildad friktion F ΜN. Strängt taget ska vi kontrollera efteråt att Μ verkligen inte klarar av att "rycka" igång hjulet i rullning, typ "kuggstång". Men vi litar på problemförfattaren och löser ut de obekanta. 56 Mekanik-Dynamik och Mathematica glider Solve Flatten ekv, y '' t 0, F ΜN x '' t , y '' t , '' t , N, F HH/ITN/BN , First g Μ m r cos Α x t g Μ cos Α g sin Α , y t 0, t ,N g m cos Α , F g Μ m cos Α J Nu är det bara att vaska fram (ODE), ode glider 1, 3 . Rule Equal g Μ m r cos Α x t g Μ cos Α g sin Α , t J applicera (BV) och lösa (BVP). x DSolve g ode, x 0 x t , 1 x t 0, x ' 0 t ,t 0, sin Α , t ' 0 Ω , 2J 2 När den inte glider längre rullar den. Då är x t Ta 0, g Μ m r t2 cos Α 2J tΩ g t2 Μ cos Α 0 First Solvex ' t r ' t . Dx r t , så svaret på delfråga a. g, t, t First J rΩ t g J Μ cos Α x . Ta g Μ m r2 cos Α J sin Α Simplify J 2 r2 Ω2 Μ cos Α J rΩ x m r2 g Μ cos Α J J sin Α g Μ cos Α J x m r2 m r2 2 g Μ cos Α J J rΩ sin Α J sin Α J r Ω2 Μ cos Α 2 J m r2 2 g Μ cos Α J 2 J sin Α mr 2 , 2 J sin Α J sin Α 2 . Ta . data g x 0.801684 0.383947, 0.801684 10.5765 För att få svaret på delfråga b måste vi börja om från början med rörelseekvationerna men med rullvillkor. rullar Solve Flatten ekv, y '' t 0, x '' t x '' t , y '' t , '' t , N, F g m r2 sin Α x t m r2 J r '' t First g m r sin Α ,y t 0, t J m r2 , g J m sin Α ,N g m cos Α , F m r2 J Nu är det bara att vaska fram (ODE) ode rullar 1, 3 . Rule g m r2 sin Α x t m r2 J Equal g m r sin Α , t J m r2 Meka ihop (BV) det vill säga de värden som gäller då glidning upphör och rullning inträder. Detta hämtas från a. Använd samma klocka, det vill säga den vi startade när spektaklet började. bv x x g, Dx g, t . Rule J rΩ g m r2 J Μ cos Α J sin Α x J 2 r2 Ω2 Μ cos Α sin Α J rΩ g m r2 J Μ cos Α J sin Α Equal . Ta 2 g m r2 J Μ cos Α J sin Α 2 J r Ω Μ cos Α sin Α m r2 J Μ cos Α J sin Α bv . data x 0.801684 x 0.801684 Lös (BVP). 0.383947 0.957852 0.801684 0.801684 10.5765 6.38568 Simplify J rΩ g m r2 J Μ cos Α J sin Α J rΩ g m r2 J Μ cos Α J sin Α J r Ω2 m r2 2 J Μ cos Α 2 J sin Α 2 g m r2 J Μ cos Α J sin Α 2 J Ω Μ cos Α sin Α m r2 J Μ cos Α J sin Α HH/ITN/BN x Mekanik-Dynamik och Mathematica DSolve r ode, bv , x t , r J g 2 m r t2 sin2 Α x t J g 2 m r t2 sin2 Α t x r mr First J r Ω2 m r2 g m r t sin Α g Μ t cos Α J 2 J Ω J sin Α , m r3 Ω2 m r2 2 J Ω g Μ t cos Α J g m r t sin Α J sin Α 1.59667 t 0.640013, 2.65617 t2 t Äntligen ges nu svaret på delfråga b av villkoret att x t Tb Simplify ExpandAll 0.398425 t2 x t 2 m r Μ cos Α J . data m r2 2 g J t Ω sin Α 2 2 g J ,t 2 g J t Ω sin Α m r2 Μ cos Α J 2 g J t 57 Solve x ' t 0 .D x r, t ,t 10.6445 t 3.75008 0 i vändläget. First JΩ t g m r sin Α x r . Tb Simplify J 2 Ω2 JΩ x g m r sin Α x r Μ 1 tan Α , m r2 2 g m Μ cos Α J J Ω2 JΩ J sin Α g m r sin Α JΜ tan Α J 2 g m r Μ cos Α J m r2 m r2 J sin Α . Tb . data x 2.00373 0.959639, 2.00373 14.4144 Slutligen delfråga c. Hur mycket är klockan när vi kommer hem och hur snurriga är vi? Tc Solve x t 1 0 .x r, g 2 J 2 Ω2 J t t Simplify m r2 Μ cos Α sin Α Μ cos Α J m r2 J sin Α g J2 Ω m r2 g J ΜΩ J sin Α tan Α Μ J g 2 m r sin Α m r2 J , tan Α 1 g 2 J 2 Ω2 J t m r2 Μ cos Α m r2 sin Α Μ cos Α J J sin Α g J2 Ω m r2 gJ ΜΩ J sin Α tan Α m r2 Μ J g 2 m r sin Α J tan Α Tc t . data 0.451771 , t 3.55569 Här duger bara den sista eftersom vår rullande modell endast gäller vid just rullning, det vill säga då klockan har passerat det klockslag som bestämdes i delfråga a. Stämmer dessutom bra eftersom vi kommer hem efter det att vi vänt Tc Tc 2 1 g 2 J 2 Ω2 J t m r2 Μ cos Α sin Α Μ cos Α J m r2 J sin Α g J2 Ω g J ΜΩ J m r2 sin Α tan Α g 2 m r sin Α m r2 Μ J J tan Α Så snurrigheten vid hemkomsten är D x r, t 1 . Tc Simplify g 2 J 2 Ω2 J x m r2 Μ cos Α sin Α Μ cos Α J m r2 J sin Α g J2 Ω m r2 g J ΜΩ J sin Α tan Α g 2 m r sin Α Μ J m r2 J tan Α r g 2 J 2 Ω2 J m r2 Μ cos Α sin Α Μ cos Α J m r2 J sin Α g sin Α J m r2 Μ J m r2 J tan Α , 58 Mekanik-Dynamik och Mathematica 1 g 2 J 2 Ω2 J m r2 Μ cos Α m r2 sin Α Μ cos Α J J sin Α HH/ITN/BN m r2 gJ ΜΩ J g J2 Ω sin Α tan Α g 2 m r sin Α m r2 Μ J J tan Α g 2 J 2 Ω2 J m r2 Μ cos Α m r2 sin Α Μ cos Α J J sin Α g sin Α J m r2 Μ J m r2 J tan Α D x r, t . Tc . data x 3.55569 x r 1.23668, . Tc 1 3.55569 Simplify g 2 J 2 Ω2 J x 8.24453 m r2 Μ cos Α m r2 sin Α Μ cos Α J J sin Α g J2 Ω gJ ΜΩ J m r2 sin Α tan Α g 2 m r sin Α Μ J m r2 0, J tan Α 1 g 2 J 2 Ω2 J m r2 Μ cos Α m r2 sin Α Μ cos Α J J sin Α g J2 Ω g J ΜΩ J sin Α tan Α g 2 m r sin Α Μ J m r2 J r Ω2 J tan Α x r m r2 2 g sin Α Μ J m r2 J tan Α . Tc . data x 3.55569 8.99505 10 16 , 3.55569 8.01684 Om man inte är intresserad(?) av de ibland omfattande symboliska uttrycken kan man naturligtvis ta sig igenom uppgiften rent numeriskt genom att redan från början applicera data, men då tappar man många lärospån på vägen Avslutnings några förtydligande(?) bilder över situationen. Eftersom det intressantaste händer under en väldigt kort tid i början använder vi logaritmisk tidsskala. Vi börjar med x t blå och x t röd. Vidare är tre vertikala linjer inlagda för att markera de tre tidszonerna: glidning, rullning till vändläge samt rullning hem. Τa , Τ b , Τc t . Ta , T b , Tc LogLinearPlotEvaluatex t x t PlotStyle GridLines .x . data; . x g, x ' t r, x ' t .D x . Dx r, t g, t 1 UnitStep t UnitStep t Τa Blue, Red , GridLinesStyle GrayLevel 0.8 , Τa , Τb , Τc , None , AxesLabel "t s ", "x t Τa . data, t, 10 3 , Τc , m , x t m s " xt m,xt ms 1.0 0.5 0.0 0.5 t s 10 3 10 Samma melodi för vinkeln, 2 0.1 t blå och LogLinearPlotEvaluate t PlotStyle GridLines .x 1 t röd. t r, .x ' t g, ' t .D x r, . Dx t g, t 1 UnitStep t UnitStep t Τa Blue, Red , GridLinesStyle GrayLevel 0.8 , Τa , Τb , Τc , None , AxesLabel "t s ", " t Τa . data, t, 10 3 , Τc , rad , t rad s " HH/ITN/BN Mekanik-Dynamik och Mathematica t rad , 59 t rad s 20 15 10 5 0 5 10 3 t s 2 10 0.1 1 Exempel: En låda skjuts iväg med farten v1 6 m s från en punkt A uppför ett lutande plan med Θ 20 . Mellan lådan och planet är friktionskoefficienten Μ 0.2. a Bestäm tidpunkt och läge då lådan vänder. b Bestäm tidpunkt och fart v2 då den återkommer till A. Lösningsförslag: Inför x-axeln uppåt längs banan och y-axeln vinkelrät däremot, snett uppåt. Frilägg och för in de "vanliga" krafterna. Teckna sedan rörelseekvationerna för resan uppåt. data 0.2, Θ Μ ekv m x '' t m y '' t mx t F 20 , v1 N g m sin Θ , m y t .. x t N g m cos Θ 0 under hela resan. Vi har dessutom glidning under hela resan, det vill säga fullt utbildad friktion Solve Flatten g Μ cos Θ 9.81 ; m g Sin Θ F, m g Cos Θ Här är nu y 0 det vill säga y F ΜN. Lös ut de obekanta. glider 6, g ekv, y '' t g sin Θ , y t 0, N 0, F g m cos Θ , F ΜN , x '' t , y '' t , N, F First g Μ m cos Θ Nu är det bara att vaska fram (ODE), applicera (BV) och lösa (BVP). upp DSolve 1 x t glider 1 g Μ t2 cos Θ . Rule g t2 sin Θ 2 Equal, x 0 Solve x ' t 0 . D upp, t , t v1 t g Μ cos Θ v1 , x t , t First 2 t v1 När den vänder är hastigheten noll, det vill säga x t Ta 0, x ' 0 0, så svaret på delfråga a. First sin Θ upp . Ta . data x 1.15409 3.46227 Innan vi går vidare måste vi först kontrollera att vi får åka hem igen eller om vi är dömda att stanna i vändläget för gott, det vill säga om lådan verkligen börjar glida nedåt igen. Detta bestäms av om givet Μ klarar av att hålla den kvar i vändläget eller ej. Bestäm därför erforderligt Μr för att stanna kvar och jämför med Μ. FN Solve ekv . x '' t F g m sin Θ , N F Abs . FN N Μr tan Θ Μr g m cos Θ . data N 0, y '' t 0 , F, N First 60 Mekanik-Dynamik och Mathematica HH/ITN/BN 0.36397 Μ . data Μr True Visst, delfråga b är sund. Tyvärr måste vi börja om från början med rörelseekvationerna, ty under hemresan har ju friktionskraften bytt tecken. glider x t Solve Flatten g Μ cos Θ ekv, y '' t g sin Θ , y t 0, N 0, F ΜN g m cos Θ , F , x '' t , y '' t , N, F First g Μ m cos Θ Meka ihop (BV), det vill säga de värden som gäller i vändpunkten. Detta hämtas från a. Använd samma klocka, det vill säga den vi startade när spektaklet började. bv x upp, D upp, t v1 g Μ cos Θ g sin Θ x . Rule Equal . Ta Simplify v21 2 g Μ cos Θ 2 g sin Θ v1 g Μ cos Θ g sin Θ 0 bv . data x 1.15409 3.46227 0 x 1.15409 Nu är det bara att vaska fram (ODE) och lösa (BVP). ned DSolve glider 1 . Rule Equal, bv , x t , t First 1 x t 2 g Μ cos Θ sin Θ g 2 Μ3 t2 cos3 Θ 2 g 2 Μ2 t2 sin Θ cos2 Θ g 2 Μ t2 sin2 Θ cos Θ Äntligen ges nu svaret på delfråga b av villkoret att x t Tb t Solve x t 0 . ned, t g 2 Μ4 v21 cos4 Θ g 2 Μ3 cos3 Θ ned, D ned, t . Tb 0. , x 3.29445 x 3.29445 2 g Μ2 t v1 cos2 Θ 2 g t v1 sin2 Θ 2 Μ v21 cos Θ 0 vid hemkomsten. First 2 g 2 Μ3 v21 sin Θ cos3 Θ g 2 Μ2 sin Θ cos2 Θ g 2 t2 sin3 Θ 2 g 2 Μ v21 sin3 Θ cos Θ g 2 Μ sin2 Θ cos Θ g 2 v21 sin4 Θ g Μ2 v1 cos2 Θ g v1 sin2 Θ g 2 sin3 Θ . data 3.23523 Avslutnings en förtydligande(?) bild över situationen. Läget x t blå och hastigheten x t röd. Vidare är två vertikala linjer inlagda för att markera de två tidszonerna utresa och hemresa. Τa , Τ b t . Ta , T b . data; PlotEvaluate x t . upp, x ' t PlotStyle AxesLabel . D upp, t 6 4 2 2 UnitStep t Τa . data , t, 0, Τb , x t . ned, x ' t . D ned, t UnitStep t Τa Blue, Red , GridLinesStyle GrayLevel 0.8 , GridLines Τa , Τb , None , "t s ", "x t m , x t m s " xt m,xt ms 0.5 1 1.0 1.5 2.0 2.5 3.0 t s HH/ITN/BN Mekanik-Dynamik och Mathematica 61 Exempel: Två bilar har samma hastighet v0 på en rak väg. Avståndet mellan dem är d. Plötsligt börjar den främre bilen bromsa med konstant retardation a f . Bestäm minsta avstånd d för att undvika kollision om föraren i den bakomvarande bilen har reaktionstiden Τ och sedan bromsar med konstant retardation ab . Undersök speciellt fallet 90 km h , a f 4k, ab 5k, k 1 ms2 , Τ 1 s. v0 25 m s Lösningsförslag: Starta klockan då den främre bilen börjar bromsa och rulla ut måttbandet från den position som den bakre bilen befinner sig. Vi får då följande (BVP) data v0 25, af 4, ab 5, Τ 1 ; xfÅbBVP DSolve xf '' t af , xf ' 0 v0 , xf 0 d, xb '' t ab UnitStep t Τ , xb ' 0 v0 , x b 0 xf t , x b t , t Simplify First t2 a f x f t ConditionalExpression t v0 , Τ d 2 0, xb t 0 , 1 ConditionalExpression ab t 2 Τ2Θ t Τ 2 t v0 , Τ 0 För att kollision inte ska ske innan den bakomvarande börjat bromsa krävs att d är minst SimplifySolvexb t 1 . xfÅbBVP . t Τ, d, Τ 0 Τ2 a f d 2 Inför tidsintervallet t xfÅb xf t Τ passar vi på att förenkla ekvationerna lite SimplifyxfÅbBVP, Τ t2 a f x f t 0, t 1 d 2 t v0 , xb t 2 ab t Τ 2 Τ 2 t v0 Om detta är uppfyllt skall kollision undvikas då båda bilarna bromsar och klockan t Τ. Minsta avstånd d ges nu av att gälla för något t och d, det vill säga "tangentkrav". dt Solvexb t ab . xfÅb, xb ' t xf ' t Τ2 a b a f Τ ab t xf t ,d af 2 ab af Speciellt har vi för den önskade situationen dt . data t 5, d 10 Bilarnas restider till stillastående. Tf Solvexf ' t v0 t 0 . DxfÅb, t, t First 0 . DxfÅb, t, t First af Solvexb ' t Tb t Τ ab v0 ab Bromssträckor. xf t d . xfÅb . Tf v20 2 af xb t . xfÅb . Tb Simplify . DxfÅb, t, t, d First xf xb xf xb ska 62 Mekanik-Dynamik och Mathematica HH/ITN/BN v20 Τ v0 2 ab Nu kan vi äntligen inspektera bilarnas resa genom vädret i skarpt läge. ShowPlotxf t . xfÅbBVP . dt 2 PlotRange Plotxb t . data, t, 0, t . Tf . data , PlotStyle 0, 100 , AxesLabel . xfÅbBVP . dt 2 80 60 40 20 1 2 3 4 5 s ", "xf t ,xb t 6 t s Red, m " , . data, t, 0, t . Tb . data , PlotStyle x f t ,xb t m 100 0 "t Blue
© Copyright 2024