dr-bill-c / MYSTRAN

MYSTRAN is a general purpose finite element analysis solver
MIT License
60 stars 28 forks source link

CBUSH non-null length element #7

Closed numenic closed 3 years ago

numenic commented 3 years ago

I was comparing some results (CBAR and CBUSH only) against NX-NASTRAN and found issues with non-null length CBUSH. Checking `EMG/EMG3/BUSH.f90 (is it the good file?) we can indeed see that CBUSH stiffness matrix is not appropriate for length > 0.

I can provide test-cases and matrix formulation if required.

Thanks

dr-bill-c commented 3 years ago

The KE matrix in BUSH.f90 is in element coordinates so you won't see the effect of the 2 grids not being coincident. It gets transformed to basic coords and then to global coords. That transformation process is programmed to take into account the actual geometry of the grids to which the BUSH element is attached. If you could send me an example of the problem you are having (along with MSC results) I would be happy to look into it

numenic commented 3 years ago

Thanks for the answer (and thanks a lot for this promising software). Here is a simple test-case I ran with NX (unfortunately, I will not be able to provide MSC results). It consists of a single CBUSH grounded at node1, free at node 2, and loaded with a single nodal transverse load:

    grid2
     o<=== P=1000  
     |
     /
     \
     /
     \             ^ Z
     |             |  
    /// grid 1     ---> X
        (grounded)

test-case-cbush.zip

You can check that reactions are not consistent: No moment are recovered from SPC forces

MYSTRAN:

                                                          S P C   F O R C E S
                                              (in global coordinate system at each grid)
           GRID     COORD      T1            T2            T3            R1            R2            R3
                     SYS
              1        0  0.0          -1.000000E+03  0.0           0.0           0.0           0.0   

NX-NASTRAN

                               F O R C E S   O F   S I N G L E - P O I N T   C O N S T R A I N T

      POINT ID.   TYPE          T1             T2             T3             R1             R2             R3
             1      G      0.0           -1.000000E+03   0.0            1.000000E+04   0.0            0.0
1    ANALYSE                                                               FEBRUARY   1, 2021  NX NASTRAN  5/23/18   PAGE    13

This can also be seen from Forces in CBUSH, where only the applied transverse load is reacted in MYSTRAN CBUSH, whereas we can expect a moment to be reacted as well:

MYSTRAN

                          E L E M E N T   E N G I N E E R I N G   F O R C E S
                              F O R   E L E M E N T   T Y P E   B U S H    
  Element      Force         Force         Force        Moment        Moment        Moment
     ID         XE            YE            ZE            XE            YE            ZE
        1 -0.000000E+00 -0.000000E+00  1.000000E+03 -0.000000E+00 -0.000000E+00 -0.000000E+00

NX-NASTRAN:

                                 F O R C E S   I N   B U S H   E L E M E N T S        ( C B U S H )

                  ELEMENT-ID        FORCE-X       FORCE-Y       FORCE-Z      MOMENT-X      MOMENT-Y      MOMENT-Z  
0                          1      0.0           0.0           1.000000E+03  0.0          -5.000000E+03  0.0
1    ANALYSE                                                               FEBRUARY   1, 2021  NX NASTRAN  5/23/18   PAGE    14
dr-bill-c commented 3 years ago

Thanks. I will take a look at it. Can you send me the Bulk Data File?

numenic commented 3 years ago

Thanks to take it into consideration!

The bulk file is included in the archive posted in my second message "test-case-cbush.zip". You'll find inside the bulk file "nx_cbush.dat" and the two f06 files for comparison: "mystran_cbush".F06 and "nx_cbush.f06"

numenic commented 3 years ago

Hello @dr-bill-c

here is the 12x12 element stiffness matrix (in element CSYS) that shows excellent correlations with NASTRAN CBUSHS (it's python, but quite easy to adress, I guess).

self.prop.k-AXIS- equals NASTRAN K-AXIS-, whereas L is element's length.

        # ---------------------------------------------------------------------
        # diagonal
        L = self.length
        k1_1 = self.prop.k1
        k2_2 = self.prop.k2
        k3_3 = self.prop.k3
        k4_4 = self.prop.k4
        k5_5 = self.prop.k5 + (self.prop.k3 * L ** 2) / 4
        k6_6 = self.prop.k6 + (self.prop.k2 * L ** 2) / 4
        # -----------------------
        k7_7 = self.prop.k1
        k8_8 = self.prop.k2
        k9_9 = self.prop.k3
        k10_10 = self.prop.k4
        k11_11 = self.prop.k5 + (self.prop.k3 * L ** 2) / 4
        k12_12 = self.prop.k6 + (self.prop.k2 * L ** 2) / 4
        # ---------------------------------------------------------------------
        # small diagonal
        k1_7 = -k1_1
        k2_8 = -k2_2
        k3_9 = -k3_3
        k4_10 = -k4_4
        k5_11 = -self.prop.k5 + (self.prop.k3 * L ** 2) / 4
        k6_12 = -self.prop.k6 + (self.prop.k2 * L ** 2) / 4
        # -------------------------------------------------------------------------
        # remaining
        k2_6 = L * self.prop.k2 / 2
        k2_12 = L * self.prop.k2 / 2
        k3_5 = -(self.prop.k3 * L) / 2
        k3_11 = -(self.prop.k3 * L) / 2
        k5_9 = (self.prop.k3 * L) / 2
        k6_8 = -L * self.prop.k2 / 2
        k8_12 = -L * self.prop.k2 / 2
        k9_11 = (self.prop.k3 * L) / 2
        # ---------------------------------------------------------------------
        # Build the local matrix ke
        # fmt: off
        ke = np.array(
        # =============================================================================
       [[k1_1, 0,    0,    0,    0,    0,    k1_7, 0,    0,    0,      0,      0     ],
        [0,    k2_2, 0,    0,    0,    k2_6, 0,    k2_8, 0,    0,      0,      k2_12 ],
        [0,    0,    k3_3, 0,    k3_5, 0,    0,    0,    k3_9, 0,      k3_11,  0     ],
        [0,    0,    0,    k4_4, 0,    0,    0,    0,    0,    k4_10,  0,      0     ],
        [0,    0,    0,    0,    k5_5, 0,    0,    0,    k5_9, 0,      k5_11,  0     ],
        [0,    0,    0,    0,    0,    k6_6, 0,    k6_8, 0,    0,      0,      k6_12 ],
        [0,    0,    0,    0,    0,    0,    k7_7, 0,    0,    0,      0,      0     ],
        [0,    0,    0,    0,    0,    0,    0,    k8_8, 0,    0,      0,      k8_12 ],
        [0,    0,    0,    0,    0,    0,    0,    0,    k9_9, 0,      k9_11,  0     ],
        [0,    0,    0,    0,    0,    0,    0,    0,    0,    k10_10, 0,      0     ],
        [0,    0,    0,    0,    0,    0,    0,    0,    0,    0,      k11_11, 0     ],
        [0,    0,    0,    0,    0,    0,    0,    0,    0,    0,      0,      k12_12]])
        # =============================================================================
        # fmt: on
        ke = symmetrize(ke)
        return ke
dr-bill-c commented 3 years ago

I have programmed the above BUSH stiffness matrix you gave into MYSTRAN. Using the example input file (nx_bush.dat) MYSTRAN now gets the same answers as were in your NX NASTRAN run. Thanks for providing the equations for the stiffness matrix. I may try to add the effect of the BUSH element not being exactly at the center of the distance between the 2 grids. BUSH-F06.txt

dr-bill-c commented 3 years ago

Sorry. I didn't mean to close the issue. That should be your prerogative

numenic commented 3 years ago

Thanks Doc. It looks perfect now!