1# ============================================================================= 2# PROJECT CHRONO - http:#projectchrono.org 3# 4# Copyright (c) 2014 projectchrono.org 5# All rights reserved. 6# 7# Use of this source code is governed by a BSD-style license that can be found 8# in the LICENSE file at the top level of the distribution and at 9# http:#projectchrono.org/license-chrono.txt. 10# 11# ============================================================================= 12# Authors: Simone Benatti 13# ============================================================================= 14# 15# FEA demo on applying loads to beams 16# 17# ============================================================================= 18 19import pychrono as chrono 20import pychrono.fea as fea 21import pychrono.irrlicht as chronoirr 22import os 23import copy 24 25out_dir = chrono.GetChronoOutputPath() + "FEA_LOADS" # Output directory 26 27 28print("Copyright (c) 2017 projectchrono.org ") 29 30# Create (if needed) output directory 31if not os.path.isdir(out_dir): 32 os.mkdir(out_dir) 33 34 35# Create the physical system 36my_system = chrono.ChSystemSMC() 37 38# Create the Irrlicht visualization 39 40application = chronoirr.ChIrrApp(my_system, "Loads on beams", chronoirr.dimension2du(800, 600)) 41application.AddTypicalLogo() 42application.AddTypicalSky() 43application.AddTypicalLights() 44application.AddTypicalCamera(chronoirr.vector3df(0.5, 0.0, -3.0), chronoirr.vector3df(0.5, 0.0, 0.0)) 45 46# Create a mesh 47mesh = fea.ChMesh() 48my_system.Add(mesh) 49 50# Create some nodes (with default mass 0) 51nodeA = fea.ChNodeFEAxyzrot(chrono.ChFrameD(chrono.ChVectorD(0, 0, 0))) 52nodeB = fea.ChNodeFEAxyzrot(chrono.ChFrameD(chrono.ChVectorD(2, 0, 0))) 53nodeA.SetMass(0.0) 54nodeB.SetMass(0.0) 55mesh.AddNode(nodeA) 56mesh.AddNode(nodeB) 57 58# Create beam section & material 59msection = fea.ChBeamSectionEulerAdvanced() 60beam_wy = 0.1 61beam_wz = 0.2 62msection.SetAsRectangularSection(beam_wy, beam_wz) 63msection.SetYoungModulus(0.01e9) 64msection.SetGshearModulus(0.01e9 * 0.3) 65msection.SetBeamRaleyghDamping(0.200) 66msection.SetDensity(1500) 67 68# Create an Eulero-Bernoulli beam with a single element 69elementA = fea.ChElementBeamEuler() 70elementA.SetNodes(nodeA, nodeB) 71elementA.SetSection(msection) 72mesh.AddElement(elementA) 73 74# Create the ground body 75ground = chrono.ChBody() 76ground.SetBodyFixed(True) 77my_system.Add(ground) 78 79# Create a constraat the end of the beam 80constrA = chrono.ChLinkMateGeneric() 81constrA.Initialize(nodeA, ground, False, nodeA.Frame(), nodeA.Frame()) 82my_system.Add(constrA) 83constrA.SetConstrainedCoords(True, True, True, # x, y, z 84 True, True, True) # Rx, Ry, Rz 85 86# ----------------------------------------------------------------- 87# Apply loads 88 89# Create the load container 90loadcontainer = chrono.ChLoadContainer() 91my_system.Add(loadcontainer) 92 93 94#define LOAD_5 95 96 97print("Applied loads: \n") 98 99 100 101# Example 6: 102# As before, create a custom load with stiff force, acting on a single node, but 103# this time we inherit directly from ChLoadCustom, i.e. a load that does not require ChLoader features. 104# This is mostly used in case one does not need the automatic surface/volume quadrature of ChLoader. 105# As a stiff load, this will automatically generate a jacobian (tangent stiffness matrix K) 106# that will be used in statics, implicit integrators, etc. 107 108print(" Custom load with stiff force, acting on a single node (VER 2).") 109 110nodeD = fea.ChNodeFEAxyz(chrono.ChVectorD(2, 10, 3)) 111mesh.AddNode(nodeD) 112 113class MyLoadCustom(chrono.ChLoadCustom): 114 def __init__(self, mloadable): 115 chrono.ChLoadCustom.__init__(self, mloadable) 116 #/ "Virtual" copy constructor (covariant return type). 117 def Clone(self): 118 newinst = copy.deepcopy(self) 119 return newinst 120 # Compute Q=Q(x,v) 121 # This is the function that you have to implement. It should return the generalized Q load 122 # (i.e.the force in generalized lagrangian coordinates). 123 # For ChNodeFEAxyz, Q loads are expected as 3-rows vectors, containing absolute force x,y,z. 124 # As this is a stiff force field, dependency from state_x and state_y must be considered. 125 def ComputeQ(self,state_x, #/< state position to evaluate Q 126 state_w): #/< state speed to evaluate Q 127 if not state_x==None and not state_w==None : 128 node_pos = chrono.ChVectorD(state_x[0], state_x[1], state_x[2]) 129 node_vel = chrono.ChVectorD(state_w[0], state_w[1], state_w[2]) 130 else: 131 mynode = fea.CastToChNodeFEAxyz( fea.CastToChNodeFEAbase( chrono.CastToChNodeBase(self.loadable) )) 132 node_pos = mynode.GetPos() 133 node_vel = mynode.GetPos_dt() 134 # Just implement a simple force+spring+damper in xy plane, 135 # for spring&damper connected to absolute reference: 136 Kx = 100 137 Ky = 400 138 Dx = 0.6 139 Dy = 0.9 140 x_offset = 2 141 y_offset = 5 142 x_force = 50 143 y_force = 0 144 # Store the computed generalized forces in this.load_Q, same x,y,z order as in state_w 145 self.load_Q[0] = x_force - Kx * (node_pos.x - x_offset) - Dx * node_vel.x 146 self.load_Q[1] = y_force - Ky * (node_pos.y - y_offset) - Dy * node_vel.y 147 self.load_Q[2] = 0 148 # Remember to set this as stiff, to enable the jacobians 149 def IsStiff(self) : 150 return True 151 152# Instance load object, applying to a node, as in previous example, and add to container: 153mloadcustom = MyLoadCustom(nodeD) 154loadcontainer.Add(mloadcustom) 155 156 157 158# Stiffness and damping matrices 159K_matrix = chrono.ChMatrixDynamicD(6, 6) 160R_matrix = chrono.ChMatrixDynamicD(6, 6) 161 162for i in range(6): 163 K_matrix[i, i] = 1e5 164 165for i in range(6): 166 for j in range(6): 167 R_matrix[i, j] = 1e-3 168 169ch_bushing = fea.ChLoadXYZROTnodeBodyBushingGeneric( 170 nodeB, 171 ground, 172 nodeB, 173 K_matrix, 174 R_matrix) 175 176loadcontainer.Add(ch_bushing) 177 178# ----------------------------------------------------------------- 179# Set visualization of the FEM mesh. 180 181mvisualizebeamA = fea.ChVisualizationFEAmesh(mesh) 182mvisualizebeamA.SetFEMdataType(fea.ChVisualizationFEAmesh.E_PLOT_ELEM_BEAM_MZ) 183mvisualizebeamA.SetColorscaleMinMax(-400, 200) 184mvisualizebeamA.SetSmoothFaces(True) 185mvisualizebeamA.SetWireframe(False) 186mesh.AddAsset(mvisualizebeamA) 187 188mvisualizebeamC = fea.ChVisualizationFEAmesh(mesh) 189mvisualizebeamC.SetFEMglyphType(fea.ChVisualizationFEAmesh.E_GLYPH_NODE_CSYS) 190mvisualizebeamC.SetFEMdataType(fea.ChVisualizationFEAmesh.E_PLOT_NONE) 191mvisualizebeamC.SetSymbolsThickness(0.006) 192mvisualizebeamC.SetSymbolsScale(0.01) 193mvisualizebeamC.SetZbufferHide(False) 194mesh.AddAsset(mvisualizebeamC) 195 196application.AssetBindAll() 197application.AssetUpdateAll() 198 199# ----------------------------------------------------------------- 200 201# Setup a MINRES solver. For FEA one cannot use the default PSOR type solver. 202solver = chrono.ChSolverMINRES() 203my_system.SetSolver(solver) 204solver.SetMaxIterations(200) 205solver.SetTolerance(1e-15) 206solver.EnableDiagonalPreconditioner(True) 207solver.SetVerbose(False) 208 209my_system.SetSolverForceTolerance(1e-13) 210 211# Set integrator 212ts = chrono.ChTimestepperEulerImplicitLinearized(my_system) 213my_system.SetTimestepper(ts) 214 215# Simulation loop 216 217application.SetTimestep(1e-3) 218 219while (application.GetDevice().run()) : 220 application.BeginScene() 221 application.DrawAll() 222 application.DoStep() 223 if not application.GetPaused() : 224 time = my_system.GetChTime() 225 posB = nodeB.GetPos() 226 application.EndScene() 227 228