Reaction force - Equilibrium problem

[jaganmail]

Describe the bug

Code: import numpy as np

Import 'FEModel3D' and 'Visualization' from 'PyNite'

from PyNite import FEModel3D

# Create a new model
frame = FEModel3D()

# Create members (all members will have the same properties in this example)
J = 0.1
Iy = 1000000
Iz = 100
A = 100

# Define a material
E = 99999999
G = 10000
nu = 0.01
rho = 100
frame.add_material('Steel', E, G, nu, rho)

# Read node coordinates from a text file
node_coordinates = np.loadtxt('MLG_node.txt', delimiter=',')

# Define the nodes
num_rows = len(node_coordinates)
for row in range(num_rows):
    frame.add_node('N'+str(node_coordinates[row,0].astype(int)),node_coordinates[row,1],node_coordinates[row,2],node_coordinates[row,3])
    print('N'+str(node_coordinates[row,0].astype(int)),node_coordinates[row,1],node_coordinates[row,2],node_coordinates[row,3])

# Read connectivity table from a text file
element_con = np.loadtxt('MLG_element.txt', dtype=int, delimiter=',')

# Define the elements
num_rows = len(element_con)
for row in range(num_rows):
    frame.add_member('M'+str(element_con[row,0]), 'N'+str(element_con[row,1]), 'N'+str(element_con[row,2]), 'Steel', Iy, Iz, J, A)
    print('M'+str(element_con[row,0]), 'N'+str(element_con[row,1]), 'N'+str(element_con[row,2]), 'Steel', Iy, Iz, J, A)

# Define the supports
frame.def_support('N1', True, True, True, False, False, False)
frame.def_support('N3', True, True, True, False, False, False)

# 20,0,1,0,1,1,1
# 21,1,0,0,1,1,1 
# beam.def_support('N1', True, True, True, False, False, False)
# beam.def_support('N2', True, True, True, True, False, False)

# Read connectivity table from a text file
force_con = np.loadtxt('MLG_force.txt', delimiter=',')
# Add nodal loads
# Define the elements
num_rows = len(force_con)
for row in range(num_rows):
    frame.add_node_load('N'+str(force_con[row,0].astype(int)),'FX',force_con[row,1])
    frame.add_node_load('N'+str(force_con[row,0].astype(int)),'FY',force_con[row,2])
    frame.add_node_load('N'+str(force_con[row,0].astype(int)),'FZ', force_con[row,3])
    frame.add_node_load('N'+str(force_con[row,0].astype(int)),'MX',force_con[row,4])
    frame.add_node_load('N'+str(force_con[row,0].astype(int)),'MY',force_con[row,5])
    frame.add_node_load('N'+str(force_con[row,0].astype(int)),'MZ', force_con[row,6])
    print('N'+str(force_con[row,0].astype(int)),'FX',force_con[row,1])

# Analyze the model
#frame.analyze(check_statics=True)
frame.analyze_PDelta(log=True, check_stability=True, max_iter=10000, tol=0.0001, sparse=True)

# Print reactions at each end of the beam
print('Left Support Reaction:', frame.Nodes['N1'].RxnFZ, 'kg')
print('Right Support Reacton:', frame.Nodes['N3'].RxnFZ, 'kg')
# Print reactions at each end of the beam
print('Left Support Reaction:', frame.Nodes['N1'].RxnFX, 'kg')
print('Right Support Reacton:', frame.Nodes['N3'].RxnFX, 'kg')
# Print reactions at each end of the beam
print('Left Support Reaction:', frame.Nodes['N1'].RxnFY, 'kg')
print('Right Support Reacton:', frame.Nodes['N3'].RxnFY, 'kg')

# Print reactions at each end of the beam
print('Left Support Reaction:', frame.Nodes['N1'].RxnMZ, 'kg')
print('Right Support Reacton:', frame.Nodes['N3'].RxnMZ, 'kg')
# Print reactions at each end of the beam
print('Left Support Reaction:', frame.Nodes['N1'].RxnMX, 'kg')
print('Right Support Reacton:', frame.Nodes['N3'].RxnMX, 'kg')
# Print reactions at each end of the beam
print('Left Support Reaction:', frame.Nodes['N1'].RxnMY, 'kg')
print('Right Support Reacton:', frame.Nodes['N3'].RxnMY, 'kg')

