angry bird.zip_bird_built6xi_小游戏

function [] = angry bird () %*********************************************************************** %- game：angry bird * %- Programed by Shu-Hui Jiang, 2014 * %- Revised by Xing-Fei Yuan, 2014 * %- Proofed by Edward C. Ting, 2011 * %*********************************************************************** clear all; %- Open file for dumping tecplot file file_plt = fopen('tec.plt','w'); file_Fc = fopen('tec.force','w'); %- Preprocess ET = 5.0; %- Termination time h = 5.e-4; %- Size of time step Tpr = 1.e-2; %- Period of dumping history data Tani = 1.e-1; %- Period of dumping tecplot data Tscr = 1.e-1; %- Period of dumping screen data zeta = 0.00; %- damping factor A = 1.0; %- Area I = 1.0; %- 2nd moment of cross section* Young's modulus E1 = 1.e4; %- Young's modulus of bird E2 = 1.e5; %- Young's modulus of frame G = 9.81; %- Gravity acceleration Mcr = 500.; %- Crack moment kc1 = 1.e5; %- Contact spring stiffness1（to element） kc2 = 1.e6; %- Contact spring stiffness2（to boundary） mu = 0.0; %- Friction Coef. of frame surface D = 0.1; %- Depth of beam section rho = 0.5; %- mass density nel = 10; %- No. of element npt = 11; %- No. of point nfrm =0; % frame no. for movie ifrm = 1.0; %- fragmentation index eL1 =1; %element length of bird eL2 =10; LL(1)=1; LL(2)=20; theta=pi/6; d2=D/2; EA1 = A*E1; EA2 = A*E2; EI1 = E1*I; EI2 = E2*I; %- allocate array memory x0 = zeros(2*nel,2); xn = zeros(2*nel,2); xn1 = zeros(2*nel,2); xn_1 = zeros(2*nel,2); dn = zeros(2*nel,2); dn1 = zeros(2*nel,2); rn = zeros(2*nel,1); rn1 = zeros(2*nel,1); rn_1 = zeros(2*nel,1); beta = zeros(2*nel,1); mass = zeros(1,2*nel); emas = zeros(1,nel); nen = zeros(nel,2); %- Elem. connectivity ned = zeros(3,2); %- boundary connectivity nje = zeros(npt,1); f = zeros(2*nel,2);m=zeros(2*nel,1);ef = zeros(nel,3);P = zeros(2*nel,2); %- Constants for time integration c1 = 1.d0 + zeta*h/2.d0; c2 = 1.d0 - zeta*h/2.d0; h2 = h*h; %- boundary a-b-c-d xa(1,1)=0;xa(1,2)=30; xa(2,1)=0;xa(2,2)=0; xa(3,1)=40;xa(3,2)=0; xa(4,1)=40;xa(4,2)=30; %- initial nodal positions CS = cos(theta); SN = sin(theta); x0(1,1) = LL(1); x0(1,2) = d2; %- pt 1 x0(2,1) = LL(1)+ eL1 ; x0(2,2) = d2; %- pt 2 x0(3,1) = x0(2,1)+ eL1*SN ; x0(3,2) = x0(2,2) + eL1*CS; %- pt 3 x0(4,1) = x0(2,1); x0(4,2) = x0(2,2) + 2*eL1*CS; %- pt 4 x0(5,1) = x0(1,1); x0(5,2) = x0(4,2); %- pt 5 x0(6,1) = x0(1,1)- eL1*SN ; x0(6,2)= x0(3,2); %- pt 6 x0(7,1)=LL(2);x0(7,2)=0; x0(8,1)=LL(2);x0(8,2)=eL2; x0(9,1)=LL(2)+eL2;x0(9,2)=0; x0(10,1)=LL(2)+eL2;x0(10,2)=eL2; x0(11,1)=LL(2)+eL2;x0(11,2)=2*eL2; %- element conectivity (for tecplot) nen(1,1) = 1; nen(1,2) = 2; nen(2,1) = 2; nen(2,2) = 3; nen(3,1) = 3; nen(3,2) = 4; nen(4,1) = 4; nen(4,2) = 5; nen(5,1) = 5; nen(5,2) = 6; nen(6,1) = 6; nen(6,2) = 1; nen(7,1) = 7; nen(7,2) = 8; nen(8,1) = 8; nen(8,2) = 10; nen(9,1) = 9; nen(9,2) = 10; nen(10,1) = 10; nen(10,2) = 11; ned(1,1) = 12; ned(1,2) = 13; ned(2,1) = 13; ned(2,2) = 14; ned(3,1) = 14; ned(3,2) = 15; %- jointed element at each pt nje=zeros(npt,1); for i = 1:1:nel n1 = nen(i,1); n2 = nen(i,2); nje(n1) = nje(n1) + 1; nje(n2) = nje(n2) + 1; end %- Dump initial configuration fprintf(file_plt,'ZONE N = %5d E = %5d F = FEPOINT ET = TRIANGLE\n',npt,nel); for i=1:1:npt fprintf(file_plt,' %12.4E %12.4E\n', x0(i,1), x0(i,2)); end for i=1:1:4 fprintf(file_plt,' %12.4E %12.4E\n', xa(i,1), xa(i,2)); end for i=1:1:nel fprintf(file_plt,' %4d %4d %4d\n',nen(i,1),nen(i,2),nen(i,2)); end for i=1:1:3 fprintf(file_plt,' %4d %4d %4d\n',ned(i,1),ned(i,2),ned(i,2)); end %- point mass emas(1:6) = rho*A*eL1; emas(7:nel)=rho*A*eL2; [mass,rmas]=ptmass(emas,nen,A,I,nel); %- initial conditions (n = 0) vn=zeros(2*nel,2); %- initial velocity for i=1:6 vn(i,1) = 12; vn(i,2) = 12; end rvn=zeros(2*nel,1); %- initial ang. velocity for i=1:6 rvn(i,1) = 0; end dn = zeros(2*nel,2); rn = zeros(2*nel,1); xn = x0 + dn; %- Initialized time variables time = 0.0; %- Analysis time hc = 0.0; %- Time counter ha = 0.0; %- Time counter hs = 0.0; %- Time counter %- Initial internal nodal force f=zeros(2*nel,2); %- internal force of point m=zeros(2*nel,1); %- internal moment of point %- Initial element internal force ef = zeros(nel,3); %- Initial external force gra=-G*mass; gra=gra'; P=zeros(2*nel,2); P(:,2)=gra; %- Evaluate x_1 for i=1:1:npt xn_1(i,:) = xn(i,:) - (h + 0.5d0*zeta*h2)*vn(i,:) + 0.5d0*h2*(P(i,:) + f(i,:))/mass(i); rn_1(i) = rn(i) - (h + 0.5d0*zeta*h2)*rvn(i) + 0.5d0*h2*(+ m(i))/rmas(i); end %- impose restraints at the ends xn_1(7,:)=x0(7,:); xn_1(9,:)=x0(9,:); xn_1(11,:)=x0(11,:); %- Dump initial data n = 1; time_dmp(n) = time; dp(n,:) = dn(1,:); %-以point1的位移追踪 vp(n,:) = vn(1,:); %- Time integration for time = 0.0: h: ET %- Central difference stencile for i=1:1:npt xn1(i,:) = (2.0*xn(i,:) - c2*xn_1(i,:) + h2*(P(i,:) + f(i,:))/mass(i))/c1; vn(i,:) = (xn1(i,:) - xn_1(i,:))/(2.0*h); rn1(i) = (2.0*rn(i) - c2*rn_1(i) + h2*(+ m(i))/rmas(i))/c1; beta(i) = rn1(i) - rn(i); dn1(i,:) = xn1(i,:) - x0(i,:); end %- impose restraints at the ends xn1(7,:)=x0(7,:); xn1(9,:)=x0(9,:); xn1(11,:)=x0(11,:); dn1(7,:) = xn1(7,:) - x0(7,:); beta(7) = rn1(7) - rn(7); dn1(9,:) = xn1(9,:) - x0(9,:); beta(9) = rn1(9) - rn(9); dn1(11,:) = xn1(11,:) - x0(11,:); beta(11) = rn1(11) - rn(11); %- Internal forces fij_i=zeros(nel,2); fij_j=zeros(nel,2); f=zeros(2*nel,2); m=zeros(2*nel,1); for i=1:1:nel n1=nen(i,1);n2=nen(i,2); if i<=6 EA=EA1;EI=EI1; else EA=EA2;EI=EI2; end [fij_i(i,:),fij_j(i,:),ef(i,:)]=beam2Df(xn(n1,:),xn1(n1,:),xn(n2,:),xn1(n2,:),beta(n1),beta(n2),EA,EI,ef(i,:)); %- Check failure if(i >= 7 && i <= 10) if(abs(ef(i,1)) >= Mcr && nje(n1) > 1) ifrm = 0.; npt = npt + 1; nen(i,1) = npt; xn1(npt,:) = xn1(n1,:); xn(npt,:) = xn(n1,:); rn1(npt) = rn1(n1); rn(npt) = rn(n1); nje(n1) = nje(n1) - 1; nje(npt) = 1; end if(abs(ef(i,2)) >= Mcr && nje(n2) > 1) ifrm = 0.; npt = npt + 1; nen(i,2) = npt; xn1(npt,:) = xn1(n2,:); xn(npt,:) = xn(n2,:); rn1(npt) = rn1(n2); rn(npt) = rn(n2); nje(n2) = nje(n2) - 1; nje(npt) = 1; end end n1 = nen(i,1); n2 = nen(i,2); %- Update point mass [mass,rmas]=ptmass(emas,nen,A,I,nel); %- Assemble force & moment to each pt f(n1,:)=f(n1,:)-fij_i(i,:); f(n2,:)=f(n2,:)-fij_j(i,:); m(n1)=m(n1)-ef(i,1); m(n2)=m(n2)-ef(i,2); end %- External force @ time + h gra=-G*mass; gra=gra'; P=zeros(2*nel,2); P(:,2)=gra; %- Contact force evaluation for i=1:1:npt for j=1:1:4 n1=nen(6+j,1);n2=nen(6+j,2); x1=xn1(n1,:); x2=xn1(n2,:); [Fc(i,:),Fcei(6+j,:),Fcej(6+j,:)] = contactf(x1,x2,xn1(i,:),vn(i,:),kc1,mu,d2); P(i,:) = P(i,:) + Fc(i,:); %- 接触后bird上接触点受到的力 P(n1,:) = P(n1,:) + Fcei(6+j,:); %- 接触后bird上接触点受到的力 P(n2,:) = P(n2,:) + Fcej(6+j,:); %- 接触后bird上接触点受到的力 end end for i=1:1:npt