Input file for Nonlinear Mass-Spring System
Input file for Nonlinear Mass-Spring System
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;