Potential bugs in Shell element

[elhaddad]

Hello Everyone,

I was having problems running linear dynamic analysis when my Opensees model has shell elements with non-zero rho value and subjected to ground excitation.
I decided to debug the code and ended up finding at least one bug.
I am posting them here to see if someone else had similar problem, or have another explanation for the source of problem.

The first bug I found was in ShellMITC4.cpp (line 860). the code shall be changed from :
------------------------
int count = 0;
for (i=0; i<4; i++) {
const Vector &Raccel = nodePointers[i]->getRV(accel);
for (int j=0; j<6; j++)
resid(count++) = Raccel(i);//this index i is not correct
}
------------------

TO:

------------------------
int count = 0;
for (i=0; i<4; i++) {
const Vector &Raccel = nodePointers[i]->getRV(accel);
for (int j=0; j<6; j++)
resid(count++) = Raccel(j);//notice i changed to j
}
------------------

Since i is the node index while j is the DOF index. it makes sense to access Raccel(j) or jth component of acceleration, as we already have the acceleration vector from node i. I am "almost sure" that this is a bug.


Another problem I noticed is in the same file (ShellMITC4.cpp), in the formInertiaTerms method. The method tries to add momentum values to the resid vector(which was populated with accelerations previously). I ended up commenting these two lines of code (981 and 982):
//for ( p = 0; p < 3; p++ )
//resid( jj+p ) += ( temp * momentum(p) ) ;

I am not sure whether that is a bug or if the code computes momentum incorrectly. But it seemed to me that it was not necessary to add momentum to resid vector or at least it was causing my linear dynamic analysis to diverge. It is possible that this is important for another type of analysis.

[fmk]

can you send me a script .. yes to 1, -MA is needed in resid in 2, just not sure he is doing it right,

[elhaddad]

Thanks Dr. McKenna for the reply.

I posted a script previously here:
http://opensees.berkeley.edu/community/ ... =2&t=61397

It is a modification of one of the shell dynamics examples but I added ground excitation to it.
It assumes El Centro earthquake is in the ElCentroEQ.thf file (which i did not post here), but any other record can be used to regenerate the problem.

The script again is quoted here:

#--------------------------------------------------------------------------------------
# ----------------------------
# Start of model generation
# ----------------------------
model basic -ndm 3 -ndf 6

# create the material
section ElasticMembranePlateSection 1 3.0e3 0.25 1.175 1.27

# set some parameters for node and element generation
set Plate ShellMITC4

set eleArgs "1"

#these should both be even
set nx 8
set ny 2

#loaded nodes
set mid [expr ( ($nx+1)*($ny+1)+1 ) / 2 ]
set side1 [expr ($nx + 2)/2 ]
set side2 [expr ($nx+1)*($ny+1) - $side1 + 1 ]

# generate the nodes and elements
block2D $nx $ny 1 1 $Plate $eleArgs {
1 -20 0 0
2 -20 0 40
3 20 0 40
4 20 0 0
5 -10 10 20
7 10 10 20
9 0 10 20
}

# TimeSeries "ElCentroEQ": tsTag dt filePath cFactor
timeSeries Path 1 -dt +2.000000E-002 -filePath ElCentroEQ.thf -factor +2.000000E+000

#EQ excitation
pattern UniformExcitation 1 2 -accel 1


# define the boundary conditions
# rotation free about x-axis (remember right-hand-rule)
fixZ 0.0 1 1 1 0 1 1
fixZ 40.0 1 1 1 0 1 1


# ----------------------------
# Start of recorder generation
# ----------------------------
recorder Node -file Node.out -time -node $mid -dof 2 disp



# ---------------------------------------
# Create and Perform the dynamic analysis
# ---------------------------------------
# Create the transient analysis
test EnergyIncr 1.0e-10 20 0
algorithm Linear
numberer Plain
constraints Plain
system SparseSPD
#integrator GeneralizedMidpoint 0.50
integrator Newmark 0.50 0.25
analysis Transient

# Perform the transient analysis (20 sec)
analyze 1000 0.02

wipe
#----------------------------------

[elhaddad]

Hello Dr. McKenna,

I think I found out how to correct this bug. It looks like it is better to avoid using resid vector for this method and use a local variable.
I created a vector called allAccel and used it to store the accelerations.
Also I have uncommented the two lines (981 and 982), that I had commented before, so that the code computes the inertia terms correctly.