Inputs:

node; 1,0.751,0,-0.0904 2,1.3244,0.0,-0.0904 3,1.681,0,-0.0904 4,1.350,-1.10,-1.025 5,1.3244,-0.4592,-0.763 6,1.3244,0.0,-0.763 7,1.3244,0.4592,-0.763 8,1.350,1.10,-1.025

element: 1,1,2 2,2,3 3,2,6 4,4,5 5,5,6 6,6,7 7,7,8 8,6,1 9,6,3 10,4,6 11,8,6 12,2,5 13,2,7

force; 4,325.0,0.0,1050.0,0,0,0 8,0.0,0.0,0.0,0,0,0 5,-19.424114840907883,23.29565645762778,-76.14771922669809,0,0,0 7,-16.131767371978825,23.29565645762778,-38.328922473255254,0,0,0

I am unable to get the reaction properly. Please see the attached files. The P-delta reaction forces are wrong/ i am unable to get proper reaction forces from the analysis

[SoundsSerious]

You are running pretty small values for your off-axis moments of inertia, and rotational inertia

# Create members (all members will have the same properties in this example)
J = 0.1 <very small
Iy = 1000000
Iz = 100 < small as well
A = 100

Can you provide MLG_node.txt so I can give this a shot?

You're running the analyze_PDelta, which applies loading after deformation. For very small stiffness values this would create an unstable structure potentially diverging from a textbook answer. I would think that I / J values should be at least as large as A and potentially closer to Iy

Could you try this with analyze_linear and see if you get the right values? Alternatively, try increasing the stiffness and see if you get the correct reaction forces?

[jaganmail]

Hi JWock82, Thank you. Actually i am trying out all possibilities. My whole interest is to get the reaction force and moments estimated at the fixing points. I am not interested in stresses/deflection etc. Please find the enclosed Files for use. Particularly Fy and Ry are not matching at the end. Actually I changed J to 0.5. Now it seems to be working for MLG. I have the same problem with another frame model. I am confused about the complexity of this Simple reaction estimate. Jagan

On Tue, Aug 15, 2023 at 12:43 PM SoundsSerious @.***> wrote:

You are running pretty small values for your off-axis moments of inertia, and rotational inertia

Create members (all members will have the same properties in this example)

J = 0.1 <very small Iy = 1000000 Iz = 100 < small as well A = 100

Can you provide MLG_node.txt so I can give this a shot?

— Reply to this email directly, view it on GitHub https://github.com/JWock82/PyNite/issues/173#issuecomment-1678515777, or unsubscribe https://github.com/notifications/unsubscribe-auth/AD23BQFK5VK7FNKELVVJT53XVMOSTANCNFSM6AAAAAA3QUM53U . You are receiving this because you authored the thread.Message ID: @.***>

4,325.0,0.0,1050.0,0,0,0 8,0.0,0.0,0.0,0,0,0 5,-19.424114840907883,23.29565645762778,-76.14771922669809,0,0,0 7,-16.131767371978825,23.29565645762778,-38.328922473255254,0,0,0 1,0.751,0,-0.0904 2,1.3244,0.0,-0.0904 3,1.681,0,-0.0904 4,1.350,-1.10,-1.025 5,1.3244,-0.4592,-0.763 6,1.3244,0.0,-0.763 7,1.3244,0.4592,-0.763 8,1.350,1.10,-1.025 1,1,2 2,2,3 3,2,6 4,4,5 5,5,6 6,6,7 7,7,8 8,6,1 9,6,3 10,4,6 11,8,6 12,2,5 13,2,7 1,4.8435,-1.6700,0.1688 2,4.8435,-1.0450,0.1688 3,4.8435,-0.5000,0.1688 4,4.8435,0.0000,0.1688 5,4.8435,0.5000,0.1688 6,4.8435,1.0450,0.1688 7,4.8435,1.6700,0.1688 8,4.6290,0.0000,-0.0904 9,4.8730,0.0000,-0.0904 1,1,2 2,2,3 3,3,4 4,4,5 5,5,6 6,6,7 7,4,8 8,4,9 1,-1.04,-2.10,-12.40,0.00,0.00,0.00 2,-0.76,-1.84,-9.15,0.00,0.00,0.00 3,-0.87,-2.51,-10.49,0.00,0.00,0.00 4,-1.09,-3.89,-13.34,0.00,0.00,0.00 5,-0.54,-2.51,-6.75,0.00,0.00,0.00 6,-0.27,-1.84,-3.44,0.00,0.00,0.00 7,-0.14,-2.10,-1.98,0.00,0.00,0.00

