/*
* =================================================================
* Nonlinear Analysis of Cantilever with Shell Finite Element
*
* -- Illustrate application of BFGS algorithm to solution of
* nonlinear equations.
*
* Written By : X.G. Chen July-August 1994
* =================================================================
*/
print "*** DEFINE PROBLEM SPECIFIC PARAMETERS \n\n";
NDimension = 3;
NDofPerNode = 5;
MaxNodesPerElement = 8;
ThicknessIntegPts = 4;
STOL = 0.5;
EPS = 1.e-5;
IMAX = 20; /* Maximum number of iteration for BFGS */
LMAX = 50; /* Maximum number of load steps */
Ef = 1E-1;
Ee = 1E-3;
COUNTER = 7;
restart_counter = 0;
SetUnitsOff;
StartMesh();
/* [a] : Generate Finite Element Mesh for Cantilever Beam */
L = 10.0; /* Length of Beam */
b = 1.0; /* Width of Beam */
t = 1.0; /* Thickness of Beam */
print "*** Generate Grid of Nodes for F.E. Model \n\n";
z = 0 ;
dL = 2;
n2 = L/dL + 1;
NN = n2;
n1 = 3*n2 + 2*(n2-1);
for ( i = 1; i <= n2; i = i + 1) {
for (j = 1; j <= 3; j = j + 1) {
node = j + 5*(i-1);
x = dL*(i-1);
y = (b/2)*(j-2);
AddNode(node, [x, y, z]);
}
}
for ( i = 1; i < n2; i = i + 1) {
node = 3 + 1 + 5*(i-1);
x = (i-1)*dL + dL/2;
y = -b/2;
AddNode(node, [x, y, z]);
node = 3 + 2 + 5*(i-1);
x = (i-1)*dL + dL/2;
y = b/2;
AddNode(node, [x, y, z]);
}
print "*** ATTACH ELEMENTS TO GRID OF NODES \n\n";
for ( i = 1; i < n2; i = i + 1) {
elmtno = i;
a = 1 + 5*(i-1);
b = 6 + 5*(i-1);
c = 8 + 5*(i-1);
d = 3 + 5*(i-1);
e = 4 + 5*(i-1);
f = 7 + 5*(i-1);
gg = 5 + 5*(i-1);
h = 2 + 5*(i-1);
node_connec = [a, b, c, d, e, f, gg, h];
AddElmt(elmtno, node_connec, "name_of_elmt_attr");
}
/*
* [b] : Define Element, Section and Material Properties
*
* Note : beta = 1 for Isotropic strain hardening
*/
ElementAttr("name_of_elmt_attr") { type = "SHELL_8N";
section = "rectangular";
material = "ELASTIC_PLASTIC";
}
MaterialAttr("ELASTIC_PLASTIC") { density = 1.0;
poisson = 0.2;
yield = 36000;
E = 2.9E+7;
Et = 2.9E+6;
beta = 1.0;
type = "Ramberg-Osgood";
n = 6;
alpha = 3/7;
}
SectionAttr("rectangular") { width = 1;
thickness = 1;
}
/* [c] : Setup boundary conditions */
dx = 1; dy = 1; dz = 1;
rx = 1; ry = 1; rz = 1;
bc = [dx,dy,dz,rx,ry,rz];
FixNode(1, bc);
FixNode(2, bc);
FixNode(3, bc);
/* [d] : Add Point Nodal Loads */
Fx = 0; Fy = 0 ; Fz = 300;
Mx = 0; My = 0 ;
NodeLoad( n1-2, [ Fx, Fy, Fz, Mx, My]);
NodeLoad( n1-1, [ Fx, Fy, Fz, Mx, My]);
NodeLoad( n1, [ Fx, Fy, Fz, Mx, My]);
/* [e] : Compile and Print Finite Element Mesh */
EndMesh();
PrintMesh();
/* [f] : Generate Load Factor Matrix for Cyclic Loading */
Load_Factor = Matrix([LMAX,1]);
for (k = 1; k <= LMAX; k = k +1) {
if(k <= 10) then {
Load_Factor[k][1] = 1.0 + 0.05*k;
} else {
if(k > 10 && k <= 15) then {
Load_Factor[k][1] = 1.5 - 0.5*(k-10);
} else {
if(k <= 28) then {
Load_Factor[k][1] = -1.0 - 0.05*(k-15);
} else {
Load_Factor[k][1] = -1.65 + 0.5*(k-28);
}
}
}
}
/* [g] : Compute Mass and Stiffness Matrices */
print "\n*** Compute and Print Stiffness, Mass, and External Load Matrices\n\n";
stiff = Stiff();
Dimen = Dimension(stiff);
size = Dimen[1][1];
lu = Decompose(stiff);
eload0 = ExternalLoad();
eload = eload0;
displ = Substitution(lu, eload);
print " ----INITAIL DISPL ----\n";
print displ[size-2][1], "\n";
iload = InternalLoad(displ);
UpdateResponse();
/* [h] : Nonlinear Analysis with BFGS Method */
K_inver = Inverse(stiff);
I = Diag([size, 1]);
ii = 1;
jj = 1;
/*
* [i] : Compute Solution to next load step K(x)*delta_displ = R-F = R1
*/
LMAX = 34;
for (jj = 2; jj <= LMAX; jj = jj + 1) {
print "jj = ", jj, "\n";
eload = Load_Factor[jj][1]*eload0;
R2 = eload - iload;
R = R2;
R1 = R2;
Iter_no = 0;
/* [i.1] : Compute Solution to next load step K(x)*delta_displ = R-F = R1
for (ii = 1; ii <= IMAX; ii = ii + 1) {
print "ii = ", ii, "\n";
Iter_no = Iter_no + 1;
beta = 1.0;
delta_displ = K_inver*R2;
iload = InternalLoad(displ, delta_displ); /* beta = 1 */
R2 = eload - iload;
/* [i.2] : Line Search */
displ_temp = displ + delta_displ;
force_crt1 = L2Norm(R2); /* Force Criterion */
force_crt2 = L2Norm(R)*Ef;
print "jj = ", jj, "\n";
print "ii = ", ii, "\n";
print " force_crt1 = ", force_crt1;
print " force_crt2 = ", force_crt2, "\n";
temp2 = QuanCast(Trans(delta_displ)*R2);
energy_crt1 = abs(temp2); /* Energy Criterion */
if(ii == 1) {
temp1 = QuanCast(Trans(delta_displ)*R);
energy_crt2 = abs(temp1*Ee);
}
print " energy_crt1 = ", energy_crt1;
print " energy_crt2 = ", energy_crt2, "\n";
if( (force_crt1 <= force_crt2) && (energy_crt1 < energy_crt2)) then {
print beta, "\n";
displ = displ_temp;
UpdateResponse();
break;
} else {
temp1 = QuanCast(Trans(delta_displ)*R1);
temp2 = QuanCast(Trans(delta_displ)*R2);
counter = 0;
print " temp2 = ", temp2;
print " temp1 = ", temp1;
print " temp1*STOL + EPS = ", temp1*STOL + EPS, "\n";
while (temp2 > temp1*STOL + EPS) {
print " temp2 = ", temp2;
print " temp1 = ", temp1;
print " temp1*STOL + EPS = ", temp1*STOL + EPS, "\n";
counter = counter + 1;
print "\n ================================================\n";
print " counter = ", counter;
print "\n ================================================\n";
if(counter > COUNTER) then {
print "\n******* Too much iterations for line search \n";
break;
} else {
beta = beta/2.0;
delta_displ_temp = beta*delta_displ;
displ_temp = displ + delta_displ_temp;
iload = InternalLoad(displ, delta_displ_temp);
print " counter = ", counter, "\n";
PrintMatrix(iload);
R2 = eload - iload;
temp2 = QuanCast(Trans(delta_displ)*R2);
}
}
displ = displ_temp;
/* [i.3] : Restart for Failed Line Search */
if( counter > COUNTER ) then {
restart_counter = restart_counter + 1;
if( restart_counter > 3) {
print "****** Too many restarts \n";
break;
}
print " \n **** Restart at new initial Value \n";
PrintMatrix(eload, iload);
ii = 1;
print " \n **** Restart flag-1 \n";
K = Stiff();
print " \n **** Restart flag-2 \n";
K_inver = Inverse(K);
print " \n **** Restart flag-3 \n";
R1 = R2;
print " \n **** Restart flag-4 \n";
} else {
/* [i.4] : BFGS Algorithm */
gamma = R1-R2;
tem1 = QuanCast(Trans(delta_displ)*R1)*beta;
tem2 = QuanCast(Trans(delta_displ)*gamma);
COND_NO = sqrt(tem2/tem1); /* condition number */
V = COND_NO*beta*R1 - gamma;
W = delta_displ/tem2;
A = V*Trans(W);
A = I + A;
K_inver = Trans(A)*K_inver*A;
print "tem1 = ", tem1, "tem2 = ", tem2, "\n";
/* [i.5] : Test Convergence Criteria */
if (COND_NO > 1E+5) {
print" condition number = ", COND_NO, " \n";
print" Large condition number: ==> Jocobian matrix inverse nearly singular\n";
break;
}
R1 = R2;
}
}
UpdateResponse();
}
print" Results: \n --------------\n";
print" Iteration Number =", Iter_no, "\n\n";
print displ[size-2][1], "\n";
}
In part [a] we specify that this will be a two-dimensional analysis. The maximum number of degrees of freedom per node will be three, and the maximum number of nodes per element will be two. The parameters NDimension, NDofPerNode, and MaxNodesPerElement are used by ALADDIN to assess memory requirements for the problem storage and solution.
We use a nested for() loop and a single if() statement to generate the planar layout of 24 finite element nodes. Before the boundary conditions are applied, the structure has 72 degrees of freedom. In Section [f] we apply full-fixity to each column at the foundation level -- this reduces degrees of freedom from 72 to 60.
Return to Home Page