J04a - Follower Force |
 
|
J04a - Follower Force |
|
This example illustrates the application of arbitrary loading via the User Solver Variables feature. The direction of the applied force is dependent on the instantaneous position of the point of application. The use of custom code enables the definition of a force which can track/follow the direction of response, via continuous adjustment of the global force vector.
The model consists of a simple horizontal beam of length L, which is rigid and massless. It is constrained at one end in all degrees of freedom except for rotation in one degree of freedom, such that in plan view, it is free to rotate about a fixed pivot point. A point mass, m, is positioned at the free end of the beam. A force of magnitude F is applied to the free end in order to induce some rotation. Based on the instantaneous orientation of the beam, the direction of this force is dynamically updated such that it always remains perpendicular to the free end. The custom code uses an appropriate local to global transformation matrix to convert the local force terms into the global axis system before augmenting the global force vector. In order to verify correct model behaviour, results are compared with an analytical solution, and also an analogous model which applied a constant rotational moment at the fixed end.
The theoretical angular acceleration is defined as follows:
(1)
where
•a is acceleration in rads/s2
•F is the applied force in N
•m is the point mass in kg
•L is the beam length
Based on the above inputs, the angular acceleration at the free end should be 2.86degs/s2.
The rotational acceleration predicted by Flexcom for the customised model is exactly 2.86degs/s2 as expected.
Angular Acceleration at Fixed End
Pre-compiled DLLs for both 32-bit and 64-bit operating systems are provided with this example. The original Fortran source code is also provided should you wish to examine it.
subroutine user_solver_variables(iter,time,ramp,nnode,nncon,ncord,nndof,nelmn,necon,nedof, &
size_force,size_mass,size_stiff,elmcon,nodcon,intoun,ietoue,cord,displacement, &
velocity,acceleration,trgb,tglu,disp_prev,pos_prev,vel_prev,acc_prev,axial_force,y_shear, &
z_shear,torque,y_bending,z_bending,eff_tension,y_curvature,z_curvature,force,mass,stiff)
!dec$ attributes dllexport, stdcall, reference :: user_solver_variables
implicit none
! Input variables - cannot be modified within this subroutine
integer, intent(in) :: iter
!dec$ attributes value :: iter
real(8), intent(in) :: time
!dec$ attributes value :: time
real(8), intent(in) :: ramp
!dec$ attributes value :: ramp
integer(4), intent(in) :: nnode
!dec$ attributes value :: nnode
integer(4), intent(in) :: nncon
!dec$ attributes value :: nncon
integer(4), intent(in) :: ncord
!dec$ attributes value :: ncord
integer(4), intent(in) :: nndof
!dec$ attributes value :: nndof
integer(4), intent(in) :: nelmn
!dec$ attributes value :: nelmn
integer(4), intent(in) :: necon
!dec$ attributes value :: necon
integer(4), intent(in) :: nedof
!dec$ attributes value :: nedof
integer(4), intent(in) :: size_force
!dec$ attributes value :: size_force
integer(4), intent(in) :: size_mass
!dec$ attributes value :: size_mass
integer(4), intent(in) :: size_stiff
!dec$ attributes value :: size_stiff
integer(4), intent(in), dimension(nncon, nelmn) :: elmcon
!dec$ attributes reference :: elmcon
integer(4), intent(in), dimension(necon,nnode) :: nodcon
!dec$ attributes reference :: nodcon
integer(4), intent(in), dimension(nnode) :: intoun
!dec$ attributes reference :: intoun
integer(4), intent(in), dimension(nelmn) :: ietoue
!dec$ attributes reference :: ietoue
real(8), intent(in), dimension(ncord,nnode) :: cord
!dec$ attributes reference :: cord
real(8), intent(in), dimension(nndof,nnode) :: displacement
!dec$ attributes reference :: displacement
real(8), intent(in), dimension(nndof,nnode) :: velocity
!dec$ attributes reference :: velocity
real(8), intent(in), dimension(nndof,nnode) :: acceleration
!dec$ attributes reference :: acceleration
real(8), intent(in), dimension(3,3,nelmn) :: trgb
!dec$ attributes reference :: trgb
real(8), intent(in), dimension(3,3,nelmn) :: tglu
!dec$ attributes reference :: tglu
real(8), intent(in), dimension(nndof,nnode) :: disp_prev
!dec$ attributes reference :: disp_prev
real(8), intent(in), dimension(nndof,nnode) :: pos_prev
!dec$ attributes reference :: pos_prev
real(8), intent(in), dimension(nndof,nnode) :: vel_prev
!dec$ attributes reference :: vel_prev
real(8), intent(in), dimension(nndof,nnode) :: acc_prev
!dec$ attributes reference :: acc_prev
real(8), intent(in), dimension(nelmn,3) :: axial_force
!dec$ attributes reference :: axial_force
real(8), intent(in), dimension(nelmn,3) :: y_shear
!dec$ attributes reference :: y_shear
real(8), intent(in), dimension(nelmn,3) :: z_shear
!dec$ attributes reference :: z_shear
real(8), intent(in), dimension(nelmn,3) :: torque
!dec$ attributes reference :: torque
real(8), intent(in), dimension(nelmn,3) :: y_bending
!dec$ attributes reference :: y_bending
real(8), intent(in), dimension(nelmn,3) :: z_bending
!dec$ attributes reference :: z_bending
real(8), intent(in), dimension(nelmn,3) :: eff_tension
!dec$ attributes reference :: eff_tension
real(8), intent(in), dimension(nelmn,3) :: y_curvature
!dec$ attributes reference :: y_curvature
real(8), intent(in), dimension(nelmn,3) :: z_curvature
!dec$ attributes reference :: z_curvature
! Output variables - can be modified within this subroutine if required
real(8), intent(inout), dimension(size_force,nedof) :: force
!dec$ attributes reference :: force
real(8), intent(inout), dimension(nedof,nedof,size_mass) :: mass
!dec$ attributes reference :: mass
real(8), intent(inout), dimension(nedof,nedof,size_stiff) :: stiff
!dec$ attributes reference :: stiff
! Variable names
! iter : Current iteration
! time : Current timestep
! ramp : Ramp
! nnode : Number of nodes in the model
! nncon : Number of nodes with connected elements
! ncord : Number of coordinates
! nndof : Number of degrees of freedom per node (6)
! nelmn : Number of elements in the model
! necon : Number of elements connected
! nedof : Number of degrees of freedom per element (14)
! size_force : Dimension of the global force vector
! size_mass : Dimension of the global mass matrix
! size_stiff : Dimension of the global stiffness matrix
! elmcon : Element connectivity array
! nodcon : Node connectivity array
! intoun : Internal node to user node numbering array
! ietoue : Internal element to user element numbering array
! cord : Initial nodal co-ordinates
! displacement : Nodal displacements at previous iteration
! velocity : Nodal velocities at previous iteration
! acceleration : Nodal accelerations at previous iteration
! trgb : Rigid body rotation (local undeformed -> convected) transformation matrix
! tglu : Global to local undeformed transformation matrix
! disp_prev : Nodal displacements at previous timestep
! pos_prev : Nodal positions at previous timestep
! vel_prev : Nodal velocities at previous timestep
! acc_prev : Nodal accelerations at previous timestep
! axial_force : Axial force in elements at previous timestep
! y_shear : Y Shear forces in elements at previous timestep
! z_shear : Z Shear forces in elements at previous timestep
! torque : Torque in elements at previous timestep
! y_bending : Y bending moments in elements at previous timestep
! z_bending : Z bending moments in elements at previous timestep
! eff_tension : Effective Tension in elements at previous timestep
! y_curvature : Y curvatures in elements at previous timestep
! z_curvature : Z curvatures in elements at previous timestep
! force : Global force vector at previous iteration
! mass : Global mass matrix
! stiff : Global stiffness matrix
! Declare local variables.
integer i,j, index, index1, index2
integer :: int_elem_no, int_node_no
integer :: node_on_elem_index
integer :: int_DOF
integer :: user_elem_no, user_node_no, user_dof
real(8), dimension(3) :: user_applied_force
real(8), dimension(3) :: global_user_applied_force
real(8),dimension(3,3) :: trgb_int_elem
real(8),dimension(3,3) :: tglu_int_elem
real(8),dimension(3,3) :: transformation_matrix
real(8),dimension(3,3) :: transpose_matrix
real(8), dimension(size_force,nedof) :: force_array
logical :: is_elem_OK, is_node_OK, is_node_on_element
! Follower Force Test
! User data
user_elem_no = 10
user_node_no = 11
! Local axes applied force
user_applied_force(1:3) = [0.0d0, 0.0d0, 1000.0d0] ! Force applied in DOF3
! Internal works
int_elem_no = 0
int_node_no = 0
transformation_matrix = 0.0d0
is_elem_OK = .false.
is_node_OK = .false.
is_node_on_element = .false.
! Check if user element exists in Flexcom list of elements and set internal element number
! if element exists.
do i = 1,nelmn
if( user_elem_no == ietoue(i) )then
int_elem_no = i
is_elem_OK = .true.
end if
end do
! Check if user node exists in Flexcom list of elements and set internal node number
! if node exists.
do i = 1, nnode
if( user_node_no == intoun(i) )then
int_node_no = i
is_node_OK = .true.
end if
end do
! Check if internal node is on internal element
do i=1,nncon
if(int_node_no == elmcon(i,int_elem_no))then
is_node_on_element = .true.
exit
end if
end do
if( is_node_on_element == .false.) then
write(*,*) "Node ", user_node_no, "is not included on any elements."
end if
! Apply user force at internal node number and internal DOF
if( is_elem_OK .and. is_node_OK .and. is_node_on_element ) then
!Extract transformation matrices at internal element number
trgb_int_elem(1:3,1:3) = trgb(1:3,1:3, int_elem_no)
tglu_int_elem(1:3,1:3) = tglu(1:3,1:3, int_elem_no)
! Perform transformation from local to global axes
transformation_matrix(1:3,1:3) = matmul(trgb_int_elem(1:3,1:3),tglu_int_elem(1:3,1:3))
! Transformation matrix transpose
transpose_matrix = transpose(transformation_matrix)
! Global user applied force
global_user_applied_force(1:3) = matmul(transpose_matrix(1:3,1:3),user_applied_force(1:3))
! Set node on element index
node_on_elem_index = 0
if ( int_node_no == elmcon(2,int_elem_no) )node_on_elem_index = 6
! Augment the force with global user applied force.
! DOF 1-3 for Start Node, DOF 7-9 for End Node
do i = 1,3
int_DOF = i + node_on_elem_index
force(int_elem_no,int_DOF) = force(int_elem_no,int_DOF) + ramp * global_user_applied_force(i)
end do
end if
! Do not alter the next line.
end subroutine user_solver_variables
Source Code