Modal analysis using ChElementBeamIGA

[Nicolas CORTES]

Hello,

I’m trying to study the dynamic of a beam with initial twist and curvatures (prebend) using ChElementBeamIGA.

First, I created a homogeneous beam with a circular section without neither twist nor curvatures, clamped at one end and free at the other. In this case, everything is perfectly working, the deformed shape in response of an applied static load looks good. Furthermore, the modal analysis realised using the matrices system.GetStiffnessMatrix(K) and system.GetMassMatrix(M) gives the good modal frequencies (w2, vr = sp.linalg.eigh(K, M)).

However, when I introduce linear or non-linear twist to this beam, the modal analysis doesn’t work anymore (w2 becomes negative).

To introduce the twist, I had given a quaternion to the frame of my nodes :
###
twist_angle = (self.twists[n_section] + self.twists[n_section+p])/2

t = (Coord_nodeB – Coord_nodeA)
t.Normalize()
tan_twist = np.tan(self.twists[i]/2)
m_rot_twist.Set_A_Rodriguez(t*tan_twist)
q_rot_twist = m_rot_twist.Get_A_quaternion()

# Assemblage
q_rot = q_rot_twist * q_rot_prebend
m_rot = chrono.ChMatrix33D()
m_rot.Set_A_quaternion(q_rot)

node = fea.ChNodeFEAxyzrot(chrono.ChFrameD(pos, m_rot))
###

I also added SectionRotation and ShearRotation to ChElasticityCosseratAdvancedGenericFPM :
###
melasticity = fea.ChElasticityCosseratAdvancedGenericFPM(mKlaw, 0, 0, 0, 0, 0, 0) melasticity.SetSectionRotation(twist_angle)
melasticity.SetShearRotation(twist_angle)
###

As my beam has a circular section and has the same properties in y and z, the twisted beam should be exactly the same that the untwisted beam. However, it appears that the matrices are changed in such a way that the modal analysis is not possible anymore.

I looked in the code and found two places where the twist affects the matrices. The first one is the rotation of the material stiffness matrix in :
###
void ChElasticityCosseratAdvancedGenericFPM::UpdateEMatrix() {
        if (!updated) {  // do it only once
                  // compute T
                   ComputeTransformMatrix();
                    // update Klaw after setting section rotation and EC/SC offset
                    this->Klaw = this->T.transpose() * this->Klaw * this->T;
                    updated = true;
                    }
            }
###

This matrix Klaw seemed to be rotated around the wrong axis for the shear part and some coupling coefficient between shearZ and Torsion was appearing. So, I changed the line 305 from “RotShearA.Set_A_Rxyz(ChVector<>(this->beta, 0, 0));” to “RotShearA.Set_A_Rxyz(ChVector<>(0, 0, this->beta));” but it did not resolve my problem of modal analysis.

The second place where the twist affects the matrices is in ChElementBeamIGA::ComputeInternalForces_impl with the use of the initial strains this->strain_e_0[ig] and this->strain_k_0[ig] but I did not look farther.

Whereas, when the initial twist is constant over the beam, the modal analysis is working. As a results, it seems that the problem is appearing when an element is confronted to variable twist.

Am I doing something wrong while giving my initial twist angle ?
Could you help me to solve this problem please ?

 

Best regards,

Nicolas