Modifying Elasticbeam2d ?

[kostast88]

Hello,

I'm trying to modify the elasticbeam2d.cpp in an attempt to account for PM Interaction without fibers. I've read and understood the cpp file, but some parts remain unclear to me. From what I've seen, the resistingforce() calculates a vector q{0,1,2} which is actually {N,M1,M2} and together with the p0{1,2,3} which is {N,V1,V2}, they form the p0vec which is send to CrdTransf (assume lineardcrdtransf2d). In there we have vectors:
q: is actually the basic forces
pl: is local forces,which is actually basic forces with the addition of V=(M1+M2)/L
pg is calculated through coordinate transformation ->global forces

pg returns to getresistingforce() in the element cpp, as the new vector P,and substracting Q (mass forces), the global forces vector is final.
Example for tests is EX1.tcl from examples manual,with LoadControl.

After many attempts, I added the following arguments in the end of getresistingforce() to test if a PM interaction algorithm would work:

[code]
P = theCoordTransf->getGlobalResistingForce(q, p0Vec);
///////////////////CHANGES//////////////////////////
if (P(1)>300) {
P(2)=P(2)*0.65;
P(5)=P(5)*0.65;
}
///////////////END OF CHANGES/////////////////////
// P = P - Q;
P.addVector(1.0, Q, -1.0);
return P;
[end of code]

This would change the global forces vector (moments only) when the global vertical (local axial) force exceeds 300 units. When I use this example everything is fine.
This change is for the P vector, so global forces only are affected, while local forces remain unchanged.
When putting 2 elastic members in a row, things change and algorithm does not work (not in the desirable way,moments change are random).

2 Questions
1) How do 2 elements (elasticbeamcolumn) exchange forces and moments?And which are those (basic,local,global)? This V=(M1+M2)/L thing in CrdTransf won't allow me to change only moments, without changes in V1,V2 at the the beginning of the code.
2)Should I keep up to that, or tottaly change the elasticbeam2d element, especially in getresisting force, writing down all the matrix equations,not just the three of them that form q (i mean q{0,1,2,3,4,5} instead of q{0,1,2}) and do the transformations in the element itself?

Any hint or tip, is most welcomed! Thank you in advance.

[fmk]

1) the interaction between the elements occurs at the global finite element level

foreach element in elements {
K.addStifness(ele.stiffness)
add ResistingForces(ele.resistingForce)
}
the element stifness and force for the 2d beam case are 6x6 matrix and a 6x1 vector. the transformation
takes the basic 3x3 matix and 3x1 vector and transforms them

2) if the transformation is messing you up .. go to the system you are used to, probably the full 6x6 and 6x1 ... the formulation
for the basic elastic element small displ and small defo is after all pretty easy and typically given in any intro textbook on matrix
structural analysis.

[kostast88]

Thank you very much professor Mckenna for your prompt response!

[kostast88]

One more question. Writing down all the equations at getResistingforce() in elasticbeam2d would require the whole displacements vector.So in :

const Vector &v= theCoordTransf->getBasicTrialDisp();

I need v(3),v(4) and v(5) as well.Does this method (gettrialdisp) include them, or do i have to change this method as well?
Watching the method in linearcrdtransf2d I say yes.
From a test,I say no, gettrialdisp returns a 3x1 vector, and new equations result to: "crt detected that the application wrote to memory after end of heap" errors occur.

So I'm confused If they exist at all... (should,since they can be given as output).

P.S. The algorithm was consistent in EX1.tcl but gave extreme numerical instabilities in a simple frame(with many solution algorithms and default convergence tolerance). Probably because global forces only change,but basic and local remain unchanged.Yes this is weird but it's the only way around...
I'm starting to wonder whether this can be done at all!

[EDIT]:I think I should create a new method in crdtransf2d (LinearCrdTransf2d::NEWgetBasicTrialDisp (void)), which will calculate ub(3),ub(4),ub(5) and then call this function from the element,instead of the original one.Then write down the whole equations.Is this correct?

[fmk]

i thought you were going to ignore the transformations .. if so go to the nodes themselves to get the nodal displacements .. if you want to keep the transformations then your element must be written in the basic system, in which case you only want the transformation to return a 3x1.

[kostast88]

Yes, I figured that. It should be like :

const Vector &disp1 = theNodes[0]->getTrialDisp();
const Vector &disp2 = theNodes[1]->getTrialDisp();

static double ug[6];
for (int i = 0; i < 3; i++) {
ug[i] = disp1(i);
ug[i+3] = disp2(i);

so I get the ug[6] with global displacements.

But I can't understand why, with the previous formulation (moments*0.65), in transient analysis, the P vector (global forces) has changes in P(0),P(4) (axial forces).
I didn't even touch them! Does this have to do with the algorithm(Newton) that has to satisfy F=k*u and since I changed F, this doesn't converge? Because tracing deep in these methods, didn't gave me any answer.Convergence was an issue in my tests.

[fmk]

Newton is going to kep changing things until convergence is acheived, if you play with the moments then yes it is going to keep iterating ..

make sure your element end forces are in equi (so if you do change the moments at the ends, do change the shear forces).