Simulasi Gempa Pada Gedung Tinggi Dengan Tuned Mass Damper (TMD) Topik Khusus Sains Komputasi (Dosen: Prof. Zaki Suud) Dian Nuraiman (20912017) Muhaza Liebenlito(20912020) April 15, 2013
1
Penda endah hulua uluan n
Gedung pencakar langit dapat berfungsi sebagai tempat untuk melakukan aktifitas bisnis ataupun ataupun sebagai sebagai simbol kedigday kedigdayaan aan suatu negara. Namun Namun untuk untuk merancang merancang gedung tersetersebut tidak mudah. Terdapat erdapat faktor-fa faktor-faktor ktor yang perlu dipertimba dipertimbangk ngkan an untuk untuk keselama keselamatan tan penggunanya seperti konstruksi bangunan, ketahanan terhadap bencana alam, dan lain-lain. Bencana alam merupakan hal yang terjadi suatu waktu dan sulit untuk diprediksi, termasuk gempa bumi. Sehingga Sehingga gedung gedung pencakar pencakar langit perlu untuk meminimalk meminimalkan resiko resiko bencana yang terjadi. Tuned Mass Damper (TMD) merupakan salah satu teknik yang mampu meminimal inimalk kan gun guncan cangan gan yang yang terjad terjadii akibat akibat gempa gempa bumi. bumi. TMD diinspir diinspirasi asi oleh oleh dua benda benda dengan massa yang berbeda yang dihubungkan dengan menggunakan pegas sehingga gerak salah satu benda akan teredam akibat dari perbedaan massa yang terjadi pada konfigurasi tersebut. Adapun Tujuan dari simulasi ini adalah untuk melihat pengaruh TMD pada gedung tinggi pada saat terjadi gempa. Asumsi yang kami gunakan dalam mendesain simulasi ini: •
Gedung yang akan disimulasikan yaitu Taipei 101.
•
Desain simulasi dua dimensi.
•
•
Lantai dan dinding diasumsikan sebagai pegas dengan dua bola bermassa yang saling berinteraksi. Pegas diposisikan menyilang antar satu lantai dengan lantai yang lain.
1
2 2.1 2.1
Teori eori Pendu enduku kung ng Shea Sh earr Force orce
Pada simulasi ini, shear force diasumsikan sebagai gaya yang dapat bekerja pada pegas yang diposisikan saling menyilang di tiap lantai. Pada Figure 1 mengilustrasikan bagaimana shear force force tersebut tersebut bekerja, bekerja, dimana dimana bidang bidang tampang tampang lintang lintang ini selalu horizontal horizontal pada saat shear force bekerja, sehingga material atau gedung tersebut dapat mengalami perenggangan secara linear. linear. Perenggan Perenggangan gan secara secara linear linear ini disimulasik disimulasikan an menggunak menggunakan an gaya gaya pegas menyilang menyilang.. Ketinggian tampang lintang awal dapat berubah menjadi lebih rendah ketika shear force bekerja. bekerja. Ini dikarenak dikarenakan an asumsi bahwa bahwa dinding dinding dari bangunan adalah b enda tegar dengan kekakuan sempurna. Perhitungan shear force dilakukan dengan cara menghitung percepatan bola bermassa di masing-masing lantai.
Figure 1: Pergerakan bola bermassa
2.2
Peak Groun Ground d Accel Accelera eratio tion n (PGA (PGA))
PGA merupakan merupakan perce p ercepatan patan terbesar dari pergerakan pergerakan tanah pada saat gempa bumi. Jika Jika pola gerakan yang diberikan adalah sinusoidal, maka percepatan terbesar dapat dirumuskan pada persamaan: persamaan: = Aω 2 amax = Aω dimana A adalah besar amplitudo getaran yang terjadi pada saat gempa dan ω adalah frekuensi getarannya.
2.3 2.3
Tun uned ed Mass Mass Dam Damper per (TMD (TMD))
Persamaan gerak untuk massa primer atau massa gedung seperti yang diperlihatkan pada Figure 2 adalah: p (1 + m ¯ )¨u + 2ωm 2ωmu˙ + ω + ω 2 u = ¯u ¨d m m −
2
¯ didefinisikan sebagai perbandingan massa antara massa gedung dengan massa damper, m ¯ = mm , u˙ adalah kecepatan, u ¨ percepatan, adalah faktor damping dari massa primer. m d
Figure 2: Sistem SDOF-TMD.
Sehingga persamaan gerak untuk TMD adalah: ¨d + 2 2 d ωd ˙ud + ω + ωd2 ud = u
u¨d
−
Tujuan menambahkan massa damper adalah untuk mengontrol vibrasi dari struktur ketika terjadi terjadi guncangan. Massa damper mempunyai parameter yaitu massa m d , konstanta pegas kd dan koefisien damping cd . Damper dibangun dibangun untuk menjadi menjadi frekuensi frekuensi poko p okok k dari struktur sehingga ωd = ω ¯ kd = mk Massa primer bergerak bergerak secara secara periodik sinusoidal sinusoidal sinΩt p = p = pˆ sinΩt maka respon yang diberikan adalah sin(Ωt + δ + δ 1 ) u = uˆ sin(Ωt ˆd sin(Ωt sin(Ωt + δ + δ 1 + δ + δ 2 ) ud = u dimana uˆ dan δ dan δ dinotasikan dinotasikan sebagai pergeseran (displacement ) benda dan phase shift , secara secara berurut. Untuk mengatasi masalah kondisi resonan, maka solusinya adalah: pˆ ˆ = ¯ 1/(1 + (2/ (2/m ¯ + 1/ 1 /2d )2) u mk ˆd = (1/ (1/2d)ˆu u tan δ 1 =
(2/ (2/m ¯ + 1/ 1/2d)
−
tan δ 2 = 3
π/2 π/2
−
(1)
Persamaan di atas menunjukan bahwa respon dari TMD adalah 90 dari fase dengan respon yang terjadi pada massa primer. Dari persamaan (1) didapat ◦
d =
m 1/1 + (2/m (2/m + + 1/ 1 /2d )2 2
(2)
Persamaan (2) menunjukan kontribusi relatif dari TMD dari total damping.
3
Hasi Ha sill da dan n Pemba embaha hasa san n
Desain gedung pada simulasi ini diasumsikan sebagai kumpulan bola bermassa yang dihubun hu bungk gkan an oleh oleh pegas pegas satu satu dengan dengan yang yang lainn lainnya ya,, dimana dimana setiap setiap lantai lantai terdir terdirii dari dari tiga tiga bola bermas b ermassa. sa. TMD diletakk diletakkan di lantai lantai paling atas yang yang dihubungk dihubungkan an dengan dua bola bermassa bermassa dari kiri dan kanan. Seperti Seperti yang terlihat terlihat pada Figure 3. Pada simulasi simulasi ini metode untuk menghitung perubahan posisi bola pada waktu t dapat menggunakan algoritma Verlet. Karena Karena berdasark berdasarkan an teori penduk p endukung ung 2.3 perubahan perubahan posisi p osisi yang terjadi terjadi pada waktu waktu t merupakan . Dengan demikian ketika didapat perubahan posisi bola bermassa pada waktu t, selanjutnya menghitung dengan shear force yang terjadi dari interaksi perubahan bola bermassa tersebut.
Figure 3: Asumsi Desain Gedung. Tabel di bawah ini adalah data-data pendukung untuk melakukan simulasi.
4
Tinggi Gedung Jumlah Lantai Massa Gedung Massa Damper Juml Jumlah ah Bola Bola per per Lan Lantai tai PGA Amplitudo Massa Konsta Konstant ntaa Pega Pegass Damper Damper Konsta Konstant ntaa Pega Pegass Dindi Dinding ng
50 9 m 19 6 700000 ton 73 0 t o n 3 0.5g, 1.7g, 2.7g (g = 9.8m/s2 ) 0.3 119047.19 kg 6.25e09 6.25e09 6.25e09 6.25e09
Dari Figure 4 terlihat bahwa gedung yang menggunakan TMD menghasilkan shear force yang lebih kecil dibandingkan dengan gedung yang tidak menggunakan TMD. Lantai atas pada gedung yang menggunakan TMD mempunyai shear force yang hampir mendekati nilai 0 diband dibanding ingk kan dengan lantai lantai 1 s/d 18. Hal ini dikare dikarenak nakan an gedung gedung pada lantai lantai 1 s/d 18 lebih dekat dengan gempa dibandingkan dengan bola TMDnya.
5
(a) PGA 0.5 (Non-TMD)
(b) PGA 0.5 (TMD)
(c) PGA 1.7 (Non-TMD)
(d) PGA 1.7 (TMD)
(e) PGA 2.7 (Non-TMD)
(f ) PGA 2.7 (Non-TMD)
Figure 4: Grafik Perbandingan Shear Force.
4
Kesimpulan ulan
Hasil Hasil simula simulasi si di atas atas menunj menunjukk ukkan an simula simulasi si dengan dengan tuned tuned mass damper damper (TMD) (TMD) mengmenghasilkan shear force yang lebih kecil dengan tanpa menggunakan TMD. Hal ini menunjukkan bahwa TMD mampu meredam getaran gempa. Simulasi dengan TMD menunjukkan bahwa gedung hanya menerima getaran di bawah lantai 20.
6
References [1] Misrha, Misrha, Rashmi. Rashmi. 2011. Application of Tuned Mass Damper For Vibration Control of Frame Structures Under Seismic Excitaions . Thesis. [2] Esteki, Esteki, Kambiz. Kambiz. 2011. Semi-Active Tuned Mass Damper for Seismic Applications . International Workshop Smart Materials, Structures and NDT in Aerospace Conference. [3] Zahrai Zahrai,, S.M. S.M. 2013. 2013. Seismic Design of Fuzzy Controller For Semi-Active Tuned Mass Damper Using Top Stories As The Mass . Asian Journal of Civil Engineering (BHRC) Vol. 14, No. 3. Semi-Active tive Tuned Tuned Mass Damper Damper Building Building System: System: Appli[4] Min-Ho Min-Ho Chey. Chey. 2010. Semi-Ac Earthquakee Engineering Engineering and Structural Structural Dynamics (EESD), 39. cation . Earthquak
7
Lampiran
Listing 1: gnu gnuplot.h plot.h // 15.2.2 15.2.2012 012 9:27 9:27 #ifndef GNUPLOT_H GNUPLOT_H 2 #ifndef #define GNUPLOT_H GNUPLOT_H 3 #define / 4 ** 1
5 6 7
[== "Gnuplot.h" "Gnuplot.h" This This class class serv serves es as easi easier er access access to gnup gnuplo lot t and its execu executa tabl ble e .txt .txt file
8 9
10
11 12 13 14 15
. The The Gnupl Gnuplot ot clas class s cons constr truc ucted ted with two .txt .txt file file crea create ted, d, one one for for storing storing data, one one for execut executin ing g comman commands ds from from the shell. shell. The defau default lt name name of those those files files are: are: plotcommand plotdata User User may may decl declar are e a Gnup Gnuplo lot t class class via via Gnuplot sine; (thi (this s will will open open a stre stream am to data data file file and and comm comman and d file file) )
16 17 18 19
. User User may inser insert t a text text to both both of the the comm comman and d and and data data file file via via objectname.co objectname.command< mmand<<" <" " objectname.da objectname.data<<" ta<<" "
20 21 22 23
. User User may may exec execut ute e the the comm comman and d in the the comm comman and d file file via via objectname.plot(); (thi (this s will close close the entir entire e data data strea stream m before before execu executi ting ng) )
24 25 26 27 28 29
. User User can do a rewr rewrit ite e to the data data file file or the comma command nd file file via via objectname.redata(); objectname.recommand(); (thi (this s will reop reopen en the the data data stream stream if they they are are closed closed) ) ==]
30
*/ class Gnuplot Gnuplot 32 class 33 { public: 34 char* plotcommand; 35 char* plotdata; 36 31
37 38 39
ofstream ofstream command; command; ofstream ofstream data;
40 41 42
bool commandFlag; commandFlag; bool dataFlag; dataFlag;
43 44
Gnuplot(char Gnuplot(char _plotcommand[ _plotcommand[]="plo ]="plotcomma tcommand.txt" nd.txt", , char _plotdata[]=" _plotdata[]=" plotdata.txt")
8
{
45
plotcommand=_plotcommand; plotdata=_plotdata; command.open(plotcommand); data.open(plotdata); command <<"unset <<"unset key\n" <<"set grid\n"; grid\n"; commandFlag=true; dataFlag=true;
46 47
48 49 50
51 52
}
53 54
virtual ˜Gnuplot() ˜Gnuplot() { }
55 56 57 58
void redata() redata() { if(dataFlag= if(dataFlag==false =false) ) { data.open(plotdata); dataFlag=true; } }
59 60 61 62 63 64 65 66
//C //Clo lose se the the data data plot plot file file (if (if it is nece necess ssar ary y in the the midd middle le of the the plot plot scheme) void closeData() closeData() 68 { 69 data.close(); 70 dataFlag=false; 71 } 72 67
73
void recommand() recommand() { if(commandFl if(commandFlag==fa ag==false) lse) { command.open(plotcommand); command <<"unset <<"unset key\n" <<"set grid\n"; grid\n"; commandFlag=true; }
74 75 76 77 78 79 80 81 82 83
}
84 85 86 87
88 89 90
91 92 93
void void plot() plot()//E //Exec xecute ute the gnuplo gnuplot t { command.close(); data.close(); string shellcommand shellcommand= = "D:\\Documen "D:\\Documents\\gn ts\\gnuplot\\ uplot\\binary binary\\gnup \\gnuplot lot e \"lo \"load ad ’"; ’"; shellcommand.append(plotcommand); shellcommand shellcommand.appen .append("’; d("’; unset out\" "); system(shellcommand.c_str());
94 95 96
commandFlag=false; dataFlag=false;
9
}
97 98
}; #endi ndif f // GNUPLOT GNUPLOT_H _H 100 #e 99
Listing 2: spher sphere.h e.h /** This is clas class s is an Sphe Sphere re clas class s using using arra array y algo algori rith thm m vers versio ion n 2 Th 15.2.2012 2 9:27 3 15.2.201 4 */ 1
5
#ifndef #ifndef SPHERE_H SPHERE_H #define SPHERE_H SPHERE_H 7 #define 6
8 9
#def #defin ine e maxNu maxNum m 1000 10000 0 //De //Defi fine ne max numb number er of Sphe Sphere res s in the the simu simula lati tion on
10 11
class class Sphere Sphere 13 { public: 14 double mass[maxNum] mass[maxNum], , radius[maxNum radius[maxNum]; ]; 15 double xpos[maxNum] xpos[maxNum], , ypos[maxNum], ypos[maxNum], zpos[maxNum]; zpos[maxNum]; 16 double xinit[maxNum xinit[maxNum], ], yinit[maxNum] yinit[maxNum], , zinit[maxNum] zinit[maxNum]; ; 17 double xvel[maxNum] xvel[maxNum], , yvel[maxNum], yvel[maxNum], zvel[maxNum]; zvel[maxNum]; 18 double xacc[maxNum] xacc[maxNum], , yacc[maxNum], yacc[maxNum], zacc[maxNum]; zacc[maxNum]; 19 double xaccTemp[max xaccTemp[maxNum], Num], yaccTemp[max yaccTemp[maxNum], Num], zaccTemp[maxN zaccTemp[maxNum]; um]; //for 20 velocity velocity verlet double xtotalF[maxN xtotalF[maxNum], um], ytotalF[maxN ytotalF[maxNum], um], ztotalF[maxNu ztotalF[maxNum]; m]; 21 12
22
int numOfSpheres numOfSpheres; ;
23 24
Sphere(int Sphere(int _numOfSphere _numOfSpheres);//co s);//const nst
25 26
virtua virtual l ˜Spher ˜Sphere() e() {}
27 28
/*the the force force appl applie ied d ussu ussual aly y not not more more that that thre three e kind kinds s of forc force e User User must must decla declare re the the force forces s outsid outside e the clas class s*/ void force0(); force0(); void force1(); force1();
29 30 31 32 33
//Th //This is funct functio ion n apply apply all of the force force funct functio ions ns void void tota totalF lF() () { force0(); force1(); }
34 35 36 37 38 39 40
//Mo //Move ve the Sphe Sphere re in one time time step step (the (the time time Step Step MUST MUST be an input input) ) void velVerletMove velVerletMove(doubl (double e timeStep); timeStep);
41 42 43 44
};
45
10
46 47
Sphere::Sphere(int Sphere::Sphere(int _numOfSphere _numOfSpheres) s) 49 { numOfSpheres numOfSpheres = _numOfSpheres _numOfSpheres; ; 50 48
51
for( for(in int t i = 0; i
52 53 54
55
mass[i]=1.0; xpos[i]=0.0; xpos[i]=0.0; ypos[i]=0.0; ypos[i]=0.0; zpos[i]=0.0; zpos[i]=0.0; xinit[i]=0.0 xinit[i]=0.0; ; yinit[i]=0.0; yinit[i]=0.0; zinit[i]=0.0; zinit[i]=0.0; xvel[i]=0.0; xvel[i]=0.0; yvel[i]=0.0; yvel[i]=0.0; zvel[i]=0.0; zvel[i]=0.0; xacc[i]=0.0; xacc[i]=0.0; yacc[i]=0.0; yacc[i]=0.0; zacc[i]=0.0; zacc[i]=0.0; xtotalF[i]=0 xtotalF[i]=0.0, .0, ytotalF[i]=0 ytotalF[i]=0.0; .0; ztotalF[i]=0. ztotalF[i]=0.0; 0; xaccTemp[i]= xaccTemp[i]=0.0; 0.0; yaccTemp[i]= yaccTemp[i]=0.0; 0.0; zaccTemp[i]=0 zaccTemp[i]=0.0;//f .0;//for or velocity velocity verlet
56 57 58 59 60 61
62
}
63 64
}
65 66
void Sphere::velVe Sphere::velVerletMo rletMove(dou ve(double ble timeStep) timeStep) 68 { //Null //Nullifi ifies es all forces forces 69 for(in for(int t i=0; i=0; i
76 77 78 79
//=====force //=====force sum===== sum===== totalF(); //===================
80 81 82 83 84 85 86
//Upda //Update te the veloci velocity ty for(in for(int t i=0; i=0; i
87
//Upda //Update te veloci velocity ty here here xvel[i xvel[i] ] += timeSte timeStep p*0.5*(xacc[i]+xaccTemp[i]); yvel[i yvel[i] ] += timeSte timeStep p*0.5*(yacc[i]+yaccTemp[i]); zvel[i zvel[i] ] += timeSte timeStep p*0.5*(zacc[i]+zaccTemp[i]);
88 89 90 91 92 93 94
}
95 96 97 98
//th //the e move move one one step step loop loop for(in for(int t i=0; i=0; i
11
xaccTe xaccTemp[ mp[i] i] = xacc[i] xacc[i]; ; yaccTe yaccTemp[ mp[i] i] = yacc[i] yacc[i]; ; zaccTe zaccTemp[ mp[i] i] = zacc[i] zacc[i]; ;
99 100 101 102
xpos[i xpos[i] ] += timeStep*xvel xvel[i [i] ] + 0.5 0.5*timeStep*timeStep*xacc[i];
ypos[i ypos[i] ] += timeStep*yvel yvel[i [i] ] + 0.5 0.5*timeStep*timeStep*yacc[i];
zpos[i zpos[i] ] += timeStep*zvel zvel[i [i] ] + 0.5 0.5*timeStep*timeStep*zacc[i];
103 104 105 106 107 108 109 110
}
111 112 113 114
}
115 116 117 118
#endi #endif f // SPHERE_ SPHERE_H H
Listing 3: main.c main.cpp pp #include
#include 3 #include 4 #include 5 #include 6 #include 7 #include 1 2
8
using using namespa namespace ce std; std; #include "Gnuplot.h" "Gnuplot.h" 10 #include #include "Sphere.h" "Sphere.h" 11 #include #define e pi 3.14159 3.14159265 265 12 #defin 9
13
//=========== //================= ============ ============ ============= ============= ========== ==== SIMULATION SIMULATION PARAMETER PARAMETER const st int int numOf numOfFl Floo oors rs = 196; 196; // juml jumlah ah lant lantai ai 15 con const st int int numOf numOfRo Room oms s = 2; //FI //FIXED XED. . DO NOT NOT CHAN CHANGE GE THIS THIS 16 con 14
17 18
cons const t int nSphere nSpheres s = numOfF numOfFloo loors rs*3;
19 20
co cons nst t int int nSpri nSpring ngs s = 1000 10000; 0; //the //there re are are 1078 1078 spri spring ngs s
21
doub double le do doub uble le double e 24 doubl double le 25 doub 22 23
timeSte timeStep p = 0.002; 0.002; time timeEnd Ends s = 5; pngPlot pngPlotPer PerSte Step p = 25; grav gravAcc Acc = 9.8; 9.8;
26
double 28 double doubl uble e 29 do double 30 double 27
springConst[n springConst[nSpring Springs]; s]; springLength[ springLength[nSprin nSprings]; gs]; wallSpr wallSpring ingCon Const st = 6.25 6.25*1e9; shearSp shearSprin ringCo gConst nst = 6.25 6.25*1e9;
12
31 32
doubl double e damping dampingCon Const st = 0.5; 0.5; double double massPer massPerSph Sphere ere = 119047 1190471.19 1.19; ;
33 34 35
doubl double e massDam massDamper per = 700000 700000; ; double double springD springDamp amperC erCons onst t = 6.25 6.25*1e9;
36
bool bool create createMovi Movie e = true; true; stri string ng movieSp movieSpeed eed = "80"; "80"; int sphere sphereInde Index[5 x[5000 000][2 ][2]; ]; //corre //correspo sponds nds to the pair pair that that uses uses the spring spring 39 int attachment int pair pairIn Inde dex x = 0; // Init Initia iate te the init initia ial l coun counte ter r of how how many many spri spring ng pairs pairs 40 int 37 38
41
bool bool torque torquePlot Plot = true; true; //oscillation //oscillation parameter parameter bool l shak shake e = true true; ; 44 boo doub uble le ampl amplitu itude de = 0.3; 0.3; 45 do doubl ble e PGA PGA = 0.5 0.5*gravAcc; 46 dou double omega omega = sqrt(P sqrt(PGA/ GA/amp amplit litude) ude); ; 47 double 42 43
48
//============================================== DECLARATION ball(nSpheres); ); 50 Sphere ball(nSpheres pball = &bal &ball; l; 51 Sphere* pball 49
52
// //I I decl declar are e the gnup gnuplo lot t to plot plot the the resu result ltin ing g posi positio tion n Gnuplot spherePlot; spherePlot; 54 Gnuplot pspherePlot t = &spherePlot; &spherePlot; 55 Gnuplot* pspherePlo 53
56 57 58
Gnuplot shearPlot("sh shearPlot("shearCom earCommand.tx mand.txt", t", "shearData.t "shearData.txt"); xt"); Gnuplot* pshearPl pshearPlot ot = &shear &shearPlo Plot; t;
59 60
//======== //=============== ============= ============ ============= ============= ============ ====== SPRING FORCE & GRAVITY GRAVITY
61
void void Sphere Sphere::fo ::force rce0() 0() //Simp //Simple le constan constant t near near surfac surface e gravit gravity y 63 { for(in for(int t i i=0; =0; i
69
vo void id atta attach chSpr Sprin ing( g(in int t i, int int j, doub double le spri spring ngL, L, doubl double e spri spring ngC) C) //to //to attac attach h spri spring ng on sphe sphere re i and and sphe sphere re j 71 { sphereIndex[pairIndex][0]=i; 72 sphereIndex[pairIndex][1]=j; 73 springConst[p spring Const[pairInd airIndex] ex] = springC; springC; 74 springLength[ springLength[pairIn pairIndex] dex] = springL; springL; 75 pairIndex++; 76 77 } 70
78
void void Sphere Sphere::fo ::force rce1() 1() //Spri //Spring ng Force, Force, works works only only on attach attached ed sphere spheres, s, with damping { 80 doubl double e R, Rx, Rx, Ry, Ry, xFji xFji=0, =0, yFji yFji=0 =0, , delt deltaR aR; ; 81 79
13
int int i,j; i,j;
82 83
//su //sums ms a all ll o of f the the forc forces es ( (be bewa ware re o of f the the sign sign) ) for for all all pair pairs s for(int for(int k=0; k=0; k
84 85 86 87 88 89 90
//dete //determi rmine ne the relati relative ve positi position on vector vector Rx = ( xpos xpos[i [i] ] - xpos xpos[j [j] ] ); Ry = ( ypos ypos[i [i] ] - ypos ypos[j [j] ] );
91 92 93 94 95
R = sqrt sqrt(R (Rx x*Rx + Ry*Ry); delt deltaR aR = R - sprin springL gLen engt gth[ h[k] k]; ;
96 97 98
xFji xFji = spring springCons Const[k t[k] ]*deltaR*Rx/R;
99 100
xtotal xtotalF[j F[j] ] += xFji - dampin dampingCo gConst nst*(xvel[j (xvel[j] ] - xvel[i xvel[i]); ]);
101 102
xtotal xtotalF[i F[i] ] += -xFji -xFji - dampin dampingCo gConst nst*(xvel[i (xvel[i] ] - xvel[j xvel[j]); ]);
103 104
yFji yFji = spring springCons Const[k t[k] ]*deltaR*Ry/R;
105 106
ytotal ytotalF[j F[j] ] += yFji - dampin dampingCo gConst nst*(yvel[j (yvel[j] ] - yvel[i yvel[i]) ]) ;
107 108 109
ytotal ytotalF[i F[i] ] += -yFji -yFji - dampin dampingCo gConst nst*(yvel[i (yvel[i] ] - yvel[j yvel[j]) ]) ;
110 111 112
}
113 114
}
115
//========== //================= ============= ============ ============= ============= ============ =========== ===== GROUND CONSTRAINT CONSTRAINT checkGround() 117 void checkGround() 118 { for(in for(int t i=0; i=0; iyp >ypos[ os[i] i] <= 0) 121 { 122 pball->yvel[ pball->yvel[i] i] = fabs(pball-> fabs(pball->yvel[i yvel[i]); ]); 123 pballpball->yp >ypos[ os[i] i] = 0; 124 } 125 } 126 } 127 116
128 129 130
//============= //=================== ============ ============= ============= ============ =========== ===== MOVER, PLOTTER, PLOTTER, MOVIE MAKER
131
void void horiOs horiOscill cillato ator(d r(doub ouble le times) times) //osci //oscilla llate te the ground ground floor floor 133 { static static double double xi = pbal pball-> l->xpos xpos[0] [0]; ; 134 132
14
static static double double xj = pbal pball-> l->xpos xpos[1] [1]; ; static static double double xk = pbal pball-> l->xpos xpos[2] [2]; ; //stat //static ic d doub ouble le x xl l = pballpball->xpo >xpos[3 s[3]; ];
135 136 137 138
//stab //stab the sphere sphere pballpball->ypo >ypos[0 s[0] ] = 0.0; 0.0; pballpball->ypo >ypos[1 s[1] ] = 0.0; 0.0; pballpball->ypo >ypos[2 s[2] ] = 0.0; 0.0;
139 140 141 142 143
pbal pballl->xp >xpos os[0 [0] ] = -4 + ampl amplitu itude de*sin(omega*times); pbal pballl->xp >xpos os[1 [1] ] = 0 + ampl amplit itude ude*sin(omega*times); pbal pballl->xp >xpos os[2 [2] ] = 4 + ampl amplit itude ude*sin(omega*times);
144 145 146 147
}
148
void plotToMonitor plotToMonitor() () 150 { pspherePlot->redata(); 151 pspherePlot->recommand(); 152 149
153
for(in for(int t l=0; l=0; ldata <xpos[l]<<" "<ypos[l]<
154 155 156 157 158 159 160
for(in for(int t k=0; k=0; k
161 162 163 164 165
166
pspherePlot->command <<"set arrow "<xp "<xpos[i]<< os[i]<<","< ball->ypos[i ypos[i ]<<" ]<<" to " <xpos <xpos[j]<<" [j]<<","<yp ll->ypos[j]< os[j]<<" <" nohead lw 6 ;"<
167
168
}
169 170 171
pspherePlot->command <<"set <<"set xzeroa xzeroaxis xis 7;set 7;set size square square;" ;" <<"set border;" border;" <<"plot [-11:11][-2: [-11:11][-2:"<
pspherePlot->plot();
172 173 174 175 176 177 178 179
}
180
voi void d plot plotTo ToPN PNG( G() ) //mo //move ve all all the the ball balls s and and put put it on a fram frame e 182 { stat static ic int imag imageN eNum um = 0.0 0.0; ; 183 pspherePlot->redata(); 184 pspherePlot->recommand(); 185 181
186 187
15
//plot //plot the motio motions ns and apply apply the perio periodic dic bound boundary ary algori algorithm thm for( for(in int t i =0; =0; idata <xpos[i]<<" "<ypos[i]<
188 189 190 191 192 193 194
for(in for(int t k=0; k=0; k
195 196 197 198 199
200
pspherePlot->command <<"set arrow "<xp "<xpos[i]<< os[i]<<","< ball->ypos[i ypos[i ]<<" ]<<" to " <xpos <xpos[j]<<" [j]<<","<yp ll->ypos[j]< os[j]<<" <" nohead lw 6 ;"<
201
202
}
203 204
205 206 207 208 209 210 211 212 213
pspherePlot->command <<"se <<"set t term term png png size size 480,7 480,750 50;s ;set et xzer xzeroa oaxi xis s 7;" <<"s <<"set et out out " <<"’images" <<"_"< <<"_"<
214 215
pspherePlot->plot();
imageNum++;
216 217 218
}
219
void void makeMovie() makeMovie() { string movieCommand= movieCommand= "C:\\Users\\u "C:\\Users\\user\\D ser\\Downloa ownloads\\ffm ds\\ffmpeg\\b peg\\bin\\ff in\\ffmpeg mpeg 222 qsca qscale le 5 -r "; movieCommand.append(movieSpeed); 223 movieC movieComma ommand. nd.app append end(" (" -b 9600 -i image\\_\ image\\_\%d. %d.png png quake.m quake.mp4 p4 \n"); \n"); 224 system(movieCommand.c_str()); 225 226 } 220 221
227
//double //double dampingRatioT dampingRatioTMD(int MD(int floor) //{ doub double le tot total alMa Mass ssBu Buil ildi ding ng = mass massPe PerS rSph pher ere e*floor; 230 // // doub do uble le mass ma ssBa Bar r = massD mas sDam amper per/t /tot otal alMa Mass ssBu Buil ildi ding ng; ; 231 double xi = 0; 232 // 233 // return ; 234 // 235 //} 228 229
236
double double meanShe meanShearF arForc orce(i e(int nt floor) floor) //meng //menghit hitung ung gaya gaya gunting gunting tiap tiap lantai lantai 238 { double f1,f2,f3,resu f1,f2,f3,result; lt; 239 237
16
f1 = pball-> pball->xac xacc[f c[floo loor r*3] * massPerSphere; f2 = pball-> pball->xac xacc[f c[floo loor r*3+1] * massPerSphere; f3 = pball-> pball->xac xacc[f c[floo loor r*3+2] * massPerSphere; result result = (f1+f2 (f1+f2+f3 +f3)/3 )/3.0; .0; return fabs(result); fabs(result);
240 241 242 243 244 245
}
246 247
double temp[numOfFlo temp[numOfFloors+1] ors+1]; ; void id shea shearF rForc orceP ePlo lot( t() ) //pl //plot ot the the shea shear r forc forces es on each each floo floor r 249 vo 250 { stat static ic int firs firstT tTim ime e = 1; 251 double tempValue; tempValue; 252 248
253
254
pshearPlot->redata();
255
if(fir if(firstTi stTime me != 0) { for(int for(int i=0; i=0; i<(num i<(numOfF OfFloor loors+1 s+1); ); i++) i++) { temp temp[i [i] ] = 0.0; 0.0; } firstTime=0; }
256 257 258 259 260 261 262 263 264
for(int for(int f=1;f
265 266 267 268
269 270
271 272 273
if(tempValue>temp[f]) { pshearPlot->data <
274
} else {
275
276 277
278
pshearPlot->data <
279
}
280 281
}
282 283 284 285
}
pshearPlot->closeData();
286
void plotShearToMo plotShearToMonitor( nitor() ) 288 { pspherePlot->recommand(); 289 287
290 291 292
pspherePlot->command <<"se <<"set t zero zeroax axis is 7 7; ; set set xlab xlabel el ’ ’Fl Floo oor r Numb Number er’; ’; s set et y yla labe bel l ’Max ’Maximu imum m Shear Shear Force Force (N)’;" (N)’;"
17
<<"set <<"set border border; ; set grid;" grid;" <<"plot <<"plot " <<"’she <<"’shearD arData ata.tx .txt’ t’ w l" <<";pau <<";pause se -1 \n"; \n";
293 294 295 296 297
298 299
pspherePlot->plot();
}
300 301
//========== //================= ============= ============ ============= ============= ========== ==== MAIN FUNCTION FUNCTION int main() main() 304 { //pl //plac ace e the y positi position on sphe spheres res 305 for( for(in int t f = 0; f<=n f<=num umOf OfRo Room oms s ; f++) f++) 306 { 307 for(i for(int nt r = 0; r<=n r<=num umOf OfFlo Floor ors; s; r++) r++) 308 { 309 pball->ypos[3*r+f] r+f] = 3.0 3.0*r; 310 } 311 } 312 302 303
313 314 315 316 317 318 319
for( for(in int t { }
r = 0; r<=n r<=num umOf OfFl Floo oors; rs; r++ r++) ) pball->xpos[r*3] = -4; pball->xpos[r*3+1] 3+1] = 0; pball->xpos[r*3+2] 3+2] = 4;
320 321 322 323 324 325 326
//wall //wall spring spring for( for(in int t r = 0; r
327 328 329 330 331 332 333
//floor //floor spring for( for(in int t r = 0; 0; r < num numOf OfFl Floo oors rs; ; r++) r++) { attachSpring(r*3, r*3+1, 4.0, wallSpringCon wallSpringConst); st); attachSpring(r*3+1, 3+1, r*3+2, 4.0, wallSpringCon wallSpringConst); st); }
334 335 336 337 338 339 340
//cross //cross for( for(in int t { }
spring r = 0; r < numO numOfF fFlo loor orss-1; 1; r++) r++) attachSpring(r*3,r*3+4, 5.0, shearSpringCo shearSpringConst); nst); attachSpring(r*3+1,r*3+5, 5.0, shearSpringCo shearSpringConst); nst);
341 342 343 344 345 346
for( for(in int t r = 0; r < numO numOfF fFlo loor orss-1; 1; r++) r++) { attachSpring(r*3+1,r*3+1+2, 5.0, shearSpringCo shearSpringConst); nst); attachSpring(r*3+2,r*3+2+2, 5.0, shearSpringCo shearSpringConst); nst); }
18
347 348
//atta //attach ch the mass mass for( for(in int t i = 0; 0; imass[i pball->mass[i] ] = massPerSphere massPerSphere; ; }
349 350 351 352 353 354
//Plot //Plot the init initial ial condit conditions ions to moni monitor tor plotToMonitor(); cout<<"calculating..."<
355 356 357
358
//Comp //Compute ute and and write write the the changi changing ng posit position ions s for(double time=0.0;time
359 360 361 362 363
364
pball->velVerletMove(timeStep);
365
if (shake==true (shake==true) ) horiOscillato horiOscillator(time r(time); );
366 367
if(torquePlot==true) shearForcePlot();
368 369 370 371
372 373
if((cr if((creat eateMo eMovie vie==t ==true) rue) and ((int) ((int) (time/t (time/time imeSte Step)% p)%25 25 == 0)) { plotToPNG(); }
374 375 376 377
378 379 380
//desc //describ ribe e the computa computatio tion n progre progress ss if( if( ((in ((int) t) (tim (time/ e/ti timeS meSte tep) p)) ) % 50 == 0) { system("cls"); cout<< cout<< (int) (int) (100 * time/ time/tim timeEn eEnds) ds) <<" % "; }
381
}
382
383 384
system("cls"); cout<<"100 cout<<"100 %"<
385 386 387 388 389
//Plot //Plot the final final condit condition ions s plotToMonitor(); cout<
390 391
if(createMovie==true) { makeMovie(); }
if(torquePlot==true) { plotShearToMonitor(); }
392 393 394 395 396 397 398 399 400
19
401 402
system("pause");
}
20
