You will need Aladdin 2.0 to run this problem.
/* * ===================================================================== * Plastic Analysis for a System of Two 1-DOF Spring Elements * * -- general matrix format for bi-linear material * -- energy calculation for balance checking * * Dynamic analysis using average acceleration method. * * Written By : Wane-Jang Lin October, 1997 * ===================================================================== */ /* plastic spring material, Bi-linear */ ks1 = 2.0 N/cm; ks2 = 1.5 N/cm; kt1 = 0.8 N/cm; kt2 = 0.5 N/cm; fy1 = 18 N; fy2 = 15 N; m1 = 1 kg; m2 = 0.5 kg; c1 = 2 N*sec/m; c2 = 1 N*sec/m; /* damping coefficient */ Ks = [ks1; ks2]; Kt = [kt1; kt2]; Fy = [fy1; fy2]; ey = Matrix([2,1]); for( ele=1; ele<=2; ele=ele+1 ) { ey[ele][1] = Fy[ele][1]/Ks[ele][1]; } /* Initial condition, unstressed */ P = [0 N; 0 N]; /* structure force */ p = [0 cm; 0 cm]; /* structure displacement */ PR = [0 N; 0 N]; /* structure resisting force */ PRe = [0 N; 0 N]; /* element resisting force */ Q_saved = [0 N; 0 N]; q_saved = [0 cm; 0 cm]; /* Initial flags and stresses and strains, initialed before any load step */ index = [0;0]; /* index for loading flag */ flag1 = [0;0]; /* yielding at load step k */ flag2 = [0;0]; /* pre_range at load step k */ flag3 = [0;0]; /* pre_load at load step k */ yielding = [0;0]; pre_range = [0;0]; pre_load = [0;0]; loading = [0;0]; sr = [0 N; 0 N]; er = [0 cm; 0 cm]; s0 = [0 N; 0 N]; e0 = [0 cm; 0 cm]; sr_saved = [0 N; 0 N]; er_saved = [0 cm; 0 cm]; s0_saved = [0 N; 0 N]; e0_saved = [0 cm; 0 cm]; sx_saved = [0 N; 0 N]; ex_saved = [0 cm; 0 cm]; /* Transformation matrix Lele, d_q = Lele*d_p[ele] for each element */ /* transform structural displacements d_p to element deformations d_q */ Rigid = [-1,1]; Transform = [1,0;0,1]; L = Rigid*Transform; /* Temporary matrix to store element tangent Ks at each step */ tangent = [ks1 ; ks2]; /* Force interpolation matrices b(x), D(x) = b(x)*Q */ /* relate section force D(x) with element force Q */ bx = 1; /* assemble initial structure tangent stiffness matrix BigK */ BigK = [ ks1+ks2, -ks2; -ks2, ks2 ]; /* assemble mass matrix and damping matrix */ BigM = [ m1, 0 kg; 0 kg, m2 ]; BigC = [ c1+c2, -c2; -c2, c2 ]; M_inv = Inverse(BigM); /* initial time = 0 sec conditions */ total_step = 400; dt = 0.01 sec; time = 0 sec; /* Setup initial velocity and acceleration */ vel = p/(1 sec); accel = vel/(1 sec); new_p = p; new_vel = vel; new_accel = accel; /* Setup initial internal force and damping force */ Fs = P; Fd = P; P_old = P; Fs_old = P; Fd_old = P; /* Allocate the matrices for storing output history */ /* column[1] : external applied load at end node (2) */ /* column[2] : total elogation at end node (2) = [3]+[4] */ /* column[3] : element elogation 1, node 1 */ /* column[4] : element elogation 2, node 2 */ /* column[5] : element force 1, node 1 */ /* column[6] : element force 2, node 2 */ /* column[7] : dynamic natrual period 1 (T1) for each step */ /* column[8] : dynamic natrual period 2 (T2) for each step */ result = ColumnUnits( Matrix([total_step+1,8]), [N,cm,cm,cm,N,N,sec,sec] ); result[1][1] = P[2][1]; result[1][2] = p[2][1]; result[1][3] = p[1][1]; result[1][4] = p[2][1] - p[1][1]; result[1][5] = P[2][1]; result[1][6] = P[2][1]; /* Compute eigen problem */ eigen = Eigen(BigK, BigM, [2]); eigenvalue = Eigenvalue(eigen); result[1][7] = 2*PI/sqrt( eigenvalue[1][1] ); result[1][8] = 2*PI/sqrt( eigenvalue[2][1] ); /* column[1] : internal strain energy for element 1 */ /* column[2] : elastic strain energy for element 1 */ /* column[3] : internal strain energy for element 2 */ /* column[4] : elastic strain energy for element 2 */ element_energy = ColumnUnits( Matrix([total_step+1,4]), [Jou] ); element_energy[1][1] = 0 Jou; element_energy[1][2] = 0 Jou; element_energy[1][3] = 0 Jou; element_energy[1][4] = 0 Jou; /* increase external load, structure determination */ for ( k=1; k<=total_step ; k=k+1 ) { time = time + dt; if( k<=20 ) { d_P = [0 N; 1 N]; } if( k>20 && k<=60 ) { d_P = -[0 N; 1 N]; } if( k>60 && k<=105 ) { d_P = [0 N; 1 N]; } if( k>105 && k<=155 ) { d_P = -[0 N; 1 N]; } if( k>155 && k<=190 ) { d_P = [0 N; 1 N]; } if( k>190 && k<=200 ) { d_P = -[0 N; 1 N]; } if( k>200 ) { d_P = [0 N; 0 N]; } P = P + d_P; /* Compute effective incremental loading */ dPeff = d_P + ((4/dt)*BigM + 2*BigC)*vel + 2*BigM*accel; element_energy[k+1][1] = element_energy[k][1]; element_energy[k+1][3] = element_energy[k][3]; index[1][1] = 1; index[2][1] = 1; /* Compute effective stiffness from tangent stiffness */ Keff = BigK + (2/dt)*BigC + (4/dt/dt)*BigM; /* Solve for d_displacement, d_velocity */ d_p = Solve( Keff, dPeff ); d_v = (2/dt)*d_p - 2*vel; /* Compute displacement, velocity */ new_p = p + d_p; new_vel = vel + d_v; /* State determination for each element -- check material yielding */ /* and compute new stiffness */ for( ele = 1; ele <= 2; ele = ele + 1 ) { if( ele==1 ) { d_pe = [ 0 m; d_p[1][1] ]; } if( ele==2 ) { d_pe = [ d_p[1][1]; d_p[2][1] ]; } Q = Q_saved[ele][1]; /* retrieve values from (j-1) */ q = q_saved[ele][1]; Dx = bx*Q; dx = bx*q; K = tangent[ele][1]; kx = tangent[ele][1]; fx = 1/kx; rx = 0 cm; DUx = 1E+7 N; d_q = QuanCast(L*d_pe); /* element deformation increment */ q = q + d_q; /* iterate for convervence of element displacements */ while( abs(DUx) > 0.00001 N ) { d_Q = K*d_q; /* element force increment */ Q = Q + d_Q; /* update element force */ /* determine the section force increments */ /* repeat for all integration points of the element */ d_Dx = bx*d_Q; /* section force increment */ d_dx = rx + fx*d_Dx; /* section deformation increment */ Dx = Dx + d_Dx; /* update section force */ dx = dx + d_dx; /* update section deformation vector */ /* get new section tangent flexibility f(x) from new d(x) */ /* and section resisting force DR(x) from material properties */ /* set flags and stresses and strains for each i-th iteration */ if( index[ele][1] == 1 ) { if( d_dx > 0 cm ) { loading[ele][1] = 1; } if( d_dx < 0 cm ) { loading[ele][1] = -1; } if( d_dx == 0 cm ) { loading[ele][1] = pre_load[ele][1]; } index[ele][1] = index[ele][1] + 1; } yielding[ele][1] = flag1[ele][1]; pre_range[ele][1] = flag2[ele][1]; pre_load[ele][1] = flag3[ele][1]; sr[ele][1] = sr_saved[ele][1]; er[ele][1] = er_saved[ele][1]; s0[ele][1] = s0_saved[ele][1]; e0[ele][1] = e0_saved[ele][1]; /* material is still in elastic range, does not have any plastic residual */ if( yielding[ele][1] == 0 ) then { if( abs(dx) <= ey[ele][1] ) then { kx = Ks[ele][1]; sx = 0 N; ex = 0 cm; yielding[ele][1] = 0; pre_range[ele][1] = 0; } else { kx = Kt[ele][1]; s0[ele][1] = Fy[ele][1]*loading[ele][1]; e0[ele][1] = ey[ele][1]*loading[ele][1]; sx = s0[ele][1]; ex = e0[ele][1]; yielding[ele][1] = 1; pre_range[ele][1] = 1; } } else { /* plastic residual occurs , yielding[ele][1] = 1 */ if( pre_load[ele][1] != loading[ele][1] ) then { if( pre_range[ele][1] == 1 ) then { sr[ele][1] = sx_saved[ele][1]; er[ele][1] = ex_saved[ele][1]; kx = Ks[ele][1]; sx = sr[ele][1]; ex = er[ele][1]; pre_range[ele][1] = 0; } else { /* pre_range[ele][1] = 0 */ if( loading[ele][1]*dx <= loading[ele][1]*er[ele][1] ) then { kx = Ks[ele][1]; sx = sr[ele][1]; ex = er[ele][1]; pre_range[ele][1] = 0; } else { if( loading[ele][1]*ex_saved[ele][1] >= loading[ele][1]*er[ele][1] ) then { kx = Ks[ele][1]; sx = sr[ele][1]; ex = er[ele][1]; pre_range[ele][1] = 0; } else { kx = Kt[ele][1]; sx = sr[ele][1]; ex = er[ele][1]; pre_range[ele][1] = 1; } } } } else { /* pre_load[ele][1] == loading[ele][1] */ if( pre_range[ele][1] == 1 ) then { kx = Kt[ele][1]; sx = s0[ele][1]; ex = e0[ele][1]; pre_range[ele][1] = 1; } else { /* pre_range[ele][1] = 0 */ if( loading[ele][1]*ex_saved[ele][1] <= loading[ele][1]*er[ele][1] ) then { /* line 2 */ if( loading[ele][1]*dx <= loading[ele][1]*er[ele][1] ) then { /* line 2 */ kx = Ks[ele][1]; sx = sr[ele][1]; ex = er[ele][1]; pre_range[ele][1] = 0; } else { /* line 2 -> line 1 */ kx = Kt[ele][1]; sx = sr[ele][1]; ex = er[ele][1]; pre_range[ele][1] = 1; } } else { /* line 4 */ if( abs(dx-er[ele][1]) <= 2*ey[ele][1] ) then { /* line 4 */ kx = Ks[ele][1]; sx = sr[ele][1]; ex = er[ele][1]; pre_range[ele][1] = 0; } else { /* line 4 -> line 1 */ s0[ele][1] = sr[ele][1] + loading[ele][1]*2*Fy[ele][1]; e0[ele][1] = er[ele][1] + loading[ele][1]*2*ey[ele][1]; kx = Kt[ele][1]; sx = s0[ele][1]; ex = e0[ele][1]; pre_range[ele][1] = 1; } } /* end of line 4 */ } } } /* calculate stress and flexibility */ fx = 1/kx; DRx = sx + kx*(dx-ex); DUx = Dx - DRx; /* section unbalanced force */ rx = fx*DUx; /* section residual deformation */ /* finish for all integration points of the element */ /* update the element flexibility and stiffness matrices */ F = bx*fx*bx; K = 1/F; /* check for element convergence */ s = bx*rx; /* element residual deformation */ d_q = -s; } /* j, while loop to check element convergence */ /* Energy calculations for each element ele */ element_energy[k+1][2*ele-1] = element_energy[k+1][2*ele-1] + 0.5*(Q_saved[ele][1]+Q)*(q-q_saved[ele][1]); if( abs(Q) > 2*Fy[ele][1] ) then { if( kx == Ks[ele][1] ) then { delta_Q = abs(sr[ele][1]) - 2*Fy[ele][1]; } else { delta_Q = abs(Q) - 2*Fy[ele][1]; } delta_x = delta_Q*(1/Kt[ele][1]-1/Ks[ele][1]); element_energy[k+1][2*ele] = 0.5*Q*Q/Ks[ele][1] + 0.5*delta_x*delta_Q; } else { if( abs(sr[ele][1]) > 2*Fy[ele][1] ) then { if( kx == Ks[ele][1] ) then { delta_Q = abs(sr[ele][1]) - 2*Fy[ele][1]; delta_x = delta_Q*(1/Kt[ele][1]-1/Ks[ele][1]); element_energy[k+1][2*ele] = 0.5*Q*Q/Ks[ele][1] + 0.5*delta_x*delta_Q; } else { if(((Q>0N)&&(sr[ele][1]>0N))||((Q<0N)&&(sr[ele][1]<0N))) then { element_energy[k+1][2*ele] = 0.5*Q*Q/Kt[ele][1]; } else { element_energy[k+1][2*ele] = 0.5*Q*Q/Ks[ele][1]; } } } else { element_energy[k+1][2*ele] = 0.5*Q*Q/Ks[ele][1]; } } Q_saved[ele][1] = Q; q_saved[ele][1] = q; tangent[ele][1] = K; /* element stiffness Ke=LT*K*L=[K,-K;-K,K] */ PRe[ele][1] = Q; /* element resistant force PRe=LT*Q=[-Q;Q] */ } /* ele */ /* assemble new structure stiffness */ BigK[1][1] = tangent[1][1] + tangent[2][1]; BigK[1][2] = -tangent[2][1]; BigK[2][1] = -tangent[2][1]; BigK[2][2] = tangent[2][1]; /* assemble new internal load Fs, damping load Fd, and acceleration */ Fs[1][1] = PRe[1][1] - PRe[2][1]; Fs[2][1] = PRe[2][1]; Fd = BigC*new_vel; new_accel = M_inv*( P-Fs-Fd ); /* Store analysis results */ result[k+1][1] = P[2][1]; result[k+1][2] = p[2][1]; result[k+1][3] = p[1][1]; result[k+1][4] = p[2][1] - p[1][1]; result[k+1][5] = Fs[1][1] + Fs[2][1]; result[k+1][6] = Fs[2][1]; /* Compute eigen problem */ eigen = Eigen(BigK, BigM, [2]); eigenvalue = Eigenvalue(eigen); result[k+1][7] = 2*PI/sqrt( eigenvalue[1][1] ); result[k+1][8] = 2*PI/sqrt( eigenvalue[2][1] ); /* Updating history: save flags and reversal points for each load step k */ for( ele=1 ; ele<=2 ; ele=ele+1 ) { pre_load[ele][1] = loading[ele][1]; flag1[ele][1] = yielding[ele][1]; flag2[ele][1] = pre_range[ele][1]; flag3[ele][1] = pre_load[ele][1]; sr_saved[ele][1] = sr[ele][1]; er_saved[ele][1] = er[ele][1]; s0_saved[ele][1] = s0[ele][1]; e0_saved[ele][1] = e0[ele][1]; sx_saved[ele][1] = Q_saved[ele][1]; ex_saved[ele][1] = q_saved[ele][1]; } /* Update information for this step */ P_old = P; Fs_old = Fs; Fd_old = Fd; p = new_p; vel = new_vel; accel = new_accel; } /* k in for() loop, increase to next time step */ PrintMatrix(result); PrintMatrix(element_energy); quit;