From 0829f84a941bb3282baeadc639e100271700567d Mon Sep 17 00:00:00 2001 From: Bonnie Jonkman Date: Mon, 14 Aug 2017 13:40:22 -0600 Subject: [PATCH] BD: fixed local sectional loads and other issues - fixed initialization of PointLoad mesh (driver only) - added variables and reorganized loops for optimization - fixed internal-loads outputs (fixes https://github.com/OpenFAST/openfast/issues/11) - Note that https://github.com/OpenFAST/openfast/issues/13 has also been fixed (in a previous commit) --- modules-local/beamdyn/src/BeamDyn.f90 | 1443 ++++++++++------- modules-local/beamdyn/src/BeamDyn_IO.f90 | 14 +- modules-local/beamdyn/src/BeamDyn_Subs.f90 | 42 +- modules-local/beamdyn/src/Driver_Beam.f90 | 4 +- .../beamdyn/src/Registry_BeamDyn.txt | 16 +- 5 files changed, 920 insertions(+), 599 deletions(-) diff --git a/modules-local/beamdyn/src/BeamDyn.f90 b/modules-local/beamdyn/src/BeamDyn.f90 index 3b1a65b445..762684d2f7 100644 --- a/modules-local/beamdyn/src/BeamDyn.f90 +++ b/modules-local/beamdyn/src/BeamDyn.f90 @@ -2,7 +2,7 @@ ! LICENSING ! Copyright (C) 2015-2016 National Renewable Energy Laboratory ! Copyright (C) 2016-2017 Envision Energy USA, LTD -! +! ! Licensed under the Apache License, Version 2.0 (the "License"); ! you may not use this file except in compliance with the License. ! You may obtain a copy of the License at @@ -533,26 +533,60 @@ SUBROUTINE BD_InitShpDerJaco( GLL_Nodes, p ) TYPE(BD_ParameterType), INTENT(INOUT) :: p !< Parameters REAL(BDKi) :: Gup0(3) - INTEGER(IntKi) :: inode - INTEGER(IntKi) :: i - INTEGER(IntKi) :: j + INTEGER(IntKi) :: i, j + INTEGER(IntKi) :: nelem, idx_qp CHARACTER(*), PARAMETER :: RoutineName = 'BD_InitShpDerJaco' CALL BD_diffmtc(p,GLL_nodes,p%Shp,p%ShpDer) - p%Jacobian(:,:) = 0.0_BDKi - DO i = 1,p%elem_total - DO j = 1, p%nqp + DO nelem = 1,p%elem_total + DO idx_qp = 1, p%nqp Gup0(:) = 0.0_BDKi - DO inode=1,p%nodes_per_elem - Gup0(:) = Gup0(:) + p%ShpDer(inode,j)*p%uuN0(1:3,inode,i) + DO i=1,p%nodes_per_elem + Gup0(:) = Gup0(:) + p%ShpDer(i,idx_qp)*p%uuN0(1:3,i,nelem) ENDDO - p%Jacobian(j,i) = TwoNorm(Gup0) + p%Jacobian(idx_qp,nelem) = TwoNorm(Gup0) ENDDO ENDDO + + ! save some variables so we don't recalculate them so often: + DO nelem=1,p%elem_total + DO idx_qp=1,p%nqp + DO j=1,p%nodes_per_elem + DO i=1,p%nodes_per_elem + p%QPtw_Shp_Shp_Jac( idx_qp,i,j,nelem) = p%Shp( i,idx_qp)*p%Shp( j,idx_qp)*p%QPtWeight(idx_qp)*p%Jacobian(idx_qp,nelem) + p%QPtw_ShpDer_ShpDer_Jac(idx_qp,i,j,nelem) = p%ShpDer(i,idx_qp)*p%ShpDer(j,idx_qp)*p%QPtWeight(idx_qp)/p%Jacobian(idx_qp,nelem) + END DO + END DO + END DO + END DO + + + DO idx_qp=1,p%nqp + DO j=1,p%nodes_per_elem + DO i=1,p%nodes_per_elem + p%QPtw_Shp_ShpDer(idx_qp,i,j) = p%Shp(i,idx_qp)*p%ShpDer(j,idx_qp)*p%QPtWeight(idx_qp) + END DO + END DO + END DO + + DO nelem=1,p%elem_total + DO i=1,p%nodes_per_elem + DO idx_qp=1,p%nqp + p%QPtw_Shp_Jac(idx_qp,i,nelem) = p%Shp( i,idx_qp)*p%QPtWeight(idx_qp)*p%Jacobian(idx_qp,nelem) + END DO + END DO + END DO + + DO i=1,p%nodes_per_elem + DO idx_qp=1,p%nqp + p%QPtw_ShpDer(idx_qp,i) = p%ShpDer(i,idx_qp)*p%QPtWeight(idx_qp) + END DO + END DO + END SUBROUTINE BD_InitShpDerJaco @@ -680,6 +714,13 @@ subroutine SetParameters(InitInp, InputFileData, p, ErrStat, ErrMsg) CALL AllocAry(p%ShpDer, p%nodes_per_elem,p%nqp, 'p%ShpDer', ErrStat2,ErrMsg2); CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) CALL AllocAry(p%Jacobian,p%nqp, p%elem_total,'p%Jacobian',ErrStat2,ErrMsg2); CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + CALL AllocAry(p%QPtw_Shp_Shp_Jac ,p%nqp,p%nodes_per_elem,p%nodes_per_elem,p%elem_total,'p%QPtw_Shp_Shp_Jac', ErrStat2,ErrMsg2); CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + CALL AllocAry(p%QPtw_Shp_ShpDer ,p%nqp,p%nodes_per_elem,p%nodes_per_elem, 'p%QPtw_Shp_ShpDer', ErrStat2,ErrMsg2); CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + CALL AllocAry(p%QPtw_ShpDer_ShpDer_Jac,p%nqp,p%nodes_per_elem,p%nodes_per_elem,p%elem_total,'p%QPtw_ShpDer_ShpDer_Jac',ErrStat2,ErrMsg2); CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + CALL AllocAry(p%QPtw_Shp_Jac ,p%nqp ,p%nodes_per_elem,p%elem_total,'p%QPtw_Shp_Jac', ErrStat2,ErrMsg2); CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + CALL AllocAry(p%QPtw_ShpDer ,p%nqp ,p%nodes_per_elem, 'p%QPtw_ShpDer', ErrStat2,ErrMsg2); CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + + CALL AllocAry(p%uuN0, p%dof_node,p%nodes_per_elem, p%elem_total,'uuN0 (initial position) array',ErrStat2,ErrMsg2); CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) CALL AllocAry(p%rrN0, (p%dof_node/2),p%nodes_per_elem, p%elem_total,'p%rrN0',ErrStat2,ErrMsg2); CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) CALL AllocAry(p%uu0, p%dof_node, p%nqp, p%elem_total,'p%uu0', ErrStat2,ErrMsg2); CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) @@ -708,8 +749,8 @@ subroutine SetParameters(InitInp, InputFileData, p, ErrStat, ErrMsg) ! Some subroutines are looking at a single element, in which case the values stored in p%nodes_elem_idx ! are used to indicate which node to start with. DO i=1,p%elem_total - p%node_elem_idx(i,1) = (i-1)*(p%nodes_per_elem-1)+1 ! First node in element - p%node_elem_idx(i,2) = i *(p%nodes_per_elem-1) ! Last node in element + p%node_elem_idx(i,1) = (i-1)*(p%nodes_per_elem-1) + 1 ! First node in element + p%node_elem_idx(i,2) = i *(p%nodes_per_elem-1) + 1 ! Last node in element ENDDO @@ -723,8 +764,8 @@ subroutine SetParameters(InitInp, InputFileData, p, ErrStat, ErrMsg) p%OutNd2NdElem(:,1) = 1 ! note this is an array indx = 2 DO i=1,p%elem_total - DO j=2,p%nodes_per_elem ! GLL nodes overlap at element end points; we will skip the first node of each element (after the first one) - p%NdIndx(indx) = (i-1)*p%nodes_per_elem + j ! Index into BldMotion mesh (to number the nodes for output without using collocated nodes) + DO j=2,p%nodes_per_elem ! GLL nodes overlap at element end points; we will skip the first node of each element (after the first one) + p%NdIndx(indx) = (i-1)*p%nodes_per_elem + j ! Index into BldMotion mesh (to number the nodes for output without using collocated nodes) p%OutNd2NdElem(1,indx) = j ! Node number. To go from an output node number to a node/elem pair p%OutNd2NdElem(2,indx) = i ! Element number. To go from an output node number to a node/elem pair indx = indx + 1 @@ -912,11 +953,8 @@ subroutine Init_u( InitInp, p, u, ErrStat, ErrMsg ) real(R8Ki) :: DCM(3,3) ! must be same type as mesh orientation fields real(ReKi) :: Pos(3) ! must be same type as mesh position fields - real(BDKi) :: temp_POS(3) - integer(intKi) :: temp_id integer(intKi) :: i,j ! loop counters - integer(intKi) :: j_start ! loop counter start for looping over quadrature points integer(intKi) :: NNodes ! number of nodes in mesh integer(intKi) :: ErrStat2 ! temporary Error status @@ -1031,7 +1069,7 @@ subroutine Init_u( InitInp, p, u, ErrStat, ErrMsg ) DO i=1,p%elem_total DO j=1,p%nodes_per_elem - temp_POS = p%GlbPos(1:3) + MATMUL(p%GlbRot,p%uuN0(1:3,j,i)) + POS = p%GlbPos(1:3) + MATMUL(p%GlbRot,p%uuN0(1:3,j,i)) ! Note: Here we can use this subroutine to get the DCM. This is under the assumption ! that there is no rotational displacement yet, so x%q is zero @@ -1072,6 +1110,7 @@ subroutine Init_u( InitInp, p, u, ErrStat, ErrMsg ) !................................. NNodes = p%nqp*p%elem_total + 2*p%qp_indx_offset + CALL MeshCreate( BlankMesh = u%DistrLoad & ,IOS = COMPONENT_INPUT & ,NNodes = NNodes & @@ -1184,7 +1223,6 @@ subroutine Init_MiscVars( p, u, y, m, ErrStat, ErrMsg ) - !if (p%analysis_type == BD_DYNAMIC_ANALYSIS) then ! allocate arrays at initialization so we don't waste so much time on memory allocation/deallocation on each step ! Arrays for use with LAPACK matrix solve routines (A*X = B solve for X) @@ -1205,11 +1243,13 @@ subroutine Init_MiscVars( p, u, y, m, ErrStat, ErrMsg ) CALL AllocAry(m%MassM, p%dof_node,p%node_total,p%dof_node,p%node_total, 'MassM', ErrStat2,ErrMsg2); CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) CALL AllocAry(m%DampG, p%dof_node,p%node_total,p%dof_node,p%node_total, 'DampG', ErrStat2,ErrMsg2); CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) - CALL AllocAry(m%BldInternalForce, p%dof_node,y%BldMotion%NNodes, 'Blade Internal Force info', ErrStat2,ErrMsg2); CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + CALL AllocAry(m%BldInternalForceFE, p%dof_node,p%node_total, 'Blade Internal Force info', ErrStat2,ErrMsg2); CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) CALL AllocAry(m%Nrrr, (p%dof_node/2),p%nodes_per_elem,p%elem_total,'Nrrr: rotation parameters relative to root', ErrStat2,ErrMsg2); CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) - CALL AllocAry(m%elf, p%dof_node,p%nodes_per_elem,p%elem_total, 'elf', ErrStat2,ErrMsg2); CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + CALL AllocAry(m%elf, p%dof_node,p%nodes_per_elem, 'elf', ErrStat2,ErrMsg2); CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + + CALL AllocAry(m%EFint, p%dof_node,p%nodes_per_elem,p%elem_total, 'Elastic Force internal', ErrStat2,ErrMsg2); CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) ! Note the index ordering here. This comes from the reshaping to other arrays used with LAPACK solvers CALL AllocAry(m%elk, p%dof_node,p%nodes_per_elem,p%dof_node,p%nodes_per_elem, 'elk', ErrStat2,ErrMsg2); CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) @@ -1221,13 +1261,14 @@ subroutine Init_MiscVars( p, u, y, m, ErrStat, ErrMsg ) ! Array for storing the position information for the quadrature points. -! CALL AllocAry(m%qp%q, p%dof_node,p%nqp,p%elem_total, 'm%qp%q -- position info at quadrature point',ErrStat2,ErrMsg2); CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) CALL AllocAry(m%qp%uuu, p%dof_node ,p%nqp,p%elem_total, 'm%qp%uuu displacement at quadrature point',ErrStat2,ErrMsg2); CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) CALL AllocAry(m%qp%uup, p%dof_node/2,p%nqp,p%elem_total, 'm%qp%uup displacement prime at quadrature point',ErrStat2,ErrMsg2); CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) CALL AllocAry(m%qp%vvv, p%dof_node ,p%nqp,p%elem_total, 'm%qp%vvv velocity at quadrature point',ErrStat2,ErrMsg2); CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) CALL AllocAry(m%qp%vvp, p%dof_node ,p%nqp,p%elem_total, 'm%qp%vvp velocity prime at quadrature point',ErrStat2,ErrMsg2); CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) CALL AllocAry(m%qp%aaa, p%dof_node ,p%nqp,p%elem_total, 'm%qp%aaa acceleration at quadrature point',ErrStat2,ErrMsg2); CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) - + ! Need to initialize this to zero especially for static cases + m%qp%aaa = 0.0_BDKi + ! E1, kappa -- used in force calculations CALL AllocAry(m%qp%E1, p%dof_node/2,p%nqp,p%elem_total, 'm%qp%E1 at quadrature point',ErrStat2,ErrMsg2); CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) CALL AllocAry(m%qp%kappa, p%dof_node/2,p%nqp,p%elem_total, 'm%qp%kappa at quadrature point',ErrStat2,ErrMsg2); CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) @@ -1266,7 +1307,9 @@ subroutine Init_MiscVars( p, u, y, m, ErrStat, ErrMsg ) CALL AllocAry(m%qp%Xd, 6,6, p%nqp,p%elem_total, 'm%qp%Xd term at quadrature point',ErrStat2,ErrMsg2); CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) CALL AllocAry(m%qp%Yd, 6,6, p%nqp,p%elem_total, 'm%qp%Yd term at quadrature point',ErrStat2,ErrMsg2); CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) - !end if + ! Initialize variables that could be output before they are set. + m%qp%Fi = 0.0_BDKi ! This could be output before it gets set. + ! create copy of u%DistrLoad at y%BldMotion locations (for WriteOutput only) if (p%OutInputs) then @@ -1574,11 +1617,11 @@ SUBROUTINE BD_CalcOutput( t, u, p, x, xd, z, OtherState, y, m, ErrStat, ErrMsg ) CALL BD_CalcForceAcc(x_tmp, m%u, p, m, ErrStat2,ErrMsg2) CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) - +!FIXME: First node pointload is missing. Same with first QP node and DistrLoad. y%ReactionForce%Force(:,1) = -MATMUL(p%GlbRot,m%RHS(1:3,1)) y%ReactionForce%Moment(:,1) = -MATMUL(p%GlbRot,m%RHS(4:6,1)) - CALL BD_GenerateInternalForceMoment( y, p, m ) + CALL BD_InternalForceMoment( x, p, m ) ELSE m%RHS = 0.0_BdKi ! accelerations are set to zero in the static case @@ -1589,10 +1632,12 @@ SUBROUTINE BD_CalcOutput( t, u, p, x, xd, z, OtherState, y, m, ErrStat, ErrMsg ) ENDDO - CALL BD_GenerateInternalForceMoment( y, p, m ) - - y%ReactionForce%Force(:,1) = -m%BldInternalForce(1:3,1) - y%ReactionForce%Moment(:,1) = -m%BldInternalForce(4:6,1) + CALL BD_InternalForceMoment( x, p, m ) + + !FIXME: Check that this works for the static case. The internal force for dynamic case uses the reaction force for the first node. + ! Get the root force from the first finite element node. + y%ReactionForce%Force(:,1) = m%BldInternalForceFE(1:3,1) + y%ReactionForce%Moment(:,1) = m%BldInternalForceFE(4:6,1) ENDIF @@ -2103,19 +2148,14 @@ END SUBROUTINE BD_QPData_mEta_rho !! !! The equations used here can be found in the NREL publication CP-2C00-60759 !! (http://www.nrel.gov/docs/fy14osti/60759.pdf) -SUBROUTINE BD_ElasticForce(nelem,idx_qp,p,m,fact) +SUBROUTINE BD_ElasticForce(nelem,p,m,fact) INTEGER(IntKi), INTENT(IN ) :: nelem !< number of current element - INTEGER(IntKi), INTENT(IN ) :: idx_qp !< Index to quadrature point currently being calculated TYPE(BD_ParameterType), INTENT(IN ) :: p !< Parameters TYPE(BD_MiscVarType), INTENT(INOUT) :: m !< Misc/optimization variables. LOGICAL, INTENT(IN ) :: fact REAL(BDKi) :: cet !< for storing the \f$ I_{yy} + I_{zz} \f$ inertia term - REAL(BDKi) :: eee(6) !< intermediate array for calculation Strain and curvature terms of Fc - REAL(BDKi) :: fff(6) !< intermediate array for calculation of the elastic force, Fc - REAL(BDKi) :: Wrk(3) - REAL(BDKi) :: e1s REAL(BDKi) :: k1s REAL(BDKi) :: Wrk33(3,3) REAL(BDKi) :: tildeE(3,3) @@ -2123,180 +2163,206 @@ SUBROUTINE BD_ElasticForce(nelem,idx_qp,p,m,fact) REAL(BDKi) :: epsi(3,3) REAL(BDKi) :: mu(3,3) + INTEGER(IntKi) :: idx_qp !< Index to quadrature point currently being calculated - !> ### Calculate the 1D strain, \f$ \underline{\epsilon} \f$, equation (5) - !! \f$ \underline{\epsilon} = \underline{x}^\prime_0 + \underline{u}^\prime - - !! \left(\underline{\underline{R}}\underline{\underline{R}}_0\right) \bar{\imath}_1 - !! = \underline{E}_1 - - !! \left(\underline{\underline{R}}\underline{\underline{R}}_0\right) \bar{\imath}_1 \f$ - !! where \f$ \bar{\imath}_1 \f$ is the unit vector along the \f$ x_1 \f$ direction in the inertial frame - !! (z in the IEC). The last term simplifies to the first column of \f$ \underline{\underline{R}}\underline{\underline{R}}_0 \f$ - !! - !! Note: \f$ \underline{\underline{R}}\underline{\underline{R}}_0 \f$ is used to go from the material basis into the inertial basis - !! and the transpose for the other direction. - eee(1:3) = m%qp%E1(1:3,idx_qp,nelem) - m%qp%RR0(1:3,3,idx_qp,nelem) ! Using RR0 z direction in IEC coords + + if (.not. fact) then + + do idx_qp=1,p%nqp + call Calc_Fc_Fd() + end do + + else + + do idx_qp=1,p%nqp + + call Calc_Fc_Fd() - !> ### Set the 1D sectional curvature, \f$ \underline{\kappa} \f$, equation (5) - !! \f$ \underline{\kappa} = \underline{k} + \underline{\underline{R}}\underline{k}_i \f$ - !! where \f$ \underline{k} = \text{axial}\left(\underline{\underline{R}}^\prime\underline{\underline{R}}^T \right) \f$ - !! and \f$ \underline{k}_i \f$ is the corresponding initial curvature vector (root reference). - !! Note that \f$ \underline{k} \f$ can be calculated with rotation parameters as - !! \f$ \underline{k} = \underline{\underline{H}}(\underline{p}) \underline{p}' \f$ - !! - !! \f$ \text{axial}\left( \underline{\underline{A}} \right) \f$ is defined as \n - !! \f$ \text{axial}\left( \underline{\underline{A}} \right) - !! = \left\{ \begin{array}{c} a_1 \\ - !! a_2 \\ - !! a_3 \end{array}\right\} - !! = \left\{ \begin{array}{c} A_{32} - A_{23} \\ - !! A_{13} - A_{31} \\ - !! A_{21} - A_{12} \end{array}\right\} - !! \f$ - !! In other words, \f$ \tilde{k} = \left(\underline{\underline{R}}^\prime\underline{\underline{R}}^T \right) \f$. - !! Note: \f$ \underline{\kappa} \f$ was already calculated in the BD_DisplacementQP routine - eee(4:6) = m%qp%kappa(1:3,idx_qp,nelem) - - -!FIXME: note that the k_i terms may not be documented correctly here. - !> ### Calculate the elastic force, \f$ \underline{\mathcal{F}}^C \f$ - !! Using equations (15), (23), (24), and (5) we write the elastic force as\n - !! \f$ - !! \underline{\mathcal{F}}^C - !! = \left\{ \begin{array}{c} \underline{F} \\ - !! \underline{M} \end{array} \right\} - !! = \left(\underline{\underline{R}}\underline{\underline{R}}_0\right) - !! \underline{\underline{C}} - !! \left(\underline{\underline{R}}\underline{\underline{R}}_0\right)^T - !! \left\{ \begin{array}{c} \underline{\epsilon} \\ - !! \underline{\kappa} \end{array} \right\} - !! \f$ - !! - !! Note then that the extension twist term is added to this so that we get - !! \f$ \underline{\mathcal{F}}^C = \underline{F}^c_{a} + \underline{F}^c_{et} \f$ - !! - !! where \f$ \underline{F}^c_{et} \f$ is the extension twist coupling term. - - !> ###Calculate the first term. - !! \f$ \underline{F}^c_a - !! = \left(\underline{\underline{R}}\underline{\underline{R}}_0\right) - !! \underline{\underline{C}} - !! \left(\underline{\underline{R}}\underline{\underline{R}}_0\right)^T - !! \left\{ \begin{array}{c} - !! \underline{\epsilon} \\ - !! \underline{k} - !! \end{array} \right\} \f$ - !! - fff(1:6) = MATMUL(m%qp%Stif(:,:,idx_qp,nelem),eee) - - - !> ###Calculate the extension twist coupling. - !! This is calculated in the material basis, not in the rotated reference frame \f$ \underline{F}^c_a \f$ was calculated in. - !! First find \f$ \epsilon_1 \f$ and \f$ \kappa_1 \f$ in the material basis\n - !! \f$ \epsilon_{m} = \left( \underline{\underline{R}}\underline{\underline{R}}_0 \right) ^T \epsilon \f$ \n - !! and \n - !! \f$ \kappa_{m} = \left( \underline{\underline{R}}\underline{\underline{R}}_0 \right) ^T \underline{k}\f$ \n - - ! Strain into the material basis (eq (39) of Dymore manual) - Wrk(:) = MATMUL(TRANSPOSE(m%qp%RR0(:,:,idx_qp,nelem)),eee(1:3)) - e1s = Wrk(3) !epsilon_{1} in material basis (for major axis of blade, which is z in the IEC formulation) - - Wrk(:) = MATMUL(TRANSPOSE(m%qp%RR0(:,:,idx_qp,nelem)),eee(4:6)) - k1s = Wrk(3) !kappa_{1} in material basis (for major axis of blade, which is z in the IEC formulation) - - - !> Add extension twist coupling terms to the \f$ \underline{F}^c_{a} \f$\n - !! \f$ \underline{F}^c = - !! \underline{F}^c_{a} + \underline{F}^c_{et} \f$\n - !! The extension twist term is calculated in the material basis, which can be simplified and rewritten as follows (_the details of this are fuzzy to me_):\n - !! \f$ \underline{F}^c_{et} = - !! \begin{bmatrix} - !! \frac{1}{2} C_{et} \kappa^2_{m}(1) \left(\underline{\underline{R}}\underline{\underline{R}}_0\right)_{:,1} \\ - !! C_{et} \kappa_{m}(1) \epsilon_{m}(1) \left(\underline{\underline{R}}\underline{\underline{R}}_0\right)_{:,1} - !! \end{bmatrix} - !! \f$ - !! - !! Where \f$ C_{et} = C_{5,5} + C_{6,6} \f$ is the inertial term along major axis, \f$x\f$, in the material coordinate system - !! Note that with coverting to the FAST / IEC coordinate system, we now are using the Ixx and Iyy terms which are located at - !! \f$ C_{et} = C_{4,4} + C_{5,5} \f$ - ! Refer Section 1.4 in "Dymore User's Manual - Formulation and finite element implementation of beam elements". - cet= p%Stif0_QP(4,4,(nelem-1)*p%nqp+idx_qp) + p%Stif0_QP(5,5,(nelem-1)*p%nqp+idx_qp) ! Dymore theory (22) - m%qp%Fc(1:3,idx_qp,nelem) = fff(1:3) + 0.5_BDKi*cet*k1s*k1s*m%qp%RR0(1:3,3,idx_qp,nelem) ! Dymore theory (25a). Note z-axis is the length of blade. - m%qp%Fc(4:6,idx_qp,nelem) = fff(4:6) + cet*e1s*k1s*m%qp%RR0(1:3,3,idx_qp,nelem) ! Dymore theory (25b). Note z-axis is the length of blade. - - !> ###Calculate \f$ \underline{\mathcal{F}}^d \f$, equation (16) - !! \f$ \underline{F}^d = - !! \underline{\underline{\Upsilon}} \underline{\mathcal{F}}^c - !! = \begin{bmatrix} \underline{\underline{0}} & \underline{\underline{0}} \\ - !! \tilde{E}_1^T & \underline{\underline{0}} \end{bmatrix} - !! \underline{\mathcal{F}}^c - !! = \begin{bmatrix} \underline{0} \\ - !! \left(\underline{\mathcal{F}}^c \times \underline{E}_1 \right)^T - !! \end{bmatrix} \f$ - m%qp%Fd(1:3,idx_qp,nelem) = 0.0_BDKi -! ADP uu0 ref: If E1 is referenced against a different curve than Stif0_QP, there will be strange coupling terms here. - m%qp%Fd(4:6,idx_qp,nelem) = cross_product(m%qp%Fc(1:3,idx_qp,nelem), m%qp%E1(:,idx_qp,nelem)) + !> ###Calculate the \f$ \underline{\underline{\mathcal{O}}} \f$ from equation (19) + !! + !! \f$ \underline{\underline{\mathcal{O}}} = + !! \begin{bmatrix} + !! \underline{\underline{0}} & \underline{\underline{C}}_{11} \tilde{E}_1 - \tilde{F} \\ + !! \underline{\underline{0}} & \underline{\underline{C}}_{21} \tilde{E}_1 - \tilde{M} + !! \end{bmatrix} + !! = \begin{bmatrix} + !! \underline{\underline{0}} & \psi_E - \tilde{F} \\ + !! \underline{\underline{0}} & \mu - \tilde{M} + !! \end{bmatrix} + !! \f$ + Wrk33(:,:) = OuterProduct(m%qp%RR0(1:3,3,idx_qp,nelem), m%qp%RR0(1:3,3,idx_qp,nelem)) ! z-direction in IEC coords + C21(:,:) = m%qp%Stif(4:6,1:3,idx_qp,nelem) + cet*k1s*Wrk33(:,:) + + tildeE = SkewSymMat(m%qp%E1(:,idx_qp,nelem)) + epsi(:,:) = MATMUL(m%qp%Stif(1:3,1:3,idx_qp,nelem),tildeE) ! Stif is RR0 * p%Stif0_QP * RR0^T + mu(:,:) = MATMUL(C21,tildeE) + + m%qp%Oe(:,:,idx_qp,nelem) = 0.0_BDKi + m%qp%Oe(1:3,4:6,idx_qp,nelem) = epsi(1:3,1:3) - SkewSymMat(m%qp%Fc(1:3,idx_qp,nelem)) + m%qp%Oe(4:6,4:6,idx_qp,nelem) = mu(1:3,1:3) - SkewSymMat(m%qp%Fc(4:6,idx_qp,nelem)) + + + !> ###Calculated \f$ \underline{\underline{\mathcal{P}}} \f$ from equation (20) + !! + !! \f$ \underline{\underline{\mathcal{P}}} = + !! \begin{bmatrix} + !! \underline{\underline{0}} & \underline{\underline{0}} \\ + !! \left(\underline{\underline{\bar{C}}}_{11} \tilde{E}_1 \right)^T + \tilde{F} + !! \left(\underline{\underline{\bar{C}}}_{11} \tilde{E}_1 \right)^T + !! \end{bmatrix} + !! = \begin{bmatrix} + !! \underline{\underline{0}} & \underline{\underline{0}} \\ + !! \psi_E^T + \tilde{F} & \mu^T + !! \end{bmatrix} \f$ + m%qp%Pe(:,:,idx_qp,nelem) = 0.0_BDKi + m%qp%Pe(4:6,1:3,idx_qp,nelem) = TRANSPOSE(epsi) + SkewSymMat(m%qp%Fc(1:3,idx_qp,nelem)) + m%qp%Pe(4:6,4:6,idx_qp,nelem) = TRANSPOSE(mu) + + !> ###Calculated \f$ \underline{\underline{\mathcal{Q}}} \f$ from equation (21) + !! + !! \f{eqnarray*}{ + !! \underline{\underline{\mathcal{Q}}} + !! & =& \underline{\underline{\Upsilon}} \underline{\underline{\mathcal{O}}} + !! = \begin{bmatrix} 0 & 0 \\ + !! \tilde{E}_1^T & 0 \end{bmatrix} + !! \underline{\underline{\mathcal{O}}} \\ + !! \begin{bmatrix} 0 & 0 \\ + !! 0 & \underline{\underline{\mathcal{Q}}}_{22} \end{bmatrix} + !! & =& \tilde{E}_1^T \underline{\underline{\mathcal{O}}}_{12} + !! = - \tilde{E}_1 \underline{\underline{\mathcal{O}}}_{12} + !! \f}\n + !! Note: \f$ \tilde{E}_1^T = - \tilde{E}_1 \f$ + m%qp%Qe(:,:,idx_qp,nelem) = 0.0_BDKi + m%qp%Qe(4:6,4:6,idx_qp,nelem) = -MATMUL(tildeE,m%qp%Oe(1:3,4:6,idx_qp,nelem)) + end do + + ENDIF +contains + subroutine Calc_Fc_Fd() + REAL(BDKi) :: e1s + REAL(BDKi) :: eee(6) !< intermediate array for calculation Strain and curvature terms of Fc + REAL(BDKi) :: fff(6) !< intermediate array for calculation of the elastic force, Fc + !REAL(BDKi) :: Wrk(3) + + + !> ### Calculate the 1D strain, \f$ \underline{\epsilon} \f$, equation (5) + !! \f$ \underline{\epsilon} = \underline{x}^\prime_0 + \underline{u}^\prime - + !! \left(\underline{\underline{R}}\underline{\underline{R}}_0\right) \bar{\imath}_1 + !! = \underline{E}_1 - + !! \left(\underline{\underline{R}}\underline{\underline{R}}_0\right) \bar{\imath}_1 \f$ + !! where \f$ \bar{\imath}_1 \f$ is the unit vector along the \f$ x_1 \f$ direction in the inertial frame + !! (z in the IEC). The last term simplifies to the first column of \f$ \underline{\underline{R}}\underline{\underline{R}}_0 \f$ + !! + !! Note: \f$ \underline{\underline{R}}\underline{\underline{R}}_0 \f$ is used to go from the material basis into the inertial basis + !! and the transpose for the other direction. + eee(1:3) = m%qp%E1(1:3,idx_qp,nelem) - m%qp%RR0(1:3,3,idx_qp,nelem) ! Using RR0 z direction in IEC coords - IF(fact) THEN - !> ###Calculate the \f$ \underline{\underline{\mathcal{O}}} \f$ from equation (19) + + !> ### Set the 1D sectional curvature, \f$ \underline{\kappa} \f$, equation (5) + !! \f$ \underline{\kappa} = \underline{k} + \underline{\underline{R}}\underline{k}_i \f$ + !! where \f$ \underline{k} = \text{axial}\left(\underline{\underline{R}}^\prime\underline{\underline{R}}^T \right) \f$ + !! and \f$ \underline{k}_i \f$ is the corresponding initial curvature vector (root reference). + !! Note that \f$ \underline{k} \f$ can be calculated with rotation parameters as + !! \f$ \underline{k} = \underline{\underline{H}}(\underline{p}) \underline{p}' \f$ !! - !! \f$ \underline{\underline{\mathcal{O}}} = - !! \begin{bmatrix} - !! \underline{\underline{0}} & \underline{\underline{C}}_{11} \tilde{E}_1 - \tilde{F} \\ - !! \underline{\underline{0}} & \underline{\underline{C}}_{21} \tilde{E}_1 - \tilde{M} - !! \end{bmatrix} - !! = \begin{bmatrix} - !! \underline{\underline{0}} & \psi_E - \tilde{F} \\ - !! \underline{\underline{0}} & \mu - \tilde{M} - !! \end{bmatrix} - !! \f$ - Wrk33(:,:) = OuterProduct(m%qp%RR0(1:3,3,idx_qp,nelem), m%qp%RR0(1:3,3,idx_qp,nelem)) ! z-direction in IEC coords - C21(:,:) = m%qp%Stif(4:6,1:3,idx_qp,nelem) + cet*k1s*Wrk33(:,:) - - tildeE = SkewSymMat(m%qp%E1(:,idx_qp,nelem)) - epsi(:,:) = MATMUL(m%qp%Stif(1:3,1:3,idx_qp,nelem),tildeE) ! Stif is RR0 * p%Stif0_QP * RR0^T - mu(:,:) = MATMUL(C21,tildeE) - - m%qp%Oe(:,:,idx_qp,nelem) = 0.0_BDKi - m%qp%Oe(1:3,4:6,idx_qp,nelem) = epsi(1:3,1:3) - SkewSymMat(m%qp%Fc(1:3,idx_qp,nelem)) - m%qp%Oe(4:6,4:6,idx_qp,nelem) = mu(1:3,1:3) - SkewSymMat(m%qp%Fc(4:6,idx_qp,nelem)) - - - !> ###Calculated \f$ \underline{\underline{\mathcal{P}}} \f$ from equation (20) + !! \f$ \text{axial}\left( \underline{\underline{A}} \right) \f$ is defined as \n + !! \f$ \text{axial}\left( \underline{\underline{A}} \right) + !! = \left\{ \begin{array}{c} a_1 \\ + !! a_2 \\ + !! a_3 \end{array}\right\} + !! = \left\{ \begin{array}{c} A_{32} - A_{23} \\ + !! A_{13} - A_{31} \\ + !! A_{21} - A_{12} \end{array}\right\} + !! \f$ + !! In other words, \f$ \tilde{k} = \left(\underline{\underline{R}}^\prime\underline{\underline{R}}^T \right) \f$. + !! Note: \f$ \underline{\kappa} \f$ was already calculated in the BD_DisplacementQP routine + eee(4:6) = m%qp%kappa(1:3,idx_qp,nelem) + + + !FIXME: note that the k_i terms may not be documented correctly here. + !> ### Calculate the elastic force, \f$ \underline{\mathcal{F}}^C \f$ + !! Using equations (15), (23), (24), and (5) we write the elastic force as\n + !! \f$ + !! \underline{\mathcal{F}}^C + !! = \left\{ \begin{array}{c} \underline{F} \\ + !! \underline{M} \end{array} \right\} + !! = \left(\underline{\underline{R}}\underline{\underline{R}}_0\right) + !! \underline{\underline{C}} + !! \left(\underline{\underline{R}}\underline{\underline{R}}_0\right)^T + !! \left\{ \begin{array}{c} \underline{\epsilon} \\ + !! \underline{\kappa} \end{array} \right\} + !! \f$ !! - !! \f$ \underline{\underline{\mathcal{P}}} = - !! \begin{bmatrix} - !! \underline{\underline{0}} & \underline{\underline{0}} \\ - !! \left(\underline{\underline{\bar{C}}}_{11} \tilde{E}_1 \right)^T + \tilde{F} - !! \left(\underline{\underline{\bar{C}}}_{11} \tilde{E}_1 \right)^T - !! \end{bmatrix} - !! = \begin{bmatrix} - !! \underline{\underline{0}} & \underline{\underline{0}} \\ - !! \psi_E^T + \tilde{F} & \mu^T - !! \end{bmatrix} \f$ - m%qp%Pe(:,:,idx_qp,nelem) = 0.0_BDKi - m%qp%Pe(4:6,1:3,idx_qp,nelem) = TRANSPOSE(epsi) + SkewSymMat(m%qp%Fc(1:3,idx_qp,nelem)) - m%qp%Pe(4:6,4:6,idx_qp,nelem) = TRANSPOSE(mu) - - !> ###Calculated \f$ \underline{\underline{\mathcal{Q}}} \f$ from equation (21) + !! Note then that the extension twist term is added to this so that we get + !! \f$ \underline{\mathcal{F}}^C = \underline{F}^c_{a} + \underline{F}^c_{et} \f$ !! - !! \f{eqnarray*}{ - !! \underline{\underline{\mathcal{Q}}} - !! & =& \underline{\underline{\Upsilon}} \underline{\underline{\mathcal{O}}} - !! = \begin{bmatrix} 0 & 0 \\ - !! \tilde{E}_1^T & 0 \end{bmatrix} - !! \underline{\underline{\mathcal{O}}} \\ - !! \begin{bmatrix} 0 & 0 \\ - !! 0 & \underline{\underline{\mathcal{Q}}}_{22} \end{bmatrix} - !! & =& \tilde{E}_1^T \underline{\underline{\mathcal{O}}}_{12} - !! = - \tilde{E}_1 \underline{\underline{\mathcal{O}}}_{12} - !! \f}\n - !! Note: \f$ \tilde{E}_1^T = - \tilde{E}_1 \f$ - m%qp%Qe(:,:,idx_qp,nelem) = 0.0_BDKi - m%qp%Qe(4:6,4:6,idx_qp,nelem) = -MATMUL(tildeE,m%qp%Oe(1:3,4:6,idx_qp,nelem)) - ENDIF - + !! where \f$ \underline{F}^c_{et} \f$ is the extension twist coupling term. + + !> ###Calculate the first term. + !! \f$ \underline{F}^c_a + !! = \left(\underline{\underline{R}}\underline{\underline{R}}_0\right) + !! \underline{\underline{C}} + !! \left(\underline{\underline{R}}\underline{\underline{R}}_0\right)^T + !! \left\{ \begin{array}{c} + !! \underline{\epsilon} \\ + !! \underline{k} + !! \end{array} \right\} \f$ + !! + fff(1:6) = MATMUL(m%qp%Stif(:,:,idx_qp,nelem),eee) + + + !> ###Calculate the extension twist coupling. + !! This is calculated in the material basis, not in the rotated reference frame \f$ \underline{F}^c_a \f$ was calculated in. + !! First find \f$ \epsilon_1 \f$ and \f$ \kappa_1 \f$ in the material basis\n + !! \f$ \epsilon_{m} = \left( \underline{\underline{R}}\underline{\underline{R}}_0 \right) ^T \epsilon \f$ \n + !! and \n + !! \f$ \kappa_{m} = \left( \underline{\underline{R}}\underline{\underline{R}}_0 \right) ^T \underline{k}\f$ \n + + ! Strain into the material basis (eq (39) of Dymore manual) + !Wrk(:) = MATMUL(TRANSPOSE(m%qp%RR0(:,:,idx_qp,nelem)),eee(1:3)) + !e1s = Wrk(3) !epsilon_{1} in material basis (for major axis of blade, which is z in the IEC formulation) + e1s = dot_product( m%qp%RR0(:,3,idx_qp,nelem), eee(1:3) ) + + !Wrk(:) = MATMUL(TRANSPOSE(m%qp%RR0(:,:,idx_qp,nelem)),eee(4:6)) + !k1s = Wrk(3) !kappa_{1} in material basis (for major axis of blade, which is z in the IEC formulation) + k1s = dot_product( m%qp%RR0(:,3,idx_qp,nelem), eee(4:6) ) + + + !> Add extension twist coupling terms to the \f$ \underline{F}^c_{a} \f$\n + !! \f$ \underline{F}^c = + !! \underline{F}^c_{a} + \underline{F}^c_{et} \f$\n + !! The extension twist term is calculated in the material basis, which can be simplified and rewritten as follows (_the details of this are fuzzy to me_):\n + !! \f$ \underline{F}^c_{et} = + !! \begin{bmatrix} + !! \frac{1}{2} C_{et} \kappa^2_{m}(1) \left(\underline{\underline{R}}\underline{\underline{R}}_0\right)_{:,1} \\ + !! C_{et} \kappa_{m}(1) \epsilon_{m}(1) \left(\underline{\underline{R}}\underline{\underline{R}}_0\right)_{:,1} + !! \end{bmatrix} + !! \f$ + !! + !! Where \f$ C_{et} = C_{5,5} + C_{6,6} \f$ is the inertial term along major axis, \f$x\f$, in the material coordinate system + !! Note that with coverting to the FAST / IEC coordinate system, we now are using the Ixx and Iyy terms which are located at + !! \f$ C_{et} = C_{4,4} + C_{5,5} \f$ + ! Refer Section 1.4 in "Dymore User's Manual - Formulation and finite element implementation of beam elements". + cet= p%Stif0_QP(4,4,(nelem-1)*p%nqp+idx_qp) + p%Stif0_QP(5,5,(nelem-1)*p%nqp+idx_qp) ! Dymore theory (22) + m%qp%Fc(1:3,idx_qp,nelem) = fff(1:3) + 0.5_BDKi*cet*k1s*k1s*m%qp%RR0(1:3,3,idx_qp,nelem) ! Dymore theory (25a). Note z-axis is the length of blade. + m%qp%Fc(4:6,idx_qp,nelem) = fff(4:6) + cet*e1s*k1s*m%qp%RR0(1:3,3,idx_qp,nelem) ! Dymore theory (25b). Note z-axis is the length of blade. + + !> ###Calculate \f$ \underline{\mathcal{F}}^d \f$, equation (16) + !! \f$ \underline{F}^d = + !! \underline{\underline{\Upsilon}} \underline{\mathcal{F}}^c + !! = \begin{bmatrix} \underline{\underline{0}} & \underline{\underline{0}} \\ + !! \tilde{E}_1^T & \underline{\underline{0}} \end{bmatrix} + !! \underline{\mathcal{F}}^c + !! = \begin{bmatrix} \underline{0} \\ + !! \left(\underline{\mathcal{F}}^c \times \underline{E}_1 \right)^T + !! \end{bmatrix} \f$ + m%qp%Fd(1:3,idx_qp,nelem) = 0.0_BDKi + ! ADP uu0 ref: If E1 is referenced against a different curve than Stif0_QP, there will be strange coupling terms here. + m%qp%Fd(4:6,idx_qp,nelem) = cross_product(m%qp%Fc(1:3,idx_qp,nelem), m%qp%E1(:,idx_qp,nelem)) + + end subroutine Calc_Fc_Fd END SUBROUTINE BD_ElasticForce @@ -2354,27 +2420,36 @@ END SUBROUTINE BD_QPDataVelocity ! Output: m%qp%aaa ! !NOTE: This routine used to be part of BD_QPDataVelocity -SUBROUTINE BD_QPDataAcceleration( nelem, idx_qp, p, OtherState, m ) +SUBROUTINE BD_QPDataAcceleration( p, OtherState, m ) - INTEGER(IntKi), INTENT(IN ) :: nelem !< index of current element - INTEGER(IntKi), INTENT(IN ) :: idx_qp !< index of current quadrature point TYPE(BD_ParameterType), INTENT(IN ) :: p !< Parameters TYPE(BD_OtherStateType), INTENT(IN ) :: OtherState !< Other states at t on input; at t+dt on outputs TYPE(BD_MiscVarType), INTENT(INOUT) :: m !< Misc/optimization variables + INTEGER(IntKi) :: nelem !< index of current element + INTEGER(IntKi) :: idx_qp !< index of current quadrature point INTEGER(IntKi) :: idx_node INTEGER(IntKi) :: elem_start - elem_start = p%node_elem_idx(nelem,1) ! Initialize to zero for summation - m%qp%aaa(:,idx_qp,nelem) = 0.0_BDKi + m%qp%aaa = 0.0_BDKi ! Calculate the and acceleration term at t+dt (OtherState%acc is at t+dt) - DO idx_node=1,p%nodes_per_elem - m%qp%aaa(:,idx_qp,nelem) = m%qp%aaa(:,idx_qp,nelem) + p%Shp(idx_node,idx_qp) * OtherState%acc(:,elem_start-1+idx_node) - ENDDO + + DO nelem=1,p%elem_total + + elem_start = p%node_elem_idx(nelem,1) + + DO idx_qp=1,p%nqp + DO idx_node=1,p%nodes_per_elem + m%qp%aaa(:,idx_qp,nelem) = m%qp%aaa(:,idx_qp,nelem) + p%Shp(idx_node,idx_qp) * OtherState%acc(:,elem_start-1+idx_node) + END DO + END DO + + END DO + END SUBROUTINE BD_QPDataAcceleration @@ -2382,10 +2457,9 @@ END SUBROUTINE BD_QPDataAcceleration !----------------------------------------------------------------------------------------------------------------------------------- !> This subroutine calculates the inertial force `m%qp%Fi` !! It also calcuates the linearized matrices `m%qp%Mi`, `m%qp%Gi`, and `m%qp%Ki` for N-R algorithm -SUBROUTINE BD_InertialForce( nelem, idx_qp, p, m, fact ) +SUBROUTINE BD_InertialForce( nelem, p, m, fact ) INTEGER(IntKi), INTENT(IN ) :: nelem !< index of current element - INTEGER(IntKi), INTENT(IN ) :: idx_qp !< index of current quadrature point TYPE(BD_ParameterType), INTENT(IN ) :: p !< Parameters TYPE(BD_MiscVarType), INTENT(INOUT) :: m !< Misc/optimization variables. LOGICAL, INTENT(IN ) :: fact @@ -2395,38 +2469,59 @@ SUBROUTINE BD_InertialForce( nelem, idx_qp, p, m, fact ) REAL(BDKi) :: nu(3) REAL(BDKi) :: epsi(3,3) REAL(BDKi) :: mu(3,3) - - - beta = cross_product(m%qp%vvv(4:6,idx_qp,nelem), m%qp%RR0mEta(:,idx_qp,nelem)) !MATMUL(SkewSymMat(ome),mEta) - gama = MATMUL(m%qp%rho(:,:,idx_qp,nelem),m%qp%vvv(4:6,idx_qp,nelem)) - nu = MATMUL(m%qp%rho(:,:,idx_qp,nelem),m%qp%aaa(4:6,idx_qp,nelem)) - - !Compute Fi (Equation 22 in Wang_2014) - m%qp%Fi(1:3,idx_qp,nelem)= p%qp%mmm(idx_qp,nelem)*m%qp%aaa(1:3,idx_qp,nelem) + cross_product(m%qp%aaa(4:6,idx_qp,nelem), m%qp%RR0mEta(:,idx_qp,nelem)) + cross_product(m%qp%vvv(4:6,idx_qp,nelem),beta) ! MATMUL(SkewSymMat(aaa(4:6)),mEta)+MATMUL(SkewSymMat(ome),beta) - m%qp%Fi(4:6,idx_qp,nelem) = cross_product(m%qp%RR0mEta(:,idx_qp,nelem), m%qp%aaa(1:3,idx_qp,nelem)) + nu + cross_product(m%qp%vvv(4:6,idx_qp,nelem), gama) ! MATMUL(SkewSymMat(mEta),aaa(1:3)) + nu + MATMUL(SkewSymMat(ome),gama) - + + REAL(BDKi) :: SSmat_vvv(3,3) + REAL(BDKi) :: SSmat_RR0mEta(3,3) + + INTEGER(IntKi) :: idx_qp !< index of current quadrature point + + ! Note, in the following there are cases where SkewSymMat(A)*SkewSymMat(B)*C has been reduced to + ! cross_product(A, cross_product(B, C)). This is correct because SkewSymMat(B)*C = cross_product(B, C) + ! but SkewSymMat(A)*SkewSymMat(B) does not give a vector so the cross products cannot be switched. + ! Note that in the matrix multiplication (A*B)*C = A*(B*C) + + ! Note that F_i is calculated in the inertial frame, not the rotating frame. + ! Note that RR0*mEta is the rotation at the current node times the (m*cm_x,m*cm_y,0) of the blade section at qp. + ! So, the omega terms are therefore correction terms? + do idx_qp=1,p%nqp + beta = cross_product(m%qp%vvv(4:6,idx_qp,nelem), m%qp%RR0mEta(:,idx_qp,nelem)) !MATMUL(SkewSymMat(ome),mEta) + gama = MATMUL(m%qp%rho(:,:,idx_qp,nelem),m%qp%vvv(4:6,idx_qp,nelem)) + nu = MATMUL(m%qp%rho(:,:,idx_qp,nelem),m%qp%aaa(4:6,idx_qp,nelem)) + + !Compute Fi (Equation 22 in Wang_2014, Equation 17.110 in Bauchau's book) + m%qp%Fi(1:3,idx_qp,nelem) = p%qp%mmm(idx_qp,nelem)*m%qp%aaa(1:3,idx_qp,nelem) & + + cross_product(m%qp%aaa(4:6,idx_qp,nelem), m%qp%RR0mEta(:,idx_qp,nelem)) & + + cross_product(m%qp%vvv(4:6,idx_qp,nelem), beta) ! MATMUL(SkewSymMat(aaa(4:6)),mEta)+MATMUL(SkewSymMat(ome),beta) + m%qp%Fi(4:6,idx_qp,nelem) = cross_product(m%qp%RR0mEta(:,idx_qp,nelem), m%qp%aaa(1:3,idx_qp,nelem)) + nu + cross_product(m%qp%vvv(4:6,idx_qp,nelem), gama) ! MATMUL(SkewSymMat(mEta),aaa(1:3)) + nu + MATMUL(SkewSymMat(ome),gama) + end do + IF(fact) THEN - CALL BD_InertialMassMatrix( nelem, idx_qp, p, m ) - - !Gyroscopic Matrix (Equation 17 in Wang_2014) - m%qp%Gi(:,1:3,idx_qp,nelem) = 0.0_BDKi - epsi = MATMUL(SkewSymMat(m%qp%vvv(4:6,idx_qp,nelem)),m%qp%rho(:,:,idx_qp,nelem)) - mu = -MATMUL(SkewSymMat(m%qp%vvv(4:6,idx_qp,nelem)),SkewSymMat(m%qp%RR0mEta(:,idx_qp,nelem))) - m%qp%Gi(1:3,4:6,idx_qp,nelem) = -SkewSymMat(beta) + mu - m%qp%Gi(4:6,4:6,idx_qp,nelem) = epsi - SkewSymMat(gama) - - !Stiffness Matrix (Equation 18 in Wang_2014) - m%qp%Ki(:,1:3,idx_qp,nelem) = 0.0_BDKi - m%qp%Ki(1:3,4:6,idx_qp,nelem) = -MATMUL(SkewSymMat(m%qp%aaa(4:6,idx_qp,nelem)),SkewSymMat(m%qp%RR0mEta(:,idx_qp,nelem))) + & - MATMUL(SkewSymMat(m%qp%vvv(4:6,idx_qp,nelem)),mu) - m%qp%Ki(4:6,4:6,idx_qp,nelem) = MATMUL(SkewSymMat(m%qp%aaa(1:3,idx_qp,nelem)),SkewSymMat(m%qp%RR0mEta(:,idx_qp,nelem))) + & - MATMUL(m%qp%rho(:,:,idx_qp,nelem),SkewSymMat(m%qp%aaa(4:6,idx_qp,nelem))) - SkewSymMat(nu) + & - MATMUL(epsi,SkewSymMat(m%qp%vvv(4:6,idx_qp,nelem))) - & - MATMUL(SkewSymMat(m%qp%vvv(4:6,idx_qp,nelem)),SkewSymMat(gama)) - ELSE - m%qp%Mi(:,:,idx_qp,nelem) = 0.0_BDKi - m%qp%Gi(:,:,idx_qp,nelem) = 0.0_BDKi - m%qp%Ki(:,:,idx_qp,nelem) = 0.0_BDKi + CALL BD_InertialMassMatrix( nelem, p, m ) ! compute Mi + do idx_qp=1,p%nqp + + !Gyroscopic Matrix (Equation 17 in Wang_2014) + m%qp%Gi(:,1:3,idx_qp,nelem) = 0.0_BDKi + SSmat_vvv = SkewSymMat(m%qp%vvv(4:6,idx_qp,nelem)) + SSmat_RR0mEta = SkewSymMat(m%qp%RR0mEta(:,idx_qp,nelem)) + + epsi = MATMUL(SSmat_vvv,m%qp%rho(:,:,idx_qp,nelem)) + mu = -MATMUL(SSmat_vvv,SSmat_RR0mEta) + nu = MATMUL(m%qp%rho(:,:,idx_qp,nelem),m%qp%aaa(4:6,idx_qp,nelem)) + beta = cross_product(m%qp%vvv(4:6,idx_qp,nelem), m%qp%RR0mEta(:,idx_qp,nelem)) !MATMUL(SkewSymMat(ome),mEta) + gama = MATMUL(m%qp%rho(:,:,idx_qp,nelem),m%qp%vvv(4:6,idx_qp,nelem)) + m%qp%Gi(1:3,4:6,idx_qp,nelem) = -SkewSymMat(beta) + mu + m%qp%Gi(4:6,4:6,idx_qp,nelem) = epsi - SkewSymMat(gama) + + !Stiffness Matrix (Equation 18 in Wang_2014) + m%qp%Ki( : ,1:3,idx_qp,nelem) = 0.0_BDKi + m%qp%Ki(1:3,4:6,idx_qp,nelem) = -MATMUL(SkewSymMat(m%qp%aaa(4:6,idx_qp,nelem)),SSmat_RR0mEta) & + + MATMUL(SSmat_vvv,mu) + m%qp%Ki(4:6,4:6,idx_qp,nelem) = MATMUL(SkewSymMat(m%qp%aaa(1:3,idx_qp,nelem)),SSmat_RR0mEta) & + + MATMUL(m%qp%rho(:,:,idx_qp,nelem),SkewSymMat(m%qp%aaa(4:6,idx_qp,nelem))) & + - SkewSymMat(nu) & + + MATMUL(epsi,SSmat_vvv) & + - MATMUL(SSmat_vvv,SkewSymMat(gama)) + end do ENDIF END SUBROUTINE BD_InertialForce @@ -2436,15 +2531,13 @@ END SUBROUTINE BD_InertialForce !> This subroutine calculates the dissipative forces and added it to Fc and Fd !! It also calcuates the linearized matrices Sd, Od, Pd and Qd !! betaC, Gd, Xd, Yd for N-R algorithm -SUBROUTINE BD_DissipativeForce( nelem, idx_qp, p,m,fact ) - INTEGER(IntKi), INTENT(IN ) :: idx_qp !< index of current quadrature point +SUBROUTINE BD_DissipativeForce( nelem, p, m,fact ) INTEGER(IntKi), INTENT(IN ) :: nelem !< index of current element in loop TYPE(BD_ParameterType), INTENT(IN ) :: p !< Parameters TYPE(BD_MiscVarType), INTENT(INOUT) :: m !< Misc/optimization variables LOGICAL, INTENT(IN ) :: fact REAL(BDKi) :: SS_ome(3,3) - REAL(BDKi) :: eed(6) REAL(BDKi) :: ffd(6) REAL(BDKi) :: D11(3,3) REAL(BDKi) :: D12(3,3) @@ -2453,98 +2546,122 @@ SUBROUTINE BD_DissipativeForce( nelem, idx_qp, p,m,fact ) REAL(BDKi) :: b11(3,3) REAL(BDKi) :: b12(3,3) REAL(BDKi) :: alpha(3,3) - REAL(BDKi) :: temp_b(6,6) - INTEGER(IntKi) :: i + INTEGER(IntKi) :: i, j - m%qp%betaC(:,:,idx_qp,nelem) = 0.0_BDKi + INTEGER(IntKi) :: idx_qp !< index of current quadrature point + + DO idx_qp=1,p%nqp + !m%qp%betaC(:,:,idx_qp,nelem) = MATMUL( diag(p%beta(i)), temp_b,m%qp%Stif(:,:,idx_qp,nelem)) + DO j=1,6 + DO i=1,6 + m%qp%betaC(i,j,idx_qp,nelem) = p%beta(i)*m%qp%Stif(i,j,idx_qp,nelem) + END DO + END DO + END DO - ! Compute strain rates - eed(1:6) = m%qp%vvp(1:6,idx_qp,nelem) - eed(1:3) = eed(1:3) + cross_product(m%qp%E1(:,idx_qp,nelem),m%qp%vvv(4:6,idx_qp,nelem)) + + IF (.NOT. fact) then ! skip all but Fc and Fd terms + + DO idx_qp=1,p%nqp + call Calc_FC_FD_ffd() ! this modifies m%qp%Fc and m%qp%Fd + END DO + + ! bjj: we don't use these values when fact is FALSE, so let's save time and ignore them here, too. + ! m%qp%Sd(:,:,:,nelem) = 0.0_BDKi + ! m%qp%Pd(:,:,:,nelem) = 0.0_BDKi + ! m%qp%Od(:,:,:,nelem) = 0.0_BDKi + ! m%qp%Qd(:,:,:,nelem) = 0.0_BDKi + ! m%qp%Gd(:,:,:,nelem) = 0.0_BDKi + ! m%qp%Xd(:,:,:,nelem) = 0.0_BDKi + ! m%qp%Yd(:,:,:,nelem) = 0.0_BDKi + + ELSE +!FIXME: sometime we can condense this with vector arithmetic and removing some variables that aren't needed. + + DO idx_qp=1,p%nqp - ! Compute damping matrix - temp_b(:,:) = 0.0_BDKi - ! set diagonal terms - DO i=1,6 - temp_b(i,i) = p%beta(i) - ENDDO - m%qp%betaC(1:6,1:6,idx_qp,nelem) = MATMUL(temp_b,m%qp%Stif(:,:,idx_qp,nelem)) - D11(1:3,1:3) = m%qp%betaC(1:3,1:3,idx_qp,nelem) - D12(1:3,1:3) = m%qp%betaC(1:3,4:6,idx_qp,nelem) - D21(1:3,1:3) = m%qp%betaC(4:6,1:3,idx_qp,nelem) - D22(1:3,1:3) = m%qp%betaC(4:6,4:6,idx_qp,nelem) + CALL Calc_FC_FD_ffd() ! this sets local variable ffd and modifies m%qp%Fc and m%qp%Fd + + D11 = m%qp%betaC(1:3,1:3,idx_qp,nelem) + D12 = m%qp%betaC(1:3,4:6,idx_qp,nelem) + D21 = m%qp%betaC(4:6,1:3,idx_qp,nelem) + D22 = m%qp%betaC(4:6,4:6,idx_qp,nelem) + + b11(1:3,1:3) = -MATMUL(SkewSymMat(m%qp%E1(:,idx_qp,nelem)),D11) + b12(1:3,1:3) = -MATMUL(SkewSymMat(m%qp%E1(:,idx_qp,nelem)),D12) + + SS_ome = SkewSymMat( m%qp%vvv(4:6,idx_qp,nelem) ) + + ! Compute stiffness matrix Sd + m%qp%Sd(1:3,1:3,idx_qp,nelem) = -MATMUL(D11,SS_ome) + m%qp%Sd(1:3,4:6,idx_qp,nelem) = -MATMUL(D12,SS_ome) + m%qp%Sd(4:6,1:3,idx_qp,nelem) = -MATMUL(D21,SS_ome) + m%qp%Sd(4:6,4:6,idx_qp,nelem) = -MATMUL(D22,SS_ome) + + ! Compute stiffness matrix Pd + m%qp%Pd(:,:,idx_qp,nelem) = 0.0_BDKi + m%qp%Pd(4:6,1:3,idx_qp,nelem) = SkewSymMat(ffd(1:3)) - MATMUL(b11,SS_ome) + m%qp%Pd(4:6,4:6,idx_qp,nelem) = -MATMUL(b12,SS_ome) + + ! Compute stiffness matrix Od + m%qp%Od(:,1:3,idx_qp,nelem) = 0.0_BDKi + alpha = SkewSymMat(m%qp%vvp(1:3,idx_qp,nelem)) - MATMUL(SS_ome,SkewSymMat(m%qp%E1(:,idx_qp,nelem))) + m%qp%Od(1:3,4:6,idx_qp,nelem) = MATMUL(D11,alpha) - SkewSymMat(ffd(1:3)) + m%qp%Od(4:6,4:6,idx_qp,nelem) = MATMUL(D21,alpha) - SkewSymMat(ffd(4:6)) + + ! Compute stiffness matrix Qd + m%qp%Qd(:,:,idx_qp,nelem) = 0.0_BDKi + m%qp%Qd(4:6,4:6,idx_qp,nelem) = -MATMUL(SkewSymMat(m%qp%E1(:,idx_qp,nelem)),m%qp%Od(1:3,4:6,idx_qp,nelem)) + ! Compute gyroscopic matrix Gd + m%qp%Gd(:,1:3,idx_qp,nelem) = 0.0_BDKi + m%qp%Gd(1:3,4:6,idx_qp,nelem) = TRANSPOSE(b11) + m%qp%Gd(4:6,4:6,idx_qp,nelem) = TRANSPOSE(b12) + + ! Compute gyroscopic matrix Xd + m%qp%Xd(:,:,idx_qp,nelem) = 0.0_BDKi + m%qp%Xd(4:6,4:6,idx_qp,nelem) = -MATMUL(SkewSymMat(m%qp%E1(:,idx_qp,nelem)),m%qp%Gd(1:3,4:6,idx_qp,nelem)) + + ! Compute gyroscopic matrix Yd + m%qp%Yd(1:3,:,idx_qp,nelem) = 0.0_BDKi + m%qp%Yd(4:6,1:3,idx_qp,nelem) = b11 + m%qp%Yd(4:6,4:6,idx_qp,nelem) = b12 + END DO + ENDIF - ! Compute dissipative force - ffd(1:6) = MATMUL(m%qp%betaC(:,:,idx_qp,nelem),eed) - m%qp%Fc(1:6,idx_qp,nelem) = m%qp%Fc(1:6,idx_qp,nelem) + ffd(1:6) - m%qp%Fd(4:6,idx_qp,nelem) = m%qp%Fd(4:6,idx_qp,nelem) + cross_product(ffd(1:3),m%qp%E1(:,idx_qp,nelem)) +CONTAINS + SUBROUTINE Calc_FC_FD_ffd() + REAL(BDKi) :: eed(6) + + ! Compute strain rates + eed = m%qp%vvp(1:6,idx_qp,nelem) + eed(1:3) = eed(1:3) + cross_product(m%qp%E1(:,idx_qp,nelem),m%qp%vvv(4:6,idx_qp,nelem)) - IF(fact) THEN - !ome(1:3) = vvv(4:6) - SS_ome = SkewSymMat( m%qp%vvv(4:6,idx_qp,nelem) ) - - ! Compute stiffness matrix Sd - m%qp%Sd(1:3,1:3,idx_qp,nelem) = -MATMUL(D11,SS_ome) - m%qp%Sd(1:3,4:6,idx_qp,nelem) = -MATMUL(D12,SS_ome) - m%qp%Sd(4:6,1:3,idx_qp,nelem) = -MATMUL(D21,SS_ome) - m%qp%Sd(4:6,4:6,idx_qp,nelem) = -MATMUL(D22,SS_ome) - - ! Compute stiffness matrix Pd - m%qp%Pd(:,:,idx_qp,nelem) = 0.0_BDKi - b11(1:3,1:3) = -MATMUL(SkewSymMat(m%qp%E1(:,idx_qp,nelem)),D11) - b12(1:3,1:3) = -MATMUL(SkewSymMat(m%qp%E1(:,idx_qp,nelem)),D12) - m%qp%Pd(4:6,1:3,idx_qp,nelem) = SkewSymMat(ffd(1:3)) - MATMUL(b11,SS_ome) - m%qp%Pd(4:6,4:6,idx_qp,nelem) = -MATMUL(b12,SS_ome) !bjj: changed this index based on how I read Eq 20 of Qi's paper. Qi confirmed this was a bug. - - ! Compute stiffness matrix Od - m%qp%Od(:,1:3,idx_qp,nelem) = 0.0_BDKi - alpha(1:3,1:3) = SkewSymMat(m%qp%vvp(1:3,idx_qp,nelem)) - MATMUL(SS_ome,SkewSymMat(m%qp%E1(:,idx_qp,nelem))) - m%qp%Od(1:3,4:6,idx_qp,nelem) = MATMUL(D11,alpha) - SkewSymMat(ffd(1:3)) - m%qp%Od(4:6,4:6,idx_qp,nelem) = MATMUL(D21,alpha) - SkewSymMat(ffd(4:6)) - - ! Compute stiffness matrix Qd - m%qp%Qd(:,:,idx_qp,nelem) = 0.0_BDKi - m%qp%Qd(4:6,4:6,idx_qp,nelem) = -MATMUL(SkewSymMat(m%qp%E1(:,idx_qp,nelem)),m%qp%Od(1:3,4:6,idx_qp,nelem)) - - ! Compute gyroscopic matrix Gd - m%qp%Gd(:,1:3,idx_qp,nelem) = 0.0_BDKi - m%qp%Gd(1:3,4:6,idx_qp,nelem) = TRANSPOSE(b11) - m%qp%Gd(4:6,4:6,idx_qp,nelem) = TRANSPOSE(b12) - - ! Compute gyroscopic matrix Xd - m%qp%Xd(:,:,idx_qp,nelem) = 0.0_BDKi - m%qp%Xd(4:6,4:6,idx_qp,nelem) = -MATMUL(SkewSymMat(m%qp%E1(:,idx_qp,nelem)),m%qp%Gd(1:3,4:6,idx_qp,nelem)) - - ! Compute gyroscopic matrix Yd - m%qp%Yd(1:3,:,idx_qp,nelem) = 0.0_BDKi - m%qp%Yd(4:6,1:3,idx_qp,nelem) = b11 - m%qp%Yd(4:6,4:6,idx_qp,nelem) = b12 - ELSE - m%qp%Sd(:,:,idx_qp,nelem) = 0.0_BDKi - m%qp%Pd(:,:,idx_qp,nelem) = 0.0_BDKi - m%qp%Od(:,:,idx_qp,nelem) = 0.0_BDKi - m%qp%Qd(:,:,idx_qp,nelem) = 0.0_BDKi - m%qp%Gd(:,:,idx_qp,nelem) = 0.0_BDKi - m%qp%Xd(:,:,idx_qp,nelem) = 0.0_BDKi - m%qp%Yd(:,:,idx_qp,nelem) = 0.0_BDKi - ENDIF + ! Compute dissipative force + ffd(1:6) = MATMUL(m%qp%betaC(:,:,idx_qp,nelem),eed) + m%qp%Fc(1:6,idx_qp,nelem) = m%qp%Fc(1:6,idx_qp,nelem) + ffd + m%qp%Fd(4:6,idx_qp,nelem) = m%qp%Fd(4:6,idx_qp,nelem) + cross_product(ffd(1:3),m%qp%E1(:,idx_qp,nelem)) + + END SUBROUTINE Calc_FC_FD_ffd END SUBROUTINE BD_DissipativeForce !----------------------------------------------------------------------------------------------------------------------------------- !> This subroutine calculates the gravity forces `m%qp%Fg` -SUBROUTINE BD_GravityForce( nelem,idx_qp,p,m,grav ) +SUBROUTINE BD_GravityForce( nelem,p,m,grav ) INTEGER(IntKi), INTENT(IN ) :: nelem !< index of current element in loop - INTEGER(IntKi), INTENT(IN ) :: idx_qp !< index of current quadrature point TYPE(BD_ParameterType), INTENT(IN ) :: p !< Parameters TYPE(BD_MiscVarType), INTENT(INOUT) :: m !< Misc/optimization variables. REAL(BDKi), INTENT(IN ) :: grav(:) !< Gravity, which is scaled in the case of Static analysis + INTEGER(IntKi) :: idx_qp !< index of current quadrature point + do idx_qp=1,p%nqp - m%qp%Fg(1:3,idx_qp,nelem) = p%qp%mmm(idx_qp,nelem) * grav(1:3) - m%qp%Fg(4:6,idx_qp,nelem) = MATMUL(SkewSymMat(m%qp%RR0mEta(:,idx_qp,nelem)),grav) + m%qp%Fg(1:3,idx_qp,nelem) = p%qp%mmm(idx_qp,nelem) * grav(1:3) + m%qp%Fg(4:6,idx_qp,nelem) = cross_product(m%qp%RR0mEta(:,idx_qp,nelem),grav) + end do + END SUBROUTINE BD_GravityForce @@ -2558,7 +2675,6 @@ SUBROUTINE BD_AssembleStiffK(nelem,p,ElemK,GlobalK) INTEGER(IntKi) :: i INTEGER(IntKi) :: j - INTEGER(IntKi) :: idx_dof1 INTEGER(IntKi) :: idx_dof2 INTEGER(IntKi) :: temp_id @@ -2566,9 +2682,7 @@ SUBROUTINE BD_AssembleStiffK(nelem,p,ElemK,GlobalK) DO j=1,p%nodes_per_elem DO idx_dof2=1,p%dof_node DO i=1,p%nodes_per_elem - DO idx_dof1=1,p%dof_node - GlobalK( idx_dof1,i+temp_id,idx_dof2,j+temp_id ) = GlobalK( idx_dof1,i+temp_id,idx_dof2,j+temp_id ) + ElemK( idx_dof1,i,idx_dof2,j ) - ENDDO + GlobalK( :,i+temp_id,idx_dof2,j+temp_id ) = GlobalK( :,i+temp_id,idx_dof2,j+temp_id ) + ElemK( :,i,idx_dof2,j ) ENDDO ENDDO ENDDO @@ -2614,46 +2728,58 @@ SUBROUTINE BD_ElementMatrixAcc( nelem, p, x, m ) INTEGER(IntKi) :: idx_qp INTEGER(IntKi) :: i INTEGER(IntKi) :: j - INTEGER(IntKi) :: idx_dof1 - INTEGER(IntKi) :: idx_dof2 + INTEGER(IntKi) :: idx_dof1, idx_dof2 CHARACTER(*), PARAMETER :: RoutineName = 'BD_ElementMatrixAcc' - m%elf(:,:,:) = 0.0_BDKi - m%elm(:,:,:,:) = 0.0_BDKi - - DO idx_qp=1,p%nqp - + CALL BD_ElasticForce( nelem, p, m, .FALSE. ) ! Calculate Fc, Fd only + IF(p%damp_flag .NE. 0) THEN + CALL BD_DissipativeForce( nelem, p, m, .FALSE. ) ! Calculate dissipative terms on Fc, Fd + ENDIF + CALL BD_GravityForce( nelem, p, m, p%gravity ) ! Calculate Fg + CALL BD_GyroForce( nelem, p, m ) ! Calculate Fb (velocity terms from InertialForce with aaa=0) - CALL BD_InertialMassMatrix( nelem, idx_qp, p, m ) - CALL BD_GyroForce( nelem, idx_qp, m ) - CALL BD_GravityForce( nelem, idx_qp, p, m, p%gravity ) - CALL BD_ElasticForce( nelem,idx_qp,p,m,.FALSE. ) ! Calculate Fc, Fd only - IF(p%damp_flag .NE. 0) THEN - CALL BD_DissipativeForce( nelem, idx_qp, p,m,.FALSE. ) ! Calculate dissipative terms on Fc, Fd - - ENDIF - - m%qp%Ftemp(:,idx_qp,nelem) = m%qp%Fd(:,idx_qp,nelem) + m%qp%Fb(:,idx_qp,nelem) - m%DistrLoad_QP(:,idx_qp,nelem) - m%qp%Fg(:,idx_qp,nelem) + + m%qp%Ftemp(:,:,nelem) = m%qp%Fd(:,:,nelem) + m%qp%Fb(:,:,nelem) - m%DistrLoad_QP(:,:,nelem) - m%qp%Fg(:,:,nelem) - DO j=1,p%nodes_per_elem - DO idx_dof2=1,p%dof_node - DO i=1,p%nodes_per_elem - DO idx_dof1=1,p%dof_node - m%elm(idx_dof1,i,idx_dof2,j) = m%elm(idx_dof1,i,idx_dof2,j) + p%Shp(i,idx_qp)*m%qp%Mi(idx_dof1,idx_dof2,idx_qp,nelem)*p%Shp(j,idx_qp)*p%Jacobian(idx_qp,nelem)*p%QPtWeight(idx_qp) - ENDDO + + CALL BD_InertialMassMatrix( nelem, p, m ) ! Calculate Mi + + DO j=1,p%nodes_per_elem + DO idx_dof2=1,p%dof_node + DO i=1,p%nodes_per_elem + DO idx_dof1=1,p%dof_node + m%elm(idx_dof1,i,idx_dof2,j) = 0.0_BDKi + DO idx_qp = 1,p%nqp + m%elm(idx_dof1,i,idx_dof2,j) = m%elm(idx_dof1,i,idx_dof2,j) + m%qp%Mi(idx_dof1,idx_dof2,idx_qp,nelem)*p%QPtw_Shp_Shp_Jac(idx_qp,i,j,nelem) + END DO ENDDO ENDDO ENDDO - DO i=1,p%nodes_per_elem - DO idx_dof1=1,p%dof_node - m%elf(idx_dof1,i,nelem) = m%elf(idx_dof1,i,nelem) - p%ShpDer(i,idx_qp)*m%qp%Fc (idx_dof1,idx_qp,nelem)*p%QPtWeight(idx_qp) - m%elf(idx_dof1,i,nelem) = m%elf(idx_dof1,i,nelem) - p%Shp(i,idx_qp) *m%qp%Ftemp(idx_dof1,idx_qp,nelem)*p%Jacobian(idx_qp,nelem)*p%QPtWeight(idx_qp) - ENDDO - ENDDO - + end do + + DO i=1,p%nodes_per_elem + DO idx_dof1=1,p%dof_node + + m%elf(idx_dof1,i) = 0.0_BDKi + DO idx_qp = 1,p%nqp ! dot_product(m%qp%Fc(idx_dof1,:,nelem),p%QPtw_ShpDer(:,i)) + m%elf(idx_dof1,i) = m%elf(idx_dof1,i) - m%qp%Fc(idx_dof1,idx_qp,nelem)*p%QPtw_ShpDer(idx_qp,i) + END DO + + DO idx_qp = 1,p%nqp ! dot_product(m%qp%Ftemp(idx_dof1,:,nelem), p%QPtw_Shp_Jac(:,i,nelem)) + m%elf(idx_dof1,i) = m%elf(idx_dof1,i) - m%qp%Ftemp(idx_dof1,idx_qp,nelem)*p%QPtw_Shp_Jac(idx_qp,i,nelem) + END DO + ! Elastic force internal + m%EFint(idx_dof1,i,nelem) = 0.0_BDKi + DO idx_qp = 1,p%nqp ! dot_product(m%qp%Fc(idx_dof1,:,nelem), p%QPtw_ShpDer( :,i)) + m%EFint(idx_dof1,i,nelem) = m%EFint(idx_dof1,i,nelem) - m%qp%Fc(idx_dof1,idx_qp,nelem)*p%QPtw_ShpDer(idx_qp,i) + END DO +! m%EFint(idx_dof1,i,nelem) = m%EFint(idx_dof1,i,nelem) - dot_product(m%qp%Fd(idx_dof1,:,nelem),p%QPtw_Shp_Jac(:,i,nelem)) + + ENDDO ENDDO + RETURN END SUBROUTINE BD_ElementMatrixAcc @@ -2661,62 +2787,64 @@ END SUBROUTINE BD_ElementMatrixAcc !> This subroutine computes the mass matrix. !! Returns value for `m%qp%Mi` ! used in BD_ElementMatrixAcc and BD_InertialForce -SUBROUTINE BD_InertialMassMatrix( nelem, idx_qp, p, m ) +SUBROUTINE BD_InertialMassMatrix( nelem, p, m ) INTEGER(IntKi), INTENT(IN ) :: nelem !< index of current element in loop - INTEGER(IntKi), INTENT(IN ) :: idx_qp !< index of current quadrature point TYPE(BD_ParameterType), INTENT(IN ) :: p !< Parameters TYPE(BD_MiscVarType), INTENT(INOUT) :: m !< Misc/optimization variables INTEGER(IntKi) :: i + INTEGER(IntKi) :: idx_qp !< index of current quadrature point - m%qp%Mi(:,:,idx_qp,nelem) = 0.0_BDKi + do idx_qp=1,p%nqp - ! Set diagonal values for mass - DO i=1,3 - m%qp%Mi(i,i,idx_qp,nelem) = p%qp%mmm(idx_qp,nelem) - ENDDO + m%qp%Mi(:,:,idx_qp,nelem) = 0.0_BDKi - ! set mass-inertia coupling terms - m%qp%Mi(1:3,4:6,idx_qp,nelem) = -SkewSymMat(m%qp%RR0mEta(:,idx_qp,nelem)) - m%qp%Mi(4:6,1:3,idx_qp,nelem) = SkewSymMat(m%qp%RR0mEta(:,idx_qp,nelem)) + ! Set diagonal values for mass + DO i=1,3 + m%qp%Mi(i,i,idx_qp,nelem) = p%qp%mmm(idx_qp,nelem) + ENDDO - ! Set inertia terms - m%qp%Mi(4:6,4:6,idx_qp,nelem) = m%qp%rho(:,:,idx_qp,nelem) + ! set mass-inertia coupling terms + m%qp%Mi(1:3,4:6,idx_qp,nelem) = -SkewSymMat(m%qp%RR0mEta(:,idx_qp,nelem)) + m%qp%Mi(4:6,1:3,idx_qp,nelem) = SkewSymMat(m%qp%RR0mEta(:,idx_qp,nelem)) + ! Set inertia terms + m%qp%Mi(4:6,4:6,idx_qp,nelem) = m%qp%rho(:,:,idx_qp,nelem) + + end do + END SUBROUTINE BD_InertialMassMatrix !----------------------------------------------------------------------------------------------------------------------------------- !> This subroutine computes gyroscopic forces -!! Returns new value for `m%qp%Fb` +!! Returns new value for `m%qp%Fb` Note that the equations here are the inertial equations with the acceleration terms set to zero. ! called by BD_ElementMatrixAcc and BD_ElementMatrixForce -SUBROUTINE BD_GyroForce( nelem, idx_qp, m ) +SUBROUTINE BD_GyroForce( nelem, p, m ) INTEGER(IntKi), INTENT(IN ) :: nelem !< index of current element in loop - INTEGER(IntKi), INTENT(IN ) :: idx_qp !< index of current quadrature point - TYPE(BD_MiscVarType), INTENT(INOUT) :: m !< misc/optimization variables + TYPE(BD_ParameterType), INTENT(IN ) :: p !< Parameters + TYPE(BD_MiscVarType), INTENT(INOUT) :: m !< misc/optimization variables - REAL(BDKi) :: Bi(6,6) - REAL(BDKi) :: ome(3) - REAL(BDKi) :: temp33(3,3) +! REAL(BDKi) :: Bi(6,6) + REAL(BDKi) :: beta(3) + REAL(BDKi) :: gama(3) + INTEGER(IntKi) :: idx_qp !< index of current quadrature point - ! Initialize Bi to zero (left quadrants stay zero) - Bi= 0.0_BDKi + m%qp%Fb(:,:,nelem) = 0.0_BDKi - ! angular velocity terms - ome(:) = m%qp%vvv(4:6,idx_qp,nelem) - ! coupling of angular velocity to mEta (?) - temp33 = -MATMUL(SkewSymMat(ome),SkewSymMat(m%qp%RR0mEta(:,idx_qp,nelem))) - Bi(1:3,4:6) = temp33 + do idx_qp=1,p%nqp + + beta = cross_product(m%qp%vvv(4:6,idx_qp,nelem), m%qp%RR0mEta(:,idx_qp,nelem)) !MATMUL(SkewSymMat(ome),mEta) + gama = MATMUL(m%qp%rho(:,:,idx_qp,nelem),m%qp%vvv(4:6,idx_qp,nelem)) - ! angular velocity times inertia tensor - temp33 = MATMUL(SkewSymMat(ome),m%qp%rho(:,:,idx_qp,nelem)) - Bi(4:6,4:6) = temp33 + !Compute Fb (Equation 22 in Wang_2014 with aaa = 0) + m%qp%Fb(1:3,idx_qp,nelem) = cross_product(m%qp%vvv(4:6,idx_qp,nelem), beta) + m%qp%Fb(4:6,idx_qp,nelem) = cross_product(m%qp%vvv(4:6,idx_qp,nelem), gama) - ! resulting gyroscopic force - m%qp%Fb(:,idx_qp,nelem) = MATMUL(Bi,m%qp%vvv(1:6,idx_qp,nelem)) + end do END SUBROUTINE BD_GyroForce @@ -3242,6 +3370,7 @@ SUBROUTINE BD_StaticSolution( x, gravity, u, p, m, piter, ErrStat, ErrMsg ) ENDDO + CALL setErrStat( ErrID_Fatal, "Solution does not converge after the maximum number of iterations", ErrStat, ErrMsg, RoutineName) RETURN @@ -3307,7 +3436,7 @@ SUBROUTINE BD_GenerateStaticElement( x, gravity, u, p, m ) CALL BD_StaticElementMatrix( nelem, gravity, p, x, m ) CALL BD_AssembleStiffK(nelem,p,m%elk,m%StifK) - CALL BD_AssembleRHS(nelem,p,m%elf(:,:,nelem),m%RHS) + CALL BD_AssembleRHS(nelem,p,m%elf,m%RHS) ENDDO @@ -3324,105 +3453,207 @@ SUBROUTINE BD_StaticElementMatrix( nelem, gravity, p, x, m ) TYPE(BD_ContinuousStateType), INTENT(IN ) :: x !< Continuous states at t TYPE(BD_MiscVarType), INTENT(INOUT) :: m !< misc/optimization variables - INTEGER(IntKi) :: idx_qp INTEGER(IntKi) :: i INTEGER(IntKi) :: j - INTEGER(IntKi) :: idx_dof1 - INTEGER(IntKi) :: idx_dof2 + INTEGER(IntKi) :: idx_dof1, idx_dof2 + INTEGER(IntKi) :: idx_qp CHARACTER(*), PARAMETER :: RoutineName = 'BD_StaticElementMatrix' - m%elk(:,:,:,:) = 0.0_BDKi - m%elf(:,:,:) = 0.0_BDKi - - - - DO idx_qp=1,p%nqp - - CALL BD_ElasticForce( nelem,idx_qp,p,m,.true. ) ! Calculate Fc, Fd [and Oe, Pe, and Qe for N-R algorithm] - CALL BD_GravityForce( nelem, idx_qp, p, m, gravity ) - m%qp%Ftemp(:,idx_qp,nelem) = m%qp%Fd(:,idx_qp,nelem) - m%qp%Fg(:,idx_qp,nelem) - m%DistrLoad_QP(:,idx_qp,nelem) + CALL BD_ElasticForce( nelem,p,m,.true. ) ! Calculate Fc, Fd [and Oe, Pe, and Qe for N-R algorithm] + CALL BD_GravityForce( nelem,p,m,gravity ) ! Calculate Fg - DO j=1,p%nodes_per_elem - DO idx_dof2=1,p%dof_node - DO i=1,p%nodes_per_elem - DO idx_dof1=1,p%dof_node - m%elk(idx_dof1,i,idx_dof2,j) = m%elk(idx_dof1,i,idx_dof2,j) + p%Shp(i,idx_qp) *m%qp%Qe(idx_dof1,idx_dof2,idx_qp,nelem) *p%Shp(j,idx_qp) *p%Jacobian(idx_qp,nelem)*p%QPtWeight(idx_qp) - m%elk(idx_dof1,i,idx_dof2,j) = m%elk(idx_dof1,i,idx_dof2,j) + p%Shp(i,idx_qp) *m%qp%Pe(idx_dof1,idx_dof2,idx_qp,nelem) *p%ShpDer(j,idx_qp)*p%QPtWeight(idx_qp) - m%elk(idx_dof1,i,idx_dof2,j) = m%elk(idx_dof1,i,idx_dof2,j) + p%ShpDer(i,idx_qp)*m%qp%Oe(idx_dof1,idx_dof2,idx_qp,nelem) *p%Shp(j,idx_qp) *p%QPtWeight(idx_qp) - m%elk(idx_dof1,i,idx_dof2,j) = m%elk(idx_dof1,i,idx_dof2,j) + p%ShpDer(i,idx_qp)*m%qp%Stif(idx_dof1,idx_dof2,idx_qp,nelem)*p%ShpDer(j,idx_qp)*p%QPtWeight(idx_qp)/p%Jacobian(idx_qp,nelem) - ENDDO + m%qp%Ftemp(:,:,nelem) = m%qp%Fd(:,:,nelem) - m%qp%Fg(:,:,nelem) - m%DistrLoad_QP(:,:,nelem) + + DO j=1,p%nodes_per_elem + DO idx_dof2=1,p%dof_node + DO i=1,p%nodes_per_elem + DO idx_dof1=1,p%dof_node + m%elk(idx_dof1,i,idx_dof2,j) = 0.0_BDKi + DO idx_qp = 1,p%nqp ! dot_product( m%qp%Qe( idx_dof1,idx_dof2,:,nelem), p%QPtw_Shp_Shp_Jac( :,i,j,nelem)) + m%elk(idx_dof1,i,idx_dof2,j) = m%elk(idx_dof1,i,idx_dof2,j) + m%qp%Qe( idx_dof1,idx_dof2,idx_qp,nelem)*p%QPtw_Shp_Shp_Jac(idx_qp,i,j,nelem) + END DO + + DO idx_qp = 1,p%nqp ! dot_product( m%qp%Pe( idx_dof1,idx_dof2,:,nelem), p%QPtw_Shp_ShpDer( :,i,j) ) + m%elk(idx_dof1,i,idx_dof2,j) = m%elk(idx_dof1,i,idx_dof2,j) + m%qp%Pe( idx_dof1,idx_dof2,idx_qp,nelem)*p%QPtw_Shp_ShpDer(idx_qp,i,j) + END DO + DO idx_qp = 1,p%nqp ! dot_product( m%qp%Oe( idx_dof1,idx_dof2,:,nelem), p%QPtw_Shp_ShpDer( :,j,i) ) + m%elk(idx_dof1,i,idx_dof2,j) = m%elk(idx_dof1,i,idx_dof2,j) + m%qp%Oe( idx_dof1,idx_dof2,idx_qp,nelem)*p%QPtw_Shp_ShpDer(idx_qp,j,i) + END DO + DO idx_qp = 1,p%nqp ! dot_product( m%qp%Stif(idx_dof1,idx_dof2,:,nelem), p%QPtw_ShpDer_ShpDer_Jac(:,i,j,nelem)) + m%elk(idx_dof1,i,idx_dof2,j) = m%elk(idx_dof1,i,idx_dof2,j) + m%qp%Stif(idx_dof1,idx_dof2,idx_qp,nelem)*p%QPtw_ShpDer_ShpDer_Jac(idx_qp,i,j,nelem) + END DO ENDDO ENDDO ENDDO + ENDDO + DO i=1,p%nodes_per_elem + DO idx_dof1=1,p%dof_node + m%elf(idx_dof1,i) = 0.0_BDKi + DO idx_qp = 1,p%nqp ! dot_product( m%qp%Fc (idx_dof1,:,nelem), p%QPtw_ShpDer( :,i)) + m%elf(idx_dof1,i) = m%elf(idx_dof1,i) - m%qp%Fc (idx_dof1,idx_qp,nelem)*p%QPtw_ShpDer(idx_qp,i) + END DO + DO idx_qp = 1,p%nqp ! dot_product(m%qp%Ftemp(idx_dof1,:,nelem), p%QPtw_Shp_Jac(:,i,nelem) ) + m%elf(idx_dof1,i) = m%elf(idx_dof1,i) - m%qp%Ftemp(idx_dof1,idx_qp,nelem)*p%QPtw_Shp_Jac(idx_qp,i,nelem) + END DO - DO i=1,p%nodes_per_elem - DO j=1,p%dof_node - m%elf(j,i,nelem) = m%elf(j,i,nelem) - p%ShpDer(i,idx_qp)*m%qp%Fc (j,idx_qp,nelem)*p%QPtWeight(idx_qp) - m%elf(j,i,nelem) = m%elf(j,i,nelem) - p%Shp(i,idx_qp) *m%qp%Ftemp(j,idx_qp,nelem)*p%Jacobian(idx_qp,nelem)*p%QPtWeight(idx_qp) - ENDDO - ENDDO + ! Elastic force internal + m%EFint(idx_dof1,i,nelem) = 0.0_BDKi + DO idx_qp = 1,p%nqp ! dot_product(m%qp%Fc(idx_dof1,:,nelem), p%QPtw_ShpDer(:,i)) + m%EFint(idx_dof1,i,nelem) = m%EFint(idx_dof1,i,nelem) - m%qp%Fc(idx_dof1,idx_qp,nelem)*p%QPtw_ShpDer(idx_qp,i) + END DO +! m%EFint(idx_dof1,i,nelem) = m%EFint(idx_dof1,i,nelem) - dot_product(m%qp%Fd(idx_dof1,:,nelem), p%QPtw_Shp_Jac(:,i,nelem)) + ENDDO ENDDO - RETURN END SUBROUTINE BD_StaticElementMatrix + !----------------------------------------------------------------------------------------------------------------------------------- ! This subroutine calculates the internal nodal forces at each finite-element -! nodes along beam axis \n +! nodes along beam axis for the static case. This is more involved than in the dynamic case because m%EFint is not calculated beforehand. ! Nodal forces = K u -SUBROUTINE BD_GenerateInternalForceMoment( y, p, m ) +!FIXME: NOTE: if we go to multiple elements for trap quadrature, we will need to double check this routine. +SUBROUTINE BD_InternalForceMoment( x, p, m ) - TYPE(BD_OutputType), INTENT(IN ) :: y !< Outputs computed at t (needed for position information) + TYPE(BD_ContinuousStateType), INTENT(IN ) :: x !< Continuous states at t TYPE(BD_ParameterType), INTENT(IN ) :: p !< Parameters TYPE(BD_MiscVarType), INTENT(INOUT) :: m !< misc/optimization variables INTEGER(IntKi) :: nelem ! number of current element INTEGER(IntKi) :: idx_node_in_elem INTEGER(IntKi) :: idx_node - CHARACTER(*), PARAMETER:: RoutineName = 'BD_GenerateInternalForceMoment' +! INTEGER(IntKi) :: idx_node_next + INTEGER(IntKi) :: idx_qp + REAL(BDKi) :: Tmp6(6) + REAL(BDKi) :: Tmp3(3) + REAL(BDKi) :: PrevNodePos(3) + INTEGER(IntKi) :: i !< generic counter + INTEGER(IntKi) :: LastNode !< Last node in element to consider in integration in FE points + INTEGER(IntKi) :: StartNode !< First node to consider in integration for QP points + CHARACTER(*), PARAMETER:: RoutineName = 'BD_InternalForceMoment' + + + INTEGER(IntKi) :: FirstQPinFEpoint(p%nodes_per_elem,p%elem_total) + REAL(BDKi) :: FEweight(p%nodes_per_elem,p%elem_total) + REAL(BDKi) :: SumShp + + + FirstQPinFEpoint = p%nqp+1 ! Set to just past the last QP in the event that the last QP is just beyond the last FE location due to numerical noise. + FEweight= 0.0_BDKi +!FIXME: move this to a parameter when it is working as expected + ! First we find which QP points are the first QP to consider for each node in each element + DO nelem=1,p%elem_total + DO i=1,p%nodes_per_elem + SumShp=0.0_BDKi + DO idx_qp=p%nqp,1,-1 ! Step inwards to find the first QP past the FE point + IF ( TwoNorm(p%uu0(1:3,idx_qp,nelem)) >= TwoNorm(p%uuN0(1:3,i,nelem))) THEN + FirstQPinFEpoint(i,nelem)=idx_qp + FEweight(i,nelem) = FEweight(i,nelem) + p%Shp(i,idx_qp) !*p%Shp(j,idx_qp) + ENDIF + SumShp=SumShp+p%Shp(i,idx_qp) !*p%Shp(j,idx_qp) + ENDDO + FEweight(i,nelem) = FEweight(i,nelem) / SumShp + ENDDO + ! Tip contribution + ! Setting FEWeight at the tip to 1. The contirbution of the tip of each element should be absolute, hence no weighting is required + FEweight(p%nodes_per_elem,nelem) = 1.0_BDKi + ENDDO ! Initialize all values to zero. - m%BldInternalForce(:,:) = 0.0_ReKi + m%BldInternalForceFE(:,:) = 0.0_BDKi + m%EFint(:,:,:) = 0.0_BDKi !FIXME: these need to be computed based on jmj's equations - ! Set values for tip - idx_node=y%BldMotion%NNodes + ! Integrate the elastic force contributions from the tip inwards. We only consider the shape function contributions at each QP beyond the current FE node. + ! Note that FE node contributions at the start of an element should be contained in the last node of the preceeding element. + DO nelem=p%elem_total,1,-1 + DO i=p%nodes_per_elem,1,-1 + ! Integrate shape functions across the quadrature points (only consider the portion of the shape function beyond the current FE point location) + DO idx_qp=p%nqp,1,-1 + ! Force contributions from current node + m%EFint(1:3,i,nelem) = m%EFint(1:3,i,nelem) & + + p%ShpDer(i,idx_qp)*m%qp%Fc(1:3,idx_qp,nelem)*p%QPtWeight(idx_qp) + + ! Moment contributions from current node + m%EFint(4:6,i,nelem) = m%EFint(4:6,i,nelem) & + + p%ShpDer(i,idx_qp)*m%qp%Fc(4:6,idx_qp,nelem)*p%QPtWeight(idx_qp) & + + p%Shp(i,idx_qp) *m%qp%Fd(4:6,idx_qp,nelem)*p%QPtWeight(idx_qp)*p%Jacobian(idx_qp,nelem) ! Fd only contains moments + ENDDO + ENDDO + ENDDO - ! Force at tip - m%BldInternalForce(1:3,idx_node) = -MATMUL(p%GlbRot,m%elf(1:3,p%nodes_per_elem,p%elem_total)) + ! Calculate the internal forces and moments at the FE nodes. + ! NOTE: the elastic force contributions are already calculated and stored, so we merely need to perform the summation. + ! NOTE: we are only counting unique points, not overlapping FE points (those are identical as the first node is not a state) + ! The p%node_elem_idx stores the first and last nodes of elements in the aggregated nodes (ignoring overlapping nodes) - ! Moment at tip - m%BldInternalForce(4:6,idx_node) = -MATMUL(p%GlbRot,m%elf(4:6,p%nodes_per_elem,p%elem_total)) + ! Working from tip to root + LastNode = p%nodes_per_elem-1 ! Already counted tip, so set the last node for iteration loop - - ! To do the integration, start at the tip and work inwards towards the root. (BJJ: why? there is no summation below....) - - ! sectional loads were already computed at the FE nodes - - DO idx_node=size(p%NdIndx)-1,1,-1 + DO nelem = p%elem_total,1,-1 - ! Set the index to the m%elf in current element - nelem = p%OutNd2NdElem(2,idx_node) - idx_node_in_elem = p%OutNd2NdElem(1,idx_node) - - ! rotate to the correct coordinate system - m%BldInternalForce(1:3,idx_node) = -MATMUL(p%GlbRot,m%elf(1:3,idx_node_in_elem,nelem)) - m%BldInternalForce(4:6,idx_node) = -MATMUL(p%GlbRot,m%elf(4:6,idx_node_in_elem,nelem)) - + m%BldInternalForceFE(1:3,nelem*LastNode+1) = FEweight(p%nodes_per_elem,nelem) * m%EFint(1:3,p%nodes_per_elem,nelem) + m%BldInternalForceFE(4:6,nelem*LastNode+1) = m%EFint(4:6,p%nodes_per_elem,nelem) + + ! Keep track of previous node for adding force contributions to moments + PrevNodePos = p%uuN0(1:3,p%nodes_per_elem,nelem) + x%q(1:3,nelem*LastNode+1) + + DO idx_node_in_elem=LastNode,2,-1 ! Skip first node on element as it corresponds to last node of previous element or the root + + ! Index to node + idx_node = p%node_elem_idx(nelem,1)-1 + idx_node_in_elem ! p%node_elem_idx(nelem,1) is the first node in the element + + ! Force term + m%BldInternalForceFE(1:3,idx_node) = FEweight(idx_node_in_elem,nelem) * m%EFint(1:3,idx_node_in_elem,nelem) + m%BldInternalForceFE(1:3,idx_node+1) + (1.0_BDKi - FEweight(idx_node_in_elem+1,nelem)) * m%EFint(1:3,idx_node_in_elem+1,nelem) + + ! Moment term including forces from next node out projected to this node + ! NOTE: this appears like double counting, but the issue is that the Fd and Fc terms are both included in EFint. + ! These terms partially cancel each other. Fd includes part of the force term as well. The exact math + ! is a bit fuzzy to me. It appears to work though. If only Fc is used, then the result is off a very + ! small amount that is unknown to me. + Tmp3 = PrevNodePos - (p%uuN0(1:3,idx_node_in_elem,nelem) + x%q(1:3,idx_node)) + m%BldInternalForceFE(4:6,idx_node) = m%EFint(4:6,idx_node_in_elem,nelem) + m%BldInternalForceFE(4:6,idx_node+1) + cross_product( Tmp3, m%BldInternalForceFE(1:3,idx_node+1) ) + + ! Keep track of node position next node in. + PrevNodePos = p%uuN0(1:3,idx_node_in_elem,nelem) + x%q(1:3,idx_node) + + ENDDO + ! only skip last node of outer element. So next element towards root needs its last node + ! included as first node of next element is skipped (they overlap) + !LastNode=p%nodes_per_elem ENDDO - + ! Now deal with the root node + ! For the root node: the root reaction force is contained in the m%RHS (root node is not a state) + IF(p%analysis_type .EQ. BD_DYNAMIC_ANALYSIS) THEN + ! Add root reaction + m%BldInternalForceFE(1:3,1) = -m%RHS(1:3,1) + m%BldInternalForceFE(4:6,1) = -m%RHS(4:6,1) + ELSE + ! Add root reaction -- This is in the first node for static case and does not need contributions from the outboard sections due + ! to how the solve is performed, but it is negative. + m%BldInternalForceFE(1:3,1) = -m%EFint(1:3,1,1) + m%BldInternalForceFE(4:6,1) = -m%EFint(4:6,1,1) + ENDIF + + ! Rotate coords + DO i=1,SIZE(m%BldInternalForceFE,DIM=2) + m%BldInternalForceFE(1:3,i) = MATMUL(p%GlbRot,m%BldInternalForceFE(1:3,i)) + m%BldInternalForceFE(4:6,i) = MATMUL(p%GlbRot,m%BldInternalForceFE(4:6,i)) + ENDDO + RETURN -END SUBROUTINE BD_GenerateInternalForceMoment +END SUBROUTINE BD_InternalForceMoment + @@ -3465,14 +3696,14 @@ SUBROUTINE BD_GA2(t,n,u,utimes,p,x,xd,z,OtherState,m,ErrStat,ErrMsg) end if - ! initialize accelerations here: + ! on first call, initialize accelerations at t=0 and put mass and stiffness matrices in the summary file: if ( n .EQ. 0 .or. .not. OtherState%InitAcc) then - + + ! initialize accelerations at t call BD_Input_extrapinterp( u, utimes, u_interp, t, ErrStat2, ErrMsg2 ) - call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName) - + call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName) - ! Transform quantities from global frame to local (blade) frame + ! Transform quantities from global frame to local (blade) frame CALL BD_InputGlobalLocal(p,u_interp) ! Copy over the DistrLoads @@ -3481,18 +3712,21 @@ SUBROUTINE BD_GA2(t,n,u,utimes,p,x,xd,z,OtherState,m,ErrStat,ErrMsg) ! Incorporate boundary conditions CALL BD_BoundaryGA2(x,p,u_interp,OtherState) - - CALL BD_InitAcc( u_interp, p, x, OtherState, m, ErrStat2, ErrMsg2) - call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName) - + ! initialize the accelerations + CALL BD_InitAcc( u_interp, p, x, m, OtherState%Acc, ErrStat2, ErrMsg2) + call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName) if (ErrStat >= AbortErrLev) then call cleanup() return end if - + + ! initialize GA2 algorithm acceleration variable (acts as a filtering value on OtherState%acc) + OtherState%Xcc(:,:) = OtherState%Acc(:,:) + ! accelerations have been initialized OtherState%InitAcc = .true. - + + ! If we are writing to a summary file if (m%Un_Sum > 0) then @@ -3507,8 +3741,8 @@ SUBROUTINE BD_GA2(t,n,u,utimes,p,x,xd,z,OtherState,m,ErrStat,ErrMsg) WRITE(m%Un_Sum,'()') CALL WrMatrix(RESHAPE(m%MassM, (/p%dof_total, p%dof_total/)), m%Un_Sum, p%OutFmt, 'Full mass matrix (IEC coordinates)') - CLOSE(m%Un_Sum) - m%Un_Sum = -1 + CLOSE(m%Un_Sum) + m%Un_Sum = -1 end if end if @@ -3525,16 +3759,16 @@ SUBROUTINE BD_GA2(t,n,u,utimes,p,x,xd,z,OtherState,m,ErrStat,ErrMsg) CALL PitchActuator_SetBC(p, u_interp, xd) ENDIF - ! GA2: prediction - CALL BD_TiSchmPredictorStep( x, OtherState, p ) ! updates x and OtherState (from values at t to predictions at t+dt?) - - ! Transform quantities from global frame to local (blade in BD coords) frame CALL BD_InputGlobalLocal(p,u_interp) ! Copy over the DistrLoads CALL BD_DistrLoadCopy( p, u_interp, m ) + + ! GA2: prediction + CALL BD_TiSchmPredictorStep( x, OtherState, p ) ! updates x and OtherState accelerations (from values at t to predictions at t+dt) + ! Incorporate boundary conditions (overwrite first node of continuous states and Acc array at t+dt) CALL BD_BoundaryGA2(x,p,u_interp,OtherState) @@ -3570,14 +3804,15 @@ SUBROUTINE BD_TiSchmPredictorStep( x, OtherState, p ) !X DO i=1,p%node_total DO i=2,p%node_total - tr = p%dt * x%dqdt(:,i) + p%coef(1) * OtherState%acc(:,i) + p%coef(2) * OtherState%xcc(:,i) - x%dqdt(:,i) = x%dqdt(:,i) + p%coef(3) * OtherState%acc(:,i) + p%coef(4) * OtherState%xcc(:,i) ! velocities at t+dt - OtherState%xcc(:,i) = p%coef(5) * OtherState%acc(:,i) + p%coef(6) * OtherState%xcc(:,i) ! xcc accelerations at t+dt - OtherState%acc(:,i) = 0.0_BDKi + tr = p%dt * x%dqdt(:,i) + p%coef(1) * OtherState%acc(:,i) + p%coef(2) * OtherState%xcc(:,i) ! displacements at t+dt + x%dqdt(:,i) = x%dqdt(:,i) + p%coef(3) * OtherState%acc(:,i) + p%coef(4) * OtherState%xcc(:,i) ! velocities at t+dt + OtherState%xcc(:,i) = p%coef(5) * OtherState%acc(:,i) + p%coef(6) * OtherState%xcc(:,i) ! xcc accelerations at t+dt: ((1-alpha_m)*xcc_(t+dt) = (1-alpha_f)*Acc_(t+dt) + alpha_f*Acc_t - alpha_m*xcc_t + OtherState%acc(:,i) = 0.0_BDKi ! acc accelerations at t+dt - x%q(1:3,i) = x%q(1:3,i) + tr(1:3) ! position at t+dt + x%q(1:3,i) = x%q(1:3,i) + tr(1:3) ! position at t+dt - CALL BD_CrvCompose(rot_temp, tr(4:6), x%q(4:6,i), FLAG_R1R2) ! rot_temp = tr(4:6) composed with x%q(4:6,i) [rotations at t], is the output + ! tr does not contain w-m parameters, yet we are treating them as such + CALL BD_CrvCompose(rot_temp, tr(4:6), x%q(4:6,i), FLAG_R1R2) ! rot_temp = tr(4:6) composed with x%q(4:6,i) [rotations at t], is the output x%q(4:6,i) = rot_temp ! rotations at t+dt ENDDO @@ -3596,21 +3831,21 @@ SUBROUTINE BD_TiSchmComputeCoefficients(p) REAL(DbKi) :: tr0 REAL(DbKi) :: tr1 REAL(DbKi) :: tr2 - REAL(DbKi) :: alfam - REAL(DbKi) :: alfaf + REAL(DbKi) :: alfam ! \alpha_M + REAL(DbKi) :: alfaf ! \alpha_F REAL(DbKi) :: gama REAL(DbKi) :: beta - REAL(DbKi) :: oalfaM - REAL(DbKi) :: deltat2 ! dt^2 + REAL(DbKi) :: oalfaM ! 1 - \alpha_M + REAL(DbKi) :: deltat2 ! {\delta t}^2 = dt^2 tr0 = p%rhoinf + 1.0_BDKi alfam = (2.0_BDKi * p%rhoinf - 1.0_BDKi) / tr0 alfaf = p%rhoinf / tr0 gama = 0.5_BDKi - alfam + alfaf - beta = 0.25 * (1.0_BDKi - alfam + alfaf) * (1.0_BDKi - alfam + alfaf) + beta = 0.25 * (1.0_BDKi - alfam + alfaf)**2 - deltat2 = p%dt * p%dt + deltat2 = p%dt**2 oalfaM = 1.0_BDKi - alfam tr0 = alfaf / oalfaM tr1 = alfam / oalfaM @@ -3664,9 +3899,9 @@ END SUBROUTINE BD_BoundaryGA2 !FIXME: note similarities to BD_StaticSolution. May be able to combine SUBROUTINE BD_DynamicSolutionGA2( x, OtherState, u, p, m, ErrStat, ErrMsg) - TYPE(BD_ContinuousStateType), INTENT(INOUT) :: x !< Continuous states at t on input at t + dt on output - TYPE(BD_OtherStateType), INTENT(INOUT) :: OtherState !< Other states at t on input; at t+dt on outputs - TYPE(BD_InputType), INTENT(IN ) :: u + TYPE(BD_ContinuousStateType), INTENT(INOUT) :: x !< Continuous states: input are the predicted values at t+dt; output are calculated values at t + dt + TYPE(BD_OtherStateType), INTENT(INOUT) :: OtherState !< Other states: input are the predicted accelerations at t+dt; output are calculated values at t + dt + TYPE(BD_InputType), INTENT(IN ) :: u !< inputs in the local coordinate system (not FAST's global system) TYPE(BD_ParameterType), INTENT(IN ) :: p !< Parameters TYPE(BD_MiscVarType), INTENT(INOUT) :: m !< misc/optimization variables INTEGER(IntKi), INTENT( OUT) :: ErrStat !< Error status of the operation @@ -3675,6 +3910,7 @@ SUBROUTINE BD_DynamicSolutionGA2( x, OtherState, u, p, m, ErrStat, ErrMsg) INTEGER(IntKi) :: ErrStat2 ! Temporary Error status CHARACTER(ErrMsgLen) :: ErrMsg2 ! Temporary Error message CHARACTER(*), PARAMETER :: RoutineName = 'BD_DynamicSolutionGA2' + REAL(DbKi) :: Eref REAL(DbKi) :: Enorm INTEGER(IntKi) :: i @@ -3688,15 +3924,11 @@ SUBROUTINE BD_DynamicSolutionGA2( x, OtherState, u, p, m, ErrStat, ErrMsg) Eref = 0.0_BDKi DO i=1,p%niter - IF(MOD(i-1,p%n_fact) .EQ. 0) THEN - fact = .TRUE. - ELSE - fact = .FALSE. - ENDIF + fact = MOD(i-1,p%n_fact) .EQ. 0 ! when true, we factor the jacobian matrix ! Apply accelerations using F=ma ? Is that what this routine does? ! Calculate Quadrature point values needed - CALL BD_QuadraturePointData( p,x,m ) ! Calculate QP values uuu, uup, RR0, kappa, E1 + CALL BD_QuadraturePointData( p,x,m ) ! Calculate QP values uuu, uup, RR0, kappa, E1 using current guess at continuous states (displacements and velocities) CALL BD_GenerateDynamicElementGA2( x, OtherState, u, p, m,fact) ! Apply additional forces / loads at GLL points (such as aerodynamic loads)? @@ -3717,7 +3949,8 @@ SUBROUTINE BD_DynamicSolutionGA2( x, OtherState, u, p, m, ErrStat, ErrMsg) ! note m%LP_indx is allocated larger than necessary (to allow us to use it in multiple places) - CALL LAPACK_getrf( p%dof_total-6, p%dof_total-6, m%LP_StifK_LU, m%LP_indx, ErrStat2, ErrMsg2); CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + CALL LAPACK_getrf( p%dof_total-6, p%dof_total-6, m%LP_StifK_LU, m%LP_indx, ErrStat2, ErrMsg2) + CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) if (ErrStat >= AbortErrLev) return ENDIF @@ -3727,11 +3960,13 @@ SUBROUTINE BD_DynamicSolutionGA2( x, OtherState, u, p, m, ErrStat, ErrMsg) m%LP_RHS_LU = m%LP_RHS(7:p%dof_total) ! Solve for X in A*X=B to get the displacement of blade - CALL LAPACK_getrs( 'N',p%dof_total-6, m%LP_StifK_LU, m%LP_indx, m%LP_RHS_LU, ErrStat2, ErrMsg2); CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + CALL LAPACK_getrs( 'N',p%dof_total-6, m%LP_StifK_LU, m%LP_indx, m%LP_RHS_LU, ErrStat2, ErrMsg2) + CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) - m%Solution(:,1) = 0.0_BDKi ! first node is not set below + m%Solution(:,1) = 0.0_BDKi ! first node is not set below. By definition, there is no displacement of the first node. m%Solution(:,2:p%node_total) = RESHAPE( m%LP_RHS_LU, (/ p%dof_node, (p%node_total - 1) /) ) + CALL BD_UpdateDynamicGA2(p,m,x,OtherState) ! Check for convergence Enorm = SQRT(abs(DOT_PRODUCT(m%LP_RHS_LU, m%LP_RHS(7:p%dof_total)))) @@ -3743,8 +3978,6 @@ SUBROUTINE BD_DynamicSolutionGA2( x, OtherState, u, p, m, ErrStat, ErrMsg) IF(Enorm .LE. Eref) RETURN ENDIF - CALL BD_UpdateDynamicGA2(p,m,x,OtherState) - ENDDO CALL setErrStat( ErrID_Fatal, "Solution does not converge after the maximum number of iterations", ErrStat, ErrMsg, RoutineName) @@ -3767,23 +4000,24 @@ SUBROUTINE BD_UpdateDynamicGA2( p, m, x, OtherState ) TYPE(BD_ContinuousStateType), INTENT(INOUT) :: x !< Continuous states at t on input at t + dt on output TYPE(BD_OtherStateType), INTENT(INOUT) :: OtherState !< Other states at t on input; at t+dt on outputs - REAL(BDKi) :: rotf_temp(3) REAL(BDKi) :: roti_temp(3) REAL(BDKi) :: rot_temp(3) INTEGER(IntKi) :: i CHARACTER(*), PARAMETER :: RoutineName = 'BD_UpdateDynamicGA2' - + ! m%Solution contains (\delta q dot dot) + ! The first node has no displacements by definition. DO i=2, p%node_total - x%q(1:3,i) = x%q(1:3,i) + p%coef(8) * m%Solution(1:3,i) - rotf_temp(1:3) = x%q(4:6,i) - roti_temp(1:3) = p%coef(8) * m%Solution(4:6,i) - CALL BD_CrvCompose(rot_temp,roti_temp,rotf_temp,FLAG_R1R2) ! rot_temp = roti_temp composed with rotf_temp - x%q(4:6,i) = rot_temp(1:3) + x%q(1:3,i) = x%q(1:3,i) + p%coef(8) * m%Solution(1:3,i) + + roti_temp = p%coef(8) * m%Solution(4:6,i) ! m%Solution(4:6,i) seems to contain accelerations (i.e., delta \omega dot), so I don't think this can be a w-m parameter + CALL BD_CrvCompose(rot_temp,roti_temp,x%q(4:6,i),FLAG_R1R2) ! rot_temp = roti_temp composed with x%q(4:6,i) + x%q(4:6,i) = rot_temp(1:3) + x%dqdt(:,i) = x%dqdt(:,i) + p%coef(7) * m%Solution(:,i) - OtherState%acc(:,i) = OtherState%acc(:,i) + m%Solution(:,i) ! Acceleration at t + dt - OtherState%xcc(:,i) = OtherState%xcc(:,i) + p%coef(9) * m%Solution(:,i) ! Acceleration at t + OtherState%acc(:,i) = OtherState%acc(:,i) + m%Solution(:,i) ! update acceleration (dqdtdt: q dot dot) to next guess for values at t+dt + OtherState%xcc(:,i) = OtherState%xcc(:,i) + p%coef(9) * m%Solution(:,i) ! update algorithm acceleration to next guess for values at t+dt: (1-alpha_m)*xcc_(n+1) = (1-alpha_f)*dqdtdt_(n+1) + alpha_f*dqdtdt_n - alpha_m*xcc_n ENDDO END SUBROUTINE BD_UpdateDynamicGA2 @@ -3808,9 +4042,13 @@ SUBROUTINE BD_GenerateDynamicElementGA2( x, OtherState, u, p, m, fact ) ! must initialize these because BD_AssembleStiffK and BD_AssembleRHS are INOUT m%RHS = 0.0_BDKi - m%StifK = 0.0_BDKi - m%MassM = 0.0_BDKi - m%DampG = 0.0_BDKi + + IF(fact) THEN + m%StifK = 0.0_BDKi + m%MassM = 0.0_BDKi + m%DampG = 0.0_BDKi + END IF + ! These values have not been set yet for the QP. @@ -3818,6 +4056,7 @@ SUBROUTINE BD_GenerateDynamicElementGA2( x, OtherState, u, p, m, fact ) CALL BD_QPData_mEta_rho( p,m ) ! Calculate the \f$ m \eta \f$ and \f$ \rho \f$ terms CALL BD_QPDataVelocity( p, x, m ) ! x%dqdt --> m%qp%vvv, m%qp%vvp + CALL BD_QPDataAcceleration( p, OtherState, m ) ! Naaa --> aaa (OtherState%Acc --> m%qp%aaa) DO nelem=1,p%elem_total @@ -3829,7 +4068,7 @@ SUBROUTINE BD_GenerateDynamicElementGA2( x, OtherState, u, p, m, fact ) CALL BD_AssembleStiffK(nelem,p,m%elm,m%MassM) CALL BD_AssembleStiffK(nelem,p,m%elg,m%DampG) ENDIF - CALL BD_AssembleRHS(nelem,p,m%elf(:,:,nelem),m%RHS) + CALL BD_AssembleRHS(nelem,p,m%elf,m%RHS) ENDDO RETURN @@ -3854,81 +4093,156 @@ SUBROUTINE BD_ElementMatrixGA2( fact, nelem, p, x, OtherState, m ) INTEGER(IntKi) :: idx_dof1 INTEGER(IntKi) :: idx_dof2 CHARACTER(*), PARAMETER :: RoutineName = 'BD_ElementMatrixGA2' - - m%elk(:,:,:,:) = 0.0_BDKi - m%elf(:,:,:) = 0.0_BDKi - m%elg(:,:,:,:) = 0.0_BDKi - m%elm(:,:,:,:) = 0.0_BDKi - - DO idx_qp=1,p%nqp + !FIXME: adp: I don't see the gyroscopic term in here. That is stored in m%qp%Fb !VA: The gyroscopic term is included in the m%qp%Gi. I think the m%qp%Fb term is used to calculate the accelerations - CALL BD_QPDataAcceleration( nelem, idx_qp, p, OtherState, m ) ! Naaa --> aaa + + + CALL BD_ElasticForce( nelem,p,m,fact ) ! Calculate Fc, Fd [and if(fact): Oe, Pe, and Qe for N-R algorithm] using m%qp%E1, m%qp%RR0, m%qp%kappa, m%qp%Stif + CALL BD_InertialForce( nelem,p,m,fact ) ! Calculate Fi [and Mi,Gi,Ki IF(fact)] + + IF(p%damp_flag .NE. 0) THEN + CALL BD_DissipativeForce( nelem,p,m,fact ) ! Calculate dissipative terms on Fc, Fd [and Sd, Od, Pd and Qd, betaC, Gd, Xd, Yd for N-R algorithm] + ENDIF + + CALL BD_GravityForce( nelem, p, m, p%gravity ) + + ! F^ext is combined with F^D (F^D = F^D-F^ext), i.e. RHS of Equation 9 in Wang_2014 + m%qp%Ftemp(:,:,nelem) = m%qp%Fd(:,:,nelem) - m%qp%Fg(:,:,nelem) - m%DistrLoad_QP(:,:,nelem) + - CALL BD_ElasticForce( nelem,idx_qp,p,m,fact ) ! Calculate Fc, Fd [and if(fact): Oe, Pe, and Qe for N-R algorithm] - CALL BD_InertialForce( nelem,idx_qp,p,m,fact ) ! Calculate Fi [and Mi,Gi,Ki IF(fact)] - IF(p%damp_flag .NE. 0) THEN - CALL BD_DissipativeForce( nelem, idx_qp, p,m,fact ) ! Calculate dissipative terms on Fc, Fd [and Sd, Od, Pd and Qd, betaC, Gd, Xd, Yd for N-R algorithm] - ENDIF - CALL BD_GravityForce( nelem, idx_qp, p, m, p%gravity ) + ! Equations 10, 11, 12 in Wang_2014 + + IF (fact) THEN + DO j=1,p%nodes_per_elem + DO idx_dof2=1,p%dof_node + DO i=1,p%nodes_per_elem + DO idx_dof1=1,p%dof_node + + m%elk(idx_dof1,i,idx_dof2,j) = 0.0_BDKi + DO idx_qp = 1,p%nqp ! dot_product(m%qp%Qe( idx_dof1,idx_dof2,:,nelem) + m%qp%Ki(idx_dof1,idx_dof2,:,nelem), p%QPtw_Shp_Shp_Jac( :,i,j,nelem) ) + m%elk(idx_dof1,i,idx_dof2,j) = m%elk(idx_dof1,i,idx_dof2,j) + (m%qp%Qe(idx_dof1,idx_dof2,idx_qp,nelem) + m%qp%Ki(idx_dof1,idx_dof2,idx_qp,nelem))*p%QPtw_Shp_Shp_Jac(idx_qp,i,j,nelem) + END DO + DO idx_qp = 1,p%nqp ! dot_product(m%qp%Pe( idx_dof1,idx_dof2,:,nelem) , p%QPtw_Shp_ShpDer( :,i,j) ) + m%elk(idx_dof1,i,idx_dof2,j) = m%elk(idx_dof1,i,idx_dof2,j) + m%qp%Pe( idx_dof1,idx_dof2,idx_qp,nelem)*p%QPtw_Shp_ShpDer(idx_qp,i,j) + END DO + DO idx_qp = 1,p%nqp ! dot_product(m%qp%Oe( idx_dof1,idx_dof2,:,nelem) , p%QPtw_Shp_ShpDer( :,j,i) ) + m%elk(idx_dof1,i,idx_dof2,j) = m%elk(idx_dof1,i,idx_dof2,j) + m%qp%Oe( idx_dof1,idx_dof2,idx_qp,nelem)*p%QPtw_Shp_ShpDer(idx_qp,j,i) + END DO + DO idx_qp = 1,p%nqp ! dot_product(m%qp%Stif(idx_dof1,idx_dof2,:,nelem) , p%QPtw_ShpDer_ShpDer_Jac(:,i,j,nelem) ) + m%elk(idx_dof1,i,idx_dof2,j) = m%elk(idx_dof1,i,idx_dof2,j) + m%qp%Stif(idx_dof1,idx_dof2,idx_qp,nelem)*p%QPtw_ShpDer_ShpDer_Jac(idx_qp,i,j,nelem) + END DO + + ENDDO + ENDDO + ENDDO + END DO - ! F^ext is combined with F^D (F^D = F^D-F^ext), i.e. RHS of Equation 9 in Wang_2014 - m%qp%Ftemp(:,idx_qp,nelem) = m%qp%Fd(:,idx_qp,nelem) - m%qp%Fg(:,idx_qp,nelem) - m%DistrLoad_QP(:,idx_qp,nelem) + DO j=1,p%nodes_per_elem + DO idx_dof2=1,p%dof_node + DO i=1,p%nodes_per_elem + DO idx_dof1=1,p%dof_node + + m%elm(idx_dof1,i,idx_dof2,j) = 0.0_BDKi + DO idx_qp = 1,p%nqp ! dot_product( m%qp%Mi(idx_dof1,idx_dof2,:,nelem), p%QPtw_Shp_Shp_Jac(:,i,j,nelem) ) + m%elm(idx_dof1,i,idx_dof2,j) = m%elm(idx_dof1,i,idx_dof2,j) + m%qp%Mi(idx_dof1,idx_dof2,idx_qp,nelem)*p%QPtw_Shp_Shp_Jac(idx_qp,i,j,nelem) + END DO + + ENDDO + ENDDO + ENDDO + END DO - ! Equations 10, 11, 12 in Wang_2014 - IF(fact) THEN + DO j=1,p%nodes_per_elem + DO idx_dof2=1,p%dof_node + DO i=1,p%nodes_per_elem + DO idx_dof1=1,p%dof_node + + m%elg(idx_dof1,i,idx_dof2,j) = 0.0_BDKi + DO idx_qp = 1,p%nqp ! dot_product( m%qp%Gi(idx_dof1,idx_dof2,:,nelem), p%QPtw_Shp_Shp_Jac(:,i,j,nelem)) + m%elg(idx_dof1,i,idx_dof2,j) = m%elg(idx_dof1,i,idx_dof2,j) + m%qp%Gi(idx_dof1,idx_dof2,idx_qp,nelem)*p%QPtw_Shp_Shp_Jac(idx_qp,i,j,nelem) + END DO + + ENDDO + ENDDO + ENDDO + END DO + + ! Dissipative terms + IF (p%damp_flag .NE. 0) THEN DO j=1,p%nodes_per_elem DO idx_dof2=1,p%dof_node DO i=1,p%nodes_per_elem DO idx_dof1=1,p%dof_node - m%elk(idx_dof1,i,idx_dof2,j) = m%elk(idx_dof1,i,idx_dof2,j) + p%Shp(i,idx_qp) *m%qp%Qe(idx_dof1,idx_dof2,idx_qp,nelem) *p%Shp(j,idx_qp) *p%Jacobian(idx_qp,nelem)*p%QPtWeight(idx_qp) - m%elk(idx_dof1,i,idx_dof2,j) = m%elk(idx_dof1,i,idx_dof2,j) + p%Shp(i,idx_qp) *m%qp%Pe(idx_dof1,idx_dof2,idx_qp,nelem) *p%ShpDer(j,idx_qp)*p%QPtWeight(idx_qp) - m%elk(idx_dof1,i,idx_dof2,j) = m%elk(idx_dof1,i,idx_dof2,j) + p%ShpDer(i,idx_qp)*m%qp%Oe(idx_dof1,idx_dof2,idx_qp,nelem) *p%Shp(j,idx_qp) *p%QPtWeight(idx_qp) - m%elk(idx_dof1,i,idx_dof2,j) = m%elk(idx_dof1,i,idx_dof2,j) + p%ShpDer(i,idx_qp)*m%qp%Stif(idx_dof1,idx_dof2,idx_qp,nelem)*p%ShpDer(j,idx_qp)*p%QPtWeight(idx_qp)/p%Jacobian(idx_qp,nelem) - m%elk(idx_dof1,i,idx_dof2,j) = m%elk(idx_dof1,i,idx_dof2,j) + p%Shp(i,idx_qp) *m%qp%Ki(idx_dof1,idx_dof2,idx_qp,nelem) *p%Shp(j,idx_qp) *p%Jacobian(idx_qp,nelem)*p%QPtWeight(idx_qp) - - m%elm(idx_dof1,i,idx_dof2,j) = m%elm(idx_dof1,i,idx_dof2,j) + p%Shp(i,idx_qp) *m%qp%Mi(idx_dof1,idx_dof2,idx_qp,nelem) *p%Shp(j,idx_qp) *p%Jacobian(idx_qp,nelem)*p%QPtWeight(idx_qp) - - m%elg(idx_dof1,i,idx_dof2,j) = m%elg(idx_dof1,i,idx_dof2,j) + p%Shp(i,idx_qp) *m%qp%Gi(idx_dof1,idx_dof2,idx_qp,nelem) *p%Shp(j,idx_qp) *p%Jacobian(idx_qp,nelem)*p%QPtWeight(idx_qp) + + DO idx_qp = 1,p%nqp ! dot_product(m%qp%Qd(idx_dof1,idx_dof2,:,nelem), p%QPtw_Shp_Shp_Jac( :,i,j,nelem)) + m%elk(idx_dof1,i,idx_dof2,j) = m%elk(idx_dof1,i,idx_dof2,j) + m%qp%Qd(idx_dof1,idx_dof2,idx_qp,nelem)*p%QPtw_Shp_Shp_Jac(idx_qp,i,j,nelem) + END DO + DO idx_qp = 1,p%nqp ! dot_product(m%qp%Pd(idx_dof1,idx_dof2,:,nelem), p%QPtw_Shp_ShpDer( :,i,j) ) + m%elk(idx_dof1,i,idx_dof2,j) = m%elk(idx_dof1,i,idx_dof2,j) + m%qp%Pd(idx_dof1,idx_dof2,idx_qp,nelem)*p%QPtw_Shp_ShpDer(idx_qp,i,j) + END DO + DO idx_qp = 1,p%nqp ! dot_product(m%qp%Od(idx_dof1,idx_dof2,:,nelem), p%QPtw_Shp_ShpDer( :,j,i) ) + m%elk(idx_dof1,i,idx_dof2,j) = m%elk(idx_dof1,i,idx_dof2,j) + m%qp%Od(idx_dof1,idx_dof2,idx_qp,nelem)*p%QPtw_Shp_ShpDer(idx_qp,j,i) + END DO + DO idx_qp = 1,p%nqp ! dot_product(m%qp%Sd(idx_dof1,idx_dof2,:,nelem), p%QPtw_ShpDer_ShpDer_Jac(:,i,j,nelem)) + m%elk(idx_dof1,i,idx_dof2,j) = m%elk(idx_dof1,i,idx_dof2,j) + m%qp%Sd(idx_dof1,idx_dof2,idx_qp,nelem)*p%QPtw_ShpDer_ShpDer_Jac(idx_qp,i,j,nelem) + END DO + ENDDO ENDDO ENDDO - ENDDO + END DO - ! Dissipative terms - IF(p%damp_flag .NE. 0) THEN - DO j=1,p%nodes_per_elem - DO idx_dof2=1,p%dof_node - DO i=1,p%nodes_per_elem - DO idx_dof1=1,p%dof_node - m%elk(idx_dof1,i,idx_dof2,j) = m%elk(idx_dof1,i,idx_dof2,j) + p%Shp(i,idx_qp) *m%qp%Qd(idx_dof1,idx_dof2,idx_qp,nelem) *p%Shp(j,idx_qp) *p%Jacobian(idx_qp,nelem)*p%QPtWeight(idx_qp) - m%elk(idx_dof1,i,idx_dof2,j) = m%elk(idx_dof1,i,idx_dof2,j) + p%Shp(i,idx_qp) *m%qp%Pd(idx_dof1,idx_dof2,idx_qp,nelem) *p%ShpDer(j,idx_qp)*p%QPtWeight(idx_qp) - m%elk(idx_dof1,i,idx_dof2,j) = m%elk(idx_dof1,i,idx_dof2,j) + p%ShpDer(i,idx_qp)*m%qp%Od(idx_dof1,idx_dof2,idx_qp,nelem) *p%Shp(j,idx_qp) *p%QPtWeight(idx_qp) - m%elk(idx_dof1,i,idx_dof2,j) = m%elk(idx_dof1,i,idx_dof2,j) + p%ShpDer(i,idx_qp)*m%qp%Sd(idx_dof1,idx_dof2,idx_qp,nelem) *p%ShpDer(j,idx_qp)*p%QPtWeight(idx_qp)/p%Jacobian(idx_qp,nelem) - - m%elg(idx_dof1,i,idx_dof2,j) = m%elg(idx_dof1,i,idx_dof2,j) + p%Shp(i,idx_qp) *m%qp%Xd(idx_dof1,idx_dof2,idx_qp,nelem) *p%Shp(j,idx_qp) *p%Jacobian(idx_qp,nelem)*p%QPtWeight(idx_qp) - m%elg(idx_dof1,i,idx_dof2,j) = m%elg(idx_dof1,i,idx_dof2,j) + p%Shp(i,idx_qp) *m%qp%Yd(idx_dof1,idx_dof2,idx_qp,nelem) *p%ShpDer(j,idx_qp)*p%QPtWeight(idx_qp) - m%elg(idx_dof1,i,idx_dof2,j) = m%elg(idx_dof1,i,idx_dof2,j) + p%ShpDer(i,idx_qp)*m%qp%Gd(idx_dof1,idx_dof2,idx_qp,nelem) *p%Shp(j,idx_qp) *p%QPtWeight(idx_qp) - m%elg(idx_dof1,i,idx_dof2,j) = m%elg(idx_dof1,i,idx_dof2,j) + p%ShpDer(i,idx_qp)*m%qp%betaC(idx_dof1,idx_dof2,idx_qp,nelem)*p%ShpDer(j,idx_qp)*p%QPtWeight(idx_qp)/p%Jacobian(idx_qp,nelem) - ENDDO + DO j=1,p%nodes_per_elem + DO idx_dof2=1,p%dof_node + DO i=1,p%nodes_per_elem + DO idx_dof1=1,p%dof_node + + DO idx_qp = 1,p%nqp ! dot_product(m%qp%Xd( idx_dof1,idx_dof2,:,nelem), p%QPtw_Shp_Shp_Jac( :,i,j,nelem)) + m%elg(idx_dof1,i,idx_dof2,j) = m%elg(idx_dof1,i,idx_dof2,j) + m%qp%Xd( idx_dof1,idx_dof2,idx_qp,nelem)*p%QPtw_Shp_Shp_Jac(idx_qp,i,j,nelem) + END DO + DO idx_qp = 1,p%nqp ! dot_product(m%qp%Yd( idx_dof1,idx_dof2,:,nelem), p%QPtw_Shp_ShpDer( :,i,j) ) + m%elg(idx_dof1,i,idx_dof2,j) = m%elg(idx_dof1,i,idx_dof2,j) + m%qp%Yd( idx_dof1,idx_dof2,idx_qp,nelem)*p%QPtw_Shp_ShpDer(idx_qp,i,j) + END DO + DO idx_qp = 1,p%nqp ! dot_product(m%qp%Gd( idx_dof1,idx_dof2,:,nelem), p%QPtw_Shp_ShpDer( :,j,i) ) + m%elg(idx_dof1,i,idx_dof2,j) = m%elg(idx_dof1,i,idx_dof2,j) + m%qp%Gd( idx_dof1,idx_dof2,idx_qp,nelem)*p%QPtw_Shp_ShpDer(idx_qp,j,i) + END DO + DO idx_qp = 1,p%nqp ! dot_product(m%qp%betaC(idx_dof1,idx_dof2,:,nelem), p%QPtw_ShpDer_ShpDer_Jac(:,i,j,nelem)) + m%elg(idx_dof1,i,idx_dof2,j) = m%elg(idx_dof1,i,idx_dof2,j) + m%qp%betaC(idx_dof1,idx_dof2,idx_qp,nelem)*p%QPtw_ShpDer_ShpDer_Jac(idx_qp,i,j,nelem) + END DO + ENDDO ENDDO ENDDO - ENDIF - - ENDIF + END DO + ENDIF ! add the dissipative terms + ENDIF + ! Equations 13 and 14 in Wang_2014. F^ext is combined with F^D (F^D = F^D-F^ext) - DO i=1,p%nodes_per_elem - DO idx_dof1=1,p%dof_node - m%elf(idx_dof1,i,nelem) = m%elf(idx_dof1,i,nelem) - p%ShpDer(i,idx_qp)*m%qp%Fc (idx_dof1,idx_qp,nelem)*p%QPtWeight(idx_qp) - m%elf(idx_dof1,i,nelem) = m%elf(idx_dof1,i,nelem) - p%Shp(i,idx_qp) *m%qp%Ftemp(idx_dof1,idx_qp,nelem)*p%Jacobian(idx_qp,nelem)*p%QPtWeight(idx_qp) - m%elf(idx_dof1,i,nelem) = m%elf(idx_dof1,i,nelem) - p%Shp(i,idx_qp) *m%qp%Fi (idx_dof1,idx_qp,nelem)*p%Jacobian(idx_qp,nelem)*p%QPtWeight(idx_qp) + + DO i=1,p%nodes_per_elem + DO idx_dof1=1,p%dof_node + + m%elf(idx_dof1,i) = 0.0_BDKi + DO idx_qp = 1,p%nqp ! dot_product( m%qp%Fc (idx_dof1,:,nelem), p%QPtw_ShpDer( :,i)) + m%elf(idx_dof1,i) = m%elf(idx_dof1,i) - m%qp%Fc(idx_dof1,idx_qp,nelem)*p%QPtw_ShpDer(idx_qp,i) + END DO + DO idx_qp = 1,p%nqp ! dot_product( m%qp%Ftemp(idx_dof1,:,nelem) + m%qp%Fi(idx_dof1,:,nelem), p%QPtw_Shp_Jac(:,i,nelem)) + m%elf(idx_dof1,i) = m%elf(idx_dof1,i) - (m%qp%Ftemp(idx_dof1,idx_qp,nelem) + m%qp%Fi(idx_dof1,idx_qp,nelem))*p%QPtw_Shp_Jac(idx_qp,i,nelem) + END DO + + ! Elastic force internal + m%EFint(idx_dof1,i,nelem) = 0.0_BDKi + DO idx_qp = 1,p%nqp ! dot_product(m%qp%Fc(idx_dof1,:,nelem), p%QPtw_ShpDer( :,i)) + m%EFint(idx_dof1,i,nelem) = m%EFint(idx_dof1,i,nelem) - m%qp%Fc(idx_dof1,idx_qp,nelem)*p%QPtw_ShpDer(idx_qp,i) ENDDO - ENDDO +! m%EFint(idx_dof1,i,nelem) = m%EFint(idx_dof1,i,nelem) - dot_product(m%qp%Fd(idx_dof1,:,nelem), p%QPtw_Shp_Jac(:,i,nelem)) + ENDDO ENDDO + RETURN END SUBROUTINE BD_ElementMatrixGA2 @@ -4082,15 +4396,15 @@ END SUBROUTINE BD_CalcIC !----------------------------------------------------------------------------------------------------------------------------------- !! Routine for computing outputs, used in both loose and tight coupling. -SUBROUTINE BD_InitAcc( u, p, x, OtherState, m, ErrStat, ErrMsg ) +SUBROUTINE BD_InitAcc( u, p, x, m, qdotdot, ErrStat, ErrMsg ) - TYPE(BD_InputType), INTENT(IN ) :: u !< Inputs at t (in BD coordinates) - TYPE(BD_ParameterType), INTENT(IN ) :: p !< Parameters - TYPE(BD_ContinuousStateType), INTENT(IN ) :: x !< Continuous states at t - TYPE(BD_OtherStateType), INTENT(INOUT) :: OtherState !< Other states at t - TYPE(BD_MiscVarType), INTENT(INOUT) :: m !< Misc/optimization variables - INTEGER(IntKi), INTENT( OUT) :: ErrStat !< Error status of the operation - CHARACTER(*), INTENT( OUT) :: ErrMsg !< Error message if ErrStat /= ErrID_None + TYPE(BD_InputType), INTENT(IN ) :: u !< Inputs at t (in BD coordinates) + TYPE(BD_ParameterType), INTENT(IN ) :: p !< Parameters + TYPE(BD_ContinuousStateType), INTENT(IN ) :: x !< Continuous states at t + TYPE(BD_MiscVarType), INTENT(INOUT) :: m !< Misc/optimization variables + REAL(BDKi), INTENT( OUT) :: qdotdot(:,:) !< accelerations + INTEGER(IntKi), INTENT( OUT) :: ErrStat !< Error status of the operation + CHARACTER(*), INTENT( OUT) :: ErrMsg !< Error message if ErrStat /= ErrID_None INTEGER(IntKi) :: ErrStat2 ! Temporary Error status CHARACTER(ErrMsgLen) :: ErrMsg2 ! Temporary Error message @@ -4110,13 +4424,10 @@ SUBROUTINE BD_InitAcc( u, p, x, OtherState, m, ErrStat, ErrMsg ) CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) ! set accelerations with inputs from the root and BD_CalcForceAcc solution - OtherState%Acc(1:3,1 ) = u%RootMotion%TranslationAcc(1:3,1) - OtherState%Acc(4:6,1 ) = u%RootMotion%RotationAcc(1:3,1) - OtherState%Acc( :,2:) = m%RHS(:,2:) + qdotdot(1:3,1 ) = u%RootMotion%TranslationAcc(:,1) + qdotdot(4:6,1 ) = u%RootMotion%RotationAcc(:,1) + qdotdot( :,2:) = m%RHS(:,2:) - ! Initialize Xcc equal to Acc: - OtherState%Xcc(:,:) = OtherState%Acc(:,:) - return END SUBROUTINE BD_InitAcc @@ -4159,7 +4470,7 @@ SUBROUTINE BD_CalcForceAcc( x_tmp, u, p, m, ErrStat, ErrMsg ) CALL BD_ElementMatrixAcc( nelem, p, x_tmp, m ) CALL BD_AssembleStiffK(nelem,p,m%elm, m%MassM) - CALL BD_AssembleRHS(nelem,p,m%elf(:,:,nelem), m%RHS) + CALL BD_AssembleRHS(nelem,p,m%elf, m%RHS) ENDDO ! ending of old BD_GenerateDynamicElementAcc diff --git a/modules-local/beamdyn/src/BeamDyn_IO.f90 b/modules-local/beamdyn/src/BeamDyn_IO.f90 index 05930615d1..316662aa4e 100644 --- a/modules-local/beamdyn/src/BeamDyn_IO.f90 +++ b/modules-local/beamdyn/src/BeamDyn_IO.f90 @@ -1,7 +1,7 @@ !********************************************************************************************************************************** ! LICENSING ! Copyright (C) 2015-2016 National Renewable Energy Laboratory -! Copyright (C) 2016-2017 Envision Energy USA, LTD +! Copyright (C) 2016-2017 Envision Energy USA, LTD ! ! Licensed under the Apache License, Version 2.0 (the "License"); ! you may not use this file except in compliance with the License. @@ -25,7 +25,7 @@ MODULE BeamDyn_IO IMPLICIT NONE - TYPE(ProgDesc), PARAMETER:: BeamDyn_Ver = ProgDesc('BeamDyn', 'v2.00.00','9-May-2017') + TYPE(ProgDesc), PARAMETER:: BeamDyn_Ver = ProgDesc('BeamDyn', 'v2.00.00','14-Aug-2017') @@ -1559,23 +1559,21 @@ SUBROUTINE Calc_WriteOutput( p, AllOuts, y, m, ErrStat, ErrMsg ) !------------------------------------ ! outputs on the nodes - do beta=1,p%NNodeOuts + DO beta=1,p%NNodeOuts j=p%OutNd(beta) j_BldMotion = p%NdIndx(j) !------------------------------------ ! Sectional force resultants at Node 1 expressed in l, given in N - !FIXME: N#Mxl & N#Myl are known to be incorrect! See https://wind.nrel.gov/forum/wind/viewtopic.php?f=3&t=1622 - temp_vec = MATMUL(y%BldMotion%Orientation(:,:,j_BldMotion), m%BldInternalForce(1:3,j_BldMotion)) + temp_vec = MATMUL(y%BldMotion%Orientation(:,:,j_BldMotion), m%BldInternalForceFE(1:3,j)) AllOuts( NFl( beta,1 ) ) = temp_vec(1) AllOuts( NFl( beta,2 ) ) = temp_vec(2) AllOuts( NFl( beta,3 ) ) = temp_vec(3) !------------------------------------ ! Sectional moment resultants at Node 1 expressed in l, given in N-m - !FIXME: N#Mxl & N#Myl are known to be incorrect! See https://wind.nrel.gov/forum/wind/viewtopic.php?f=3&t=1622 - temp_vec = MATMUL(y%BldMotion%Orientation(:,:,j_BldMotion), m%BldInternalForce(4:6,j_BldMotion)) + temp_vec = MATMUL(y%BldMotion%Orientation(:,:,j_BldMotion), m%BldInternalForceFE(4:6,j)) AllOuts( NMl( beta,1 ) ) = temp_vec(1) AllOuts( NMl( beta,2 ) ) = temp_vec(2) AllOuts( NMl( beta,3 ) ) = temp_vec(3) @@ -1585,7 +1583,7 @@ SUBROUTINE Calc_WriteOutput( p, AllOuts, y, m, ErrStat, ErrMsg ) ! Sectional translational deflection (relative to the undeflected position) expressed in r d = y%BldMotion%TranslationDisp(:, j_BldMotion) - m%u2%RootMotion%TranslationDisp(:,1) d_ref = y%BldMotion%Position( :, j_BldMotion) - m%u2%RootMotion%Position( :,1) - temp_vec2 = d + d_ref - matmul( RootRelOrient, d_ref ) ! tip displacement + temp_vec2 = d + d_ref - MATMUL( RootRelOrient, d_ref ) ! tip displacement temp_vec = MATMUL(m%u2%RootMotion%Orientation(:,:,1),temp_vec2) AllOuts( NTDr( beta,1 ) ) = temp_vec(1) diff --git a/modules-local/beamdyn/src/BeamDyn_Subs.f90 b/modules-local/beamdyn/src/BeamDyn_Subs.f90 index 631cb8aa56..86d3826c90 100644 --- a/modules-local/beamdyn/src/BeamDyn_Subs.f90 +++ b/modules-local/beamdyn/src/BeamDyn_Subs.f90 @@ -315,7 +315,9 @@ SUBROUTINE BD_CrvCompose( rr, pp, qq, flag) END SUBROUTINE BD_CrvCompose !----------------------------------------------------------------------------------------------------------------------------------- !> This subroutine computes the CRV parameters given -!! the rotation matrix +!! the rotation matrix. +!! NOTE: values for the DCM passed into this routine must be accurate at least ~12 digits of precision on the diagonal for this +!! routine to work properly. The received DMC is assumed to have all the properties of a well formed accurate DCM. SUBROUTINE BD_CrvExtractCrv(Rr,cc) REAL(BDKi), INTENT(IN ):: Rr(3,3) !< Rotation Matrix @@ -628,27 +630,31 @@ SUBROUTINE Set_BldMotion_Mesh(p, u, x, m, y) ! set positions and velocities (not accelerations) call Set_BldMotion_NoAcc(p, u, x, m, y) - - ! now set the accelerations: - - ! The first node on the mesh is just the root location: - y%BldMotion%TranslationAcc(:,1) = u%RootMotion%TranslationAcc(:,1) - y%BldMotion%RotationAcc(:,1) = u%RootMotion%RotationAcc(:,1) + + ! Only need this bit for dynamic cases + IF (p%analysis_type == BD_DYNAMIC_ANALYSIS) THEN + + ! now set the accelerations: + + ! The first node on the mesh is just the root location: + y%BldMotion%TranslationAcc(:,1) = u%RootMotion%TranslationAcc(:,1) + y%BldMotion%RotationAcc(:,1) = u%RootMotion%RotationAcc(:,1) - j_start = 2 ! we'll skip the first node on the first element; otherwise we will start at the first node on the other elements - DO i=1,p%elem_total - DO j=j_start,p%nodes_per_elem - temp_id = (i-1)*(p%nodes_per_elem-1)+j - temp_id2= (i-1)*p%nodes_per_elem+j + j_start = 2 ! we'll skip the first node on the first element; otherwise we will start at the first node on the other elements + DO i=1,p%elem_total + DO j=j_start,p%nodes_per_elem + temp_id = (i-1)*(p%nodes_per_elem-1)+j + temp_id2= (i-1)*p%nodes_per_elem+j - y%BldMotion%TranslationAcc(1:3,temp_id2) = MATMUL(p%GlbRot, m%RHS(1:3,temp_id) ) + y%BldMotion%TranslationAcc(1:3,temp_id2) = MATMUL(p%GlbRot, m%RHS(1:3,temp_id) ) - y%BldMotion%RotationAcc(1:3,temp_id2) = MATMUL(p%GlbRot, m%RHS(4:6,temp_id) ) - ENDDO - j_start = 1 - ENDDO - + y%BldMotion%RotationAcc(1:3,temp_id2) = MATMUL(p%GlbRot, m%RHS(4:6,temp_id) ) + ENDDO + j_start = 1 + ENDDO + END IF + END SUBROUTINE Set_BldMotion_Mesh !----------------------------------------------------------------------------------------------------------------------------------- !> This subroutine finds the (initial) nodal position and curvature based on a relative position, eta, along the blade. diff --git a/modules-local/beamdyn/src/Driver_Beam.f90 b/modules-local/beamdyn/src/Driver_Beam.f90 index df4b821ca5..221c37509a 100644 --- a/modules-local/beamdyn/src/Driver_Beam.f90 +++ b/modules-local/beamdyn/src/Driver_Beam.f90 @@ -68,7 +68,7 @@ PROGRAM BeamDyn_Driver_Program - TYPE(ProgDesc), PARAMETER :: version = ProgDesc( 'BeamDyn Driver', 'v2.00.00', '9-May-2017' ) ! The version number of this program. + TYPE(ProgDesc), PARAMETER :: version = ProgDesc( 'BeamDyn Driver', 'v2.00.00', '14-Aug-2017' ) ! The version number of this program. ! ------------------------------------------------------------------------- @@ -84,7 +84,7 @@ PROGRAM BeamDyn_Driver_Program ! Display the copyright notice CALL DispCopyrightLicense( version ) ! Tell our users what they're running - CALL WrScr( ' Running '//GetNVD( version )//NewLine//' linked with '//TRIM( GetNVD( NWTC_Ver ))//NewLine ) + CALL WrScr( ' Running '//GetNVD( version )//' compiled on '//__DATE__//NewLine//' linked with '//TRIM( GetNVD( NWTC_Ver ))//NewLine ) ! ------------------------------------------------------------------------- ! Initialization of glue-code time-step variables diff --git a/modules-local/beamdyn/src/Registry_BeamDyn.txt b/modules-local/beamdyn/src/Registry_BeamDyn.txt index 92328e2747..eb9fef671e 100644 --- a/modules-local/beamdyn/src/Registry_BeamDyn.txt +++ b/modules-local/beamdyn/src/Registry_BeamDyn.txt @@ -110,8 +110,8 @@ typedef ^ ConstraintStateType ReKi DummyConstrState - - - "A variable, # at a given time, etc.) #vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv #the following are BDKi = R8Ki -typedef ^ OtherStateType R8Ki acc {:}{:} - - "Acceleration in GA2" -typedef ^ OtherStateType ^ xcc {:}{:} - - "Algorithm acceleration in GA2" +typedef ^ OtherStateType R8Ki acc {:}{:} - - "Acceleration (dqdtdt)" +typedef ^ OtherStateType ^ xcc {:}{:} - - "Algorithm acceleration in GA2: (1-alpha_m)*xcc_(n+1) = (1-alpha_f)*Acc_(n+1) + alpha_f*Acc_n - alpha_m*xcc_n" #end of BDKi-type variables #^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ typedef ^ OtherStateType Logical InitAcc - - - "flag to determine if accerlerations have been initialized in updateStates" @@ -194,6 +194,12 @@ typedef ^ ParameterType ReKi pitchC - - - typedef ^ ParameterType ReKi torqM {2}{2} - - "Pitch actuator matrix: (I-hA)^-1" - typedef ^ ParameterType qpParam qp - - - "Quadrature point info that does not change during simulation" - typedef ^ ParameterType IntKi qp_indx_offset - - - "Offset for computing index of the quadrature arrays (gauss skips the first [end-point] node)" - +# .... arrays for optimization ........................................................................................................ +typedef ^ ParameterType R8Ki QPtw_Shp_Shp_Jac {:}{:}{:}{:} - - "optimization variable: QPtw_Shp_Shp_Jac(idx_qp,i,j,nelem) = p%Shp(i,idx_qp)*p%Shp(j,idx_qp)*p%QPtWeight(idx_qp)*p%Jacobian(idx_qp,nelem)" - +typedef ^ ParameterType ^ QPtw_Shp_ShpDer {:}{:}{:} - - "optimization variable: QPtw_Shp_ShpDer(idx_qp,i,j) = p%Shp(i,idx_qp)*p%ShpDer(j,idx_qp)*p%QPtWeight(idx_qp)" - +typedef ^ ParameterType ^ QPtw_ShpDer_ShpDer_Jac {:}{:}{:}{:} - - "optimization variable: QPtw_ShpDer_ShpDer_Jac(idx_qp,i,j,nelem) = p%ShpDer(i,idx_qp)*p%ShpDer(j,idx_qp)*p%QPtWeight(idx_qp)/p%Jacobian(idx_qp,nelem)" - +typedef ^ ParameterType ^ QPtw_Shp_Jac {:}{:}{:} - - "optimization variable: QPtw_Shp_Jac(idx_qp,i,nelem) = p%Shp(i,idx_qp)*p%QPtWeight(idx_qp)*p%Jacobian(idx_qp,nelem)" - +typedef ^ ParameterType ^ QPtw_ShpDer {:}{:} - - "optimization variable: QPtw_ShpDer(idx_qp,i) = p%ShpDer(i,idx_qp)*p%QPtWeight(idx_qp)" - # ..... Inputs # .................................................................................................................... @@ -283,7 +289,8 @@ typedef ^ MiscVarType EqMotionQP qp - - - "Quadra #vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv #the following are BDKi = R8Ki typedef ^ MiscVarType R8Ki Nrrr ::: - - "Rotation parameters relative to root node, from Nuuu (at GLL points)" - -typedef ^ MiscVarType ^ elf {:}{:}{:} - - "" - +typedef ^ MiscVarType ^ elf {:}{:} - - "" - +typedef ^ MiscVarType ^ EFint {:}{:}{:} - - "Elastic forces internal to blade. For output calculations only." - typedef ^ MiscVarType ^ elk {:}{:}{:}{:} - - "" - typedef ^ MiscVarType ^ elg {:}{:}{:}{:} - - "" - typedef ^ MiscVarType ^ elm {:}{:}{:}{:} - - "" - @@ -293,7 +300,7 @@ typedef ^ MiscVarType ^ StifK {:}{:}{:}{:} - - "S typedef ^ MiscVarType ^ MassM {:}{:}{:}{:} - - "Mass Matrix" - typedef ^ MiscVarType ^ DampG {:}{:}{:}{:} - - "Damping Matrix" - typedef ^ MiscVarType ^ RHS {:}{:} - - "Right-hand-side vector" - -typedef ^ MiscVarType ^ BldInternalForce {:}{:} - - "temporary force array used in dynamic calculations" - +typedef ^ MiscVarType ^ BldInternalForceFE {:}{:} - - "Force/Moment array for internal force calculations at FE" - typedef ^ MiscVarType ^ Solution {:}{:} - - "Result from LAPACK solve (X from A*X = B solve)" - # arrays for lapack routines typedef ^ MiscVarType ^ LP_StifK {:}{:} - - "Stiffness Matrix" - @@ -308,4 +315,3 @@ typedef ^ MiscVarType IntKi LP_indx {:} - - "I typedef ^ MiscVarType BD_InputType u - - - "Inputs converted to the internal BD coordinate system" - typedef ^ MiscVarType BD_InputType u2 - - - "Inputs in the FAST coordinate system, possibly modified by pitch actuator" - -