[jaganmail]

Tried all your suggestions. No Sucess. Jagan

[JWock82]

Are you mixing up the strong and weak axes? Iy = weak axis. Iz = strong axis.

[jaganmail]

Hi JWock82, Thank you for your time. I am working on a project to create Shear force diagram, BMD for a component in an aircraft structure. Iy and Iz may probably be same in this case (Like consider a wing idealized as a beam element with element between each station). Given an external 3D load and moment at various node. my whole interest is to estimate the reaction and pass it on to the main component for SFD and BMD. Up to now (At least until now), My understanding is given a force and B.C, and only the proper geometry with any rigid assumption (E, I) a linear frame element must easily produce the reaction forces as we expected to calculate from the simple beam/frame analysis. I tried to play with Iy, Iz, J and A , also E and G to create a rigid structure. but no success. Is any assumption you suggest solving the above problem.

[JWock82]

Can you send me a simple script that reproduces the issue? I'll have a look at it. I use the reactions all the time with no issues. There is one known issue related to a stiffness approximation for plate drilling that changes the reactions - but it's usually insignificant.

[jaganmail]

Hi Craig, Thank you. I attached all the script and input in earlier mail too; it looks like GitHub mail id may not be able to deliver attachments. Again I am attaching the same for your reference and One docx sheet, where I tried various input changes to make a rigid frame.

Regards Jagan

On Tue, Aug 15, 2023 at 10:44 PM Craig @.***> wrote:

Can you send me a simple script that reproduces the issue? I'll have a look at it. I use the reactions all the time with no issues. There is one known issue related to a stiffness approximation for plate drilling that changes the reactions - but it's usually insignificant.

— Reply to this email directly, view it on GitHub https://github.com/JWock82/PyNite/issues/173#issuecomment-1679313440, or unsubscribe https://github.com/notifications/unsubscribe-auth/AD23BQEJJA47UZ3AAUIJPPLXVOU6FANCNFSM6AAAAAA3QUM53U . You are receiving this because you authored the thread.Message ID: @.***>

1,4.8435,-1.6700,0.1688 2,4.8435,-1.0450,0.1688 3,4.8435,-0.5000,0.1688 4,4.8435,0.0000,0.1688 5,4.8435,0.5000,0.1688 6,4.8435,1.0450,0.1688 7,4.8435,1.6700,0.1688 8,4.6290,0.0000,-0.0904 9,4.8730,0.0000,-0.0904 1,-1.04,-2.10,-12.40,0.00,0.00,0.00 2,-0.76,-1.84,-9.15,0.00,0.00,0.00 3,-0.87,-2.51,-10.49,0.00,0.00,0.00 4,-1.09,-3.89,-13.34,0.00,0.00,0.00 5,-0.54,-2.51,-6.75,0.00,0.00,0.00 6,-0.27,-1.84,-3.44,0.00,0.00,0.00 7,-0.14,-2.10,-1.98,0.00,0.00,0.00

1,1,2 2,2,3 3,3,4 4,4,5 5,5,6 6,6,7 7,4,8 8,4,9 9,8,1 10,8,2 11,8,3 12,8,4 13,8,5 14,8,6 15,8,7 16,9,1 17,9,2 18,9,3 19,9,4 20,9,5 21,9,6 22,9,7

