1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112
| double Iy = M_PI*(pow(p_pipe->OD,4)-pow(p_pipe->OD-2*p_pipe->WT,4))/64; //Moment quadratique
double Iz = M_PI*(pow(p_pipe->OD,4)-pow(p_pipe->OD-2*p_pipe->WT,4))/64; //Moment quadratique
double J = M_PI*(pow(p_pipe->OD,4)-pow(p_pipe->OD-2*p_pipe->WT,4))/32; //Moment quadratique polaire
double S = M_PI*(pow(p_pipe->OD,2)-pow(p_pipe->OD-2*p_pipe->WT,2))/4; //Section du pipe
double G = p_pipe->E/(2*(1+p_pipe->poissonRatio)); //Module de cisaillement
double E = p_pipe->E; //Module d'Young
meshPipe(); //Remplissage de "nodes" (création des noeuds)
//###################################### CONSTRUCTION DE LA MATRICE DE RAIDEUR ###################################################
mat.setSize(6*nodes.size(),6*nodes.size());
Matrix *matElement = new Matrix; //On définit une matrice d'un élément (12x12) (adaptée pour chaque élément)
matElement->setSize(12,12);
double L = 0; //Longueur de l'élément (adaptée pour chaque élément)
for(int k=0;k<nodes.size()-1;k++)//On itère sur les éléments (nombre de noeuds-1)
{
L = nodes.at(k).data.distanceToPoint(nodes.at(k+1).data); //Longueur de l'élément
matElement->resetToZero(); //Remplit matElement de 0
//Remplissage du triangle supérieur de la matrice
(*matElement)(0,0) = E*S/L;
(*matElement)(0,6) = -E*S/L;
(*matElement)(1,1) = 12*E*Iy/pow(L,3);
(*matElement)(1,5) = 6*E*Iz/pow(L,2);
(*matElement)(1,7) = -12*E*Iz/pow(L,3);
(*matElement)(1,11) = 6*E*Iz/pow(L,2);
(*matElement)(2,2) = 12*E*Iy/pow(L,3);
(*matElement)(2,4) = -6*E*Iy/pow(L,2);
(*matElement)(2,8) = -12*E*Iy/pow(L,3);
(*matElement)(2,10) = -6*E*Iy/pow(L,2);
(*matElement)(3,3) = G*J/L;
(*matElement)(3,9) = -G*J/L;
(*matElement)(4,4) = 4*E*Iy/L;
(*matElement)(4,8) = 6*E*Iy/pow(L,2);
(*matElement)(4,10) = 2*E*Iy/L;
(*matElement)(5,5) = 4*E*Iz/L;
(*matElement)(5,7) = -6*E*Iz/pow(L,2);
(*matElement)(5,11) = 2*E*Iz/L;
(*matElement)(6,6) = E*S/L;
(*matElement)(7,7) = 12*E*Iz/pow(L,3);
(*matElement)(7,11) = -6*E*Iz/pow(L,2);
(*matElement)(8,8) = 12*E*Iy/pow(L,3);
(*matElement)(8,10) = 6*E*Iy/pow(L,2);
(*matElement)(9,9) = G*J/L;
(*matElement)(10,10) = 4*E*Iy/L;
(*matElement)(11,11) = 4*E*Iz/L;
//Remplissage de la moitié inférieure par symétrie
for(int i=0;i<matElement->rowCount();i++)
{
for(int j=0;j<i;j++)
{
(*matElement)(i,j) = (*matElement)(j,i);
}
}
//On crée une matrice rotation pour mettre la matrice de rigidité de l'élément dans le repère global.
//On utilisera [K]=[R]t[K'][R] pour changer de repère (avec [R]t = [R]-1, transposée egal à l'inverse car transformation
//orthogonale).
Matrix rotation;
rotation.setSize(12,12);
rotation(0,0) = (nodes.at(k+1).data.x-nodes.at(k).data.x)/nodes.at(k).data.distanceToPoint(nodes.at(k+1).data);
rotation(0,1) = (nodes.at(k+1).data.y-nodes.at(k).data.y)/nodes.at(k).data.distanceToPoint(nodes.at(k+1).data);
rotation(0,2) = (nodes.at(k+1).data.z-nodes.at(k).data.z)/nodes.at(k).data.distanceToPoint(nodes.at(k+1).data);
rotation(1,0) = -rotation(0,1)/(pow(pow(rotation(0,1),2)+pow(rotation(0,0),2),0.5));
rotation(1,1) = rotation(0,0)/(pow(pow(rotation(0,1),2)+pow(rotation(0,0),2),0.5));
rotation(1,2) = 0;
rotation(2,0) = -rotation(0,2)*rotation(1,1);
rotation(2,1) = -rotation(0,2)*rotation(1,0);
rotation(2,2) = rotation(0,0)*rotation(1,1)-rotation(1,0)*rotation(0,1);
//On remplit le reste de la matrice pour obtenir une matrice de rotation 12x12 (taille de la matrice rigidité)
for(int l=0;l<4;l++)
{
for(int j=0;j<3;j++)
{
for(int i=0 ; i<3 ; i++)
{
rotation(l*3+i,l*3+j) = rotation(i,j);
}
}
}
Matrix rotationTranspose;
rotationTranspose.setSize(12,12);
for(int i=0;i<12;i++)
{
for(int j=0;j<12;j++)
{
rotationTranspose(i,j) = rotation(i,j);
}
}
rotationTranspose.transpose();
//Changement de repère de la matrice grâce à la matrice de rotation
(*matElement) = rotationTranspose*(*matElement)*rotation;
//Assemblage avec la matrice globale de rigidité
for(int i=0;i<matElement->rowCount();i++)
{
for(int j=0;j<matElement->columnCount();j++)
{
mat(k*6+i,k*6+j) += (*matElement)(i,j);
}
}
} |
Partager