Skip to content

Commit 5b61795

Browse files
authored
fixed the indefinite mass matrix for beam element. Thanks to John Neiferd for identifying the issue. (#54)
1 parent 6eb2e6a commit 5b61795

File tree

1 file changed

+25
-20
lines changed

1 file changed

+25
-20
lines changed

src/property_cards/solid_1d_section_element_property_card.cpp

Lines changed: 25 additions & 20 deletions
Original file line numberDiff line numberDiff line change
@@ -928,14 +928,16 @@ InertiaMatrix::operator() (const libMesh::Point& p,
928928
m(3,3) = Ip;
929929

930930
// rotational velocities
931-
m(0,4) = Ay; m(4,0) = Ay; // w-displacement
932-
m(0,5) = -Az; m(5,0) = -Az; // v-displacement
933-
934-
// bending rotation inertia
935-
for (unsigned int i=0; i<2; i++)
936-
for (unsigned int j=0; j<2; j++)
937-
m(4+i,4+j) = I(i,j)*1.e-6;
938-
931+
// theta-y rotation
932+
m(0,4) = Ay; m(4,0) = Ay;
933+
m(4,4) = I(1,1); // I11 is defined about y-y axis for theta-y
934+
935+
// theta-z rotation
936+
m(0,5) = Az; m(5,0) = Az;
937+
m(5,5) = I(0,0); // I00 is defined about z-z axis for theta-z
938+
939+
m(4,5) = m(5,4) = I(0,1);
940+
939941
m *= rho;
940942
}
941943

@@ -972,18 +974,21 @@ InertiaMatrix::derivative ( const MAST::FunctionBase& f,
972974
dm(3,3) = dIp;
973975

974976
// rotational velocities
975-
m(0,4) = Ay; m(4,0) = Ay; // w-displacement
976-
dm(0,4) = dAy; dm(4,0) = dAy; // w-displacement
977-
m(0,5) = -Az; m(5,0) = -Az; // v-displacement
978-
dm(0,5) = -dAz; m(5,0) = -dAz; // v-displacement
979-
980-
// bending rotation inertia
981-
for (unsigned int i=0; i<2; i++)
982-
for (unsigned int j=0; j<2; j++) {
983-
m(4+i,4+j) = I(i,j)*1.e-6;
984-
dm(4+i,4+j) = dI(i,j)*1.e-6;
985-
}
986-
977+
// theta-y rotation
978+
m(0,4) = Ay; m(4,0) = Ay;
979+
m(4,4) = I(1,1);
980+
dm(0,4) = dAy; dm(4,0) = dAy;
981+
dm(4,4) = dI(1,1);
982+
983+
// theta-z rotation
984+
m(0,5) = Az; m(5,0) = Az;
985+
m(5,5) = I(0,0);
986+
dm(0,5) = dAz; dm(5,0) = dAz; // v-displacement
987+
dm(5,5) = dI(0,0);
988+
989+
m(4,5) = m(5,4) = I(0,1);
990+
dm(4,5) = dm(5,4) = dI(0,1);
991+
987992
m *= drho;
988993
m += rho*dm;
989994
}

0 commit comments

Comments
 (0)