Incorrect bending about y-axis

[JohnDN90]

It seems that the bending about the y-axis is incorrect. It appears to be using the using the I for bending about the z-axis instead of the I for bending about the y-axis.

I believe the problem is is located here and here, and the indices need to be changed from (0,1) to (1,1) and from (0,5) to (1,5) respectively. I tested this change and it seemed to work on my end.

If you agree this is the appropriate fix, I can create a pull request for the fix.

[manavbhatia]

(0,1) and (0,5) are the correct locations for the shape functions for the Bmat_v operator. The intent of this operator is to provide a matrix, such that multiplication with the displacement vector gives the strain. The strain is

epsilon_xx =  du/dx - y d2v/dx2 - z d2w/dx2

Since both v and w displacements contribute to the epsilon_xx term, they both should be in the first row of the Bmat matrix.

The second row in this matrix provides the torsional component of the stiffness matrix.

I think the issue is in the definition of the section property matrix for bending here. The area inertia terms are defined in a 2x2 matrix with the diagonal elements set as the I_zz and I_yy terms in its (0,0) and (1,1) locations, respectively.

Then, this matrix is multiplied by the modulus of elasticity here, which is requested by the element here and then used to compute the element stiffness terms here and here . So, the problem is that we only ever use the (0,0) term from this matrix, which is I_zz needed for bending with v. Even bending with w will also use I_zz, which is likely what you are seeing.

One fix is to modify this to use the (0,0) term in the material_D_mat and this to use the (1,1) term in the material_D_mat.

One possible issue with this is that the off-diagonal terms I_yz will not get used in this modification, but that we can make explicit modifications in the internal_residual computations to account for that.

[JohnDN90]

One fix is to modify this to use the (0,0) term in the material_D_mat and this to use the (1,1) term in the material_D_mat. \ One possible issue with this is that the off-diagonal terms I_yz will not get used in this modification, but that we can make explicit modifications in the internal_residual computations to account for that.

I haven't tried this yet, but looking through the code. Couldn't we just leave this as is and then modify this to use the transpose of material_D_mat? I think that way they will both use the appropriate I term, and since material_D_mat is symmetric, the correct off-diagonal terms I_yz, should still be used.

I'll give that a shot and report back.

[manavbhatia]

This would require using (1, 0) and (1,5) as the location of entries for the Bmat_w. Is that what you had in mind?

[JohnDN90]

Hmm, no it wasn't. I'll take some time to look through the structural 1D code. I'll have to get more familiar with how things are set up since I'll have to modify some stuff to get accurate stress/strains in non-rectangular cross-section beams anyway.

[JohnDN90]

I haven't tried this yet, but looking through the code. Couldn't we just leave this as is and then modify this to use the transpose of material_D_mat? I think that way they will both use the appropriate I term, and since material_D_mat is symmetric, the correct off-diagonal terms I_yz, should still be used. \ I'll give that a shot and report back.

I read through the comments in this issue again, but still don't quite understand your last comment as to why using the transpose of material_D_mat isn't appropriate.

I have a simple one element cantilever beam model and a more complex model of an aircraft wing represented using only beams (I removed the plate elements for testing purposes). Using the material_D_mat.transpose() as I suggested above (without any modifications to the B_mat_v or B_mat_w matrices), both models match Nastran results.

Could you explain why using material_D_mat.transpose() isn't the correct method?

You said the bottom row of B_mat is for torsion. So, for the strains, we should have: epsZ = B_mat_w * u_e = [epsilon_bendZ; epsilon_torsion] epsY = B_mat_v * u_e = [epsilon_bendY; epsilon_torsion]

then multiplying by material_D_mat should yield

sigmaZ = material_D_mat * epsZ = [E*Izz*epsilon_bendZ + off_diag*epsilon_torsion; off_diag*epsilon_bendZ + E*Iyy*epsilon_torsion] sigmaY = material_D_mat.transpose() * epsY = [E*Iyy*epsilon_bendY + off_diag*epsilon_torsion; off_diag*epsilon_bendY + E*Izz*epsilon_torsion]

Since the residual goes to zero and the result matches the results from Nastran, I'm having trouble understanding why this method isn't appropriate.

[JohnDN90]

Also, could you clarify what the y- and z-directions are when defining a beam in MAST? If looking at the beam's cross section, is the y-direction up/down the screen and z-direction left/right across the screen (like Nastran defines it) or is it vice versa?

[manavbhatia]

material_D_mat = E *[ Izz  Izy;  Iyz  Iyy]

and its transpose is

material_D_mat.transpose() = E *[ Izz  Iyz;  Izy  Iyy] 

So, the Izz and Iyy diagonal terms do not get moved by the transpose operations. If you use (0, 0) and (0, 5) as the entries of the shape functions in the B_mat_w matrix, then the direct strain will be multiplied by Izz irrespective of transpose.

However, if you modify this to (1, 0) and (1,5), then this will be fine. This is what I was noting in my last message.

[JohnDN90]

My apologies, I feel a little ridiculous for making that simple mistake. That make complete sense now. So, the reason why it is still working and matching Nastran is that I'm not actually taking the transpose. I'm manually re-arranging material_D_mat, by: material_D_mat_2 = RealMatrixX::Zero(2,2); material_D_mat_2(0,0) = material_D_mat(1,1) material_D_mat_2(1,1) = material_D_mat(0,0) material_D_mat_2(0,1) = material_D_mat(1,0) material_D_mat_2(1,0) = material_D_mat(0,1)

Then it works out. Does that make more sense now?

[manavbhatia]

epsZ = B_mat_w u_e = [epsilon_bendZ; epsilon_torsion] epsY = B_mat_v u_e = [epsilon_bendY; epsilon_torsion]

Note that B_mat_v and B_mat_w don't have any shape functions for torsion. So, the correct definition should be

 epsZ = B_mat_w * u_e = [epsilon_bendZ; 0]
 epsY = B_mat_v * u_e = [epsilon_bendY; 0]

Then, if you multiple this with the material_D matrix this becomes

material_D * epsZ = [E Izz epsilon_bendZ; 0]
material_D * epsY = [E Izz epsilon_bendY; 0]

So, with a (0,0) and (0,5) location of shape function in Bmat_z this does not use Iyy.

Also, for the sake of accuracy, we should clarify that these are the section integrated moments stress, which gives us the inertia terms Izz and Iyy. For stress these should be y and z only.

[manavbhatia]

My apologies, I feel a little ridiculous for making that simple mistake. That make complete sense now. So, the reason why it is still working and matching Nastran is that I'm not actually taking the transpose. I'm manually re-arranging material_D_mat, by: material_D_mat_2 = RealMatrixX::Zero(2,2); material_D_mat_2(0,0) = material_D_mat(1,1) material_D_mat_2(1,1) = material_D_mat(0,0) material_D_mat_2(0,1) = material_D_mat(1,0) material_D_mat_2(1,0) = material_D_mat(0,1)

Then it works out. Does that make more sense now?

Yes, this does seem ok.

[manavbhatia]

Also, could you clarify what the y- and z-directions are when defining a beam in MAST? If looking at the beam's cross section, is the y-direction up/down the screen and z-direction left/right across the screen (like Nastran defines it) or is it vice versa?

x- is along the element axis. y- is created using the y-vector provided by the use in the section card, which defines the x-y plane. Then, z- is obtained from the cross-product.

[JohnDN90]

This should be fixed in the upcoming catch2 unit tests capability pull request. Once we confirm its been merged into master we can then close this issue.

[jdeaton]

Closed by PR #75