I have exported the mass and stiffness matrices into text files, and use them to verify the analysis results using Matlab.
My final corrected code looks like this:
//-----------------------
int
ShellMITC4::addInertiaLoadToUnbalance(const Vector &accel)
{
int tangFlag = 1 ;

int i;

int allRhoZero = 0;
for (i=0; i<4; i++) {
if (materialPointers[i]->getRho() != 0.0)
allRhoZero = 1;
}

if (allRhoZero == 0)
return 0;

int count = 0;
Vector allAccel(24) ;//new vector
for (i=0; i<4; i++) {
const Vector &Raccel = nodePointers[i]->getRV(accel);
for (int j=0; j<6; j++)
allAccel(count++) = Raccel(j);
//resid(count++) = Raccel(j);//this line was replaced
}

formInertiaTerms( tangFlag ) ;
if (load == 0)
load = new Vector(24);
load->addMatrixVector(1.0, mass, allAccel, -1.0);
//load->addMatrixVector(1.0, mass, resid, -1.0);//this line was replaced
return 0;
}
//-----------------------


I noticed the same bug also exist in the nonlinear shell element (shellNL.cpp). I haven't used that element, but the code looks the same.

[fmk]

i noticed the other as well. guess it was not tested with ele mass. i have sent an email to the codes author. a quick solution is also just to let the base class method handle it.

[erahmani]

Dear Frank,

I have plan to use the Shell element (ShellMITC4) with J2 Plasticity (or Drucker Prager) Material to model nonlinear wood and gypsum shear walls. I wonder if the potential bugs in this element have been resolved. If not, do you have any suggestion to use a different shell element?

[MIRomero]

Elhaddad,
I am having the same troubles that you faced in March, were you able to somehow include the mass in your shell element?
Thank you for the response.
I am trying to run properly a linear dynamic analysis in a rectangular tunnel with Shellelements.
I need my tunnel model to work to include it in a 3D soil model afterwards.
Thank you

[elhaddad]

Hi MIRomero,

Yes, I got this fixed in the source code I have and compiled it to an exe
Do you have a copy of the source that you can compile? or are you using the executable files provided by the website?
If you can modify the source and compile it, then you will need to do three modifications (like shown above), as follows:

1- add the line: Vector allAccel(24) ;//new vector

2- replace the line: resid(count++) = Raccel(i);
with the line:
allAccel(count++) = Raccel(j);

3- replace the line: load->addMatrixVector(1.0, mass, resid, -1.0);
with the line: load->addMatrixVector(1.0, mass, allAccel, -1.0);

I have verified the results using matlab and it works correctly

If you do not have the source to compile, let me know and I can share the file I compiled with you

Good Luck
Wael Elhaddad

[MIRomero]

Hi Wael Elhaddad,

Thank you for your answer.
I do not have the source to compile. I will really appreciate if you can share the file you compiled with me.
I have created a box folder, where you can upload it.
Here the link:
https://uofi.box.com/s/bzu8pfi95yykhd6jhn6z

Thank you in advance for your time.

Maria Ines Romero

[mmusgrove]

Hi Frank,

I was wondering if you had heard back from the code authors yet regarding the changes proposed by elhaddad and/or your suggestion of allowing the base class to handle it? I ask because I have successfully recompiled OpenSees with the changes proposed by elhadded, but have had no luck trying to get OpenSeesSP to compile (which is what MIRomero and I would prefer to use given the size/complexity of our models).

Thanks,
Michael Musgrove

[suiwenwu]

Hello, Michael, do you solve the problem with ShellMITC4 in OpenSees?

[mmusgrove]

Hi suiwenwu,

I never got OpenSeesSP to compile correctly, so we ended up switching to using BBarbrick instead of the shell elements. I saw in the SVN logs that changes had been made to the shell element code, but I haven't tried recompiling.

Michael

[fmk]

the Shell elements in the current version, 2.4.5 (rel 5855) work for transient analysis with uniform excitation .. the first modification suggested by elhadded was needed, the second modification suggested was not correct, an alternative modification was needed .. the revised code has been checked in and is available for downloading.

[mmusgrove]

Thanks for the update, Frank. I'll give it a try.

Michael