[jaganmail]

HT_element.txt HT_force_inertia.txt HT_node_inertia.txt Program.txt Hi Craig, Please see the attached files. for your use Jagan

[jaganmail]

Program.txt Sorry Attached the python code too

[jaganmail]

Hi Craig, Sorry to bother you. Is there any success? Any mistakes from my side to be corrected? Thank you, Jagan.

[JWock82]

I haven't had a chance to dive into it just yet. It's been a busy week. I'll have a look this weekend.

[jaganmail]

Hi Craig, Thank you. Jagan

[SoundsSerious]

@jaganmail

Just catching up here. I am having trouble understanding something.

In your reply about modifying section properties you say: Actually I changed J to 0.5. Now it seems to be working for MLG. Now i'm not sure what MLG is, but if changing J 5x from 0.1 to 0.5 which is still a very small value compared to your stiffest dimension of 1M then there's a good case the structure as defined wasn't physically feasible.

I'm not sure what your background is with FEA but PyNite works based on the principle of Virtual Work. This essentially minimizes the stress x strain product in a structure. If a structure isn't stiff this method doesn't work well. This is what I suspect is happening in your case.

A couple of requests / questions so that we can debug your case:

1. If fixing the stiffness actually fixed your first case then is it possible it could fix your second case as well?
2. Please put all your code per case in one single, compilable python file so we don't have to merge your texts files.
3. Please provide the expected reaction forces you think are correct per case

[jaganmail]

Hi Craig, Thank you for your time. I had an introductory course on FEM and solid mechanics and have a basic idea of beam theory (I did not have any idea about the P-delta effect until now, I am reading to understand more). Sorry I will try to put my problems more clearly here. 1) I have different structures connected by beam elements. 2) I have the overall dimension-like length (Easy to measure using a ruler) of the structure, but I do not know the cross-sectional properties. 3) I also know the fixing point of this structure with another structure 4) Given the known length, I have idealized and created the structure as a connection of beam and created the model, and applied the BC and forces from theoretical calculations at some nodal points on the structure 5) I have created the code in such a way that, if you keep all the files in a folder and open the program file and run it, it will automatically fetch the inputs and provide the output, please try it, if difficulties, I will try to add, nodal, elemental and force connectivity inside the program and share with you 6) My whole aim is to get the reaction forces, - I expect the sum of forces and moments created by the external forces has to match with support force and moment sum ( Equilibrium). That's not happening- if I use the Statistics check; 7) It is perfectly Ok if you please try and suggest some approximate values for A, J, and ... other parameters in such a way that it could produce equilibrated boundary reaction forces. Sorry for the lengthy reply

[JWock82]

I've run your model, rendered it, and it looks unstable to me. There are a few issues that stick out to me:

1. What's to prevent the structure from rotating about the red axis passing through nodes 8 and 9 that I've sketched below?
2. The bigger question for me is how did this get through the analysis without flagging any instabilities? Pynite should have picked up on this. I'll have to have a deeper look at what's causing this modeling error to go unchecked.

[jaganmail]

Hi Craig, Thank you. By adding a rotational constraint, I managed to reproduce the expected boundary reactions, and the entire system now works perfectly. I wanted to express my gratitude for your assistance and the wonderful code you've provided. I initially overlooked the issue, but thanks to your guidance, I was able to identify the problem and implement the necessary corrections. I truly appreciate the time you took to help me and your generosity in sharing such a valuable resource as freeware. Your code has not only helped me resolve a critical issue but has also been instrumental in enhancing my understanding of the subject. Your dedication to creating high-quality and freely accessible tools is commendable. Once again, thank you for your support and for making this experience a positive one. Your efforts have made a significant difference in my project, and I'm truly grateful. Jagan

[JWock82]

I'm glad your problem is solved. I'm sorry it took so long for me to respond. I do this in my free time, and work has needed a lot of attention from me for the last several months, so it's been difficult for me to get around to it.

[jaganmail]

Thank you Craig Jagan