diff --git a/src/vpmCommon/supElTypeModule.f90 b/src/vpmCommon/supElTypeModule.f90 index e865804b..0a4f863a 100644 --- a/src/vpmCommon/supElTypeModule.f90 +++ b/src/vpmCommon/supElTypeModule.f90 @@ -144,6 +144,8 @@ module SupElTypeModule real(dp) :: dmpScl(2) !< Current and previous scaling of mDmp0 and kDmp0 integer :: dmpSclIdx !< Index to structural damping scaling function + real(dp) :: massScl(2) !< Current and previous mass scaling factors + integer :: massSclIdx !< Index to mass scaling function real(dp) :: stifScl(2) !< Current and previous stiffness scaling factors integer :: stifSclIdx !< Index to stiffness scaling function @@ -356,6 +358,7 @@ subroutine WriteSupElType (sup,io,complexity) write(io,3) 'mDmpFactor =', sup%mDmpFactor write(io,3) 'kDmpFactor =', sup%kDmpFactor write(io,*) 'dmpScale(id) =', sup%dmpSclIdx + write(io,*) 'massScale(id) =', sup%massSclIdx write(io,*) 'stifScale(id) =', sup%stifSclIdx write(io,*) 'savePos =', sup%savePos write(io,*) 'saveVar =', sup%saveVar @@ -666,8 +669,10 @@ subroutine NullifySupEl (sup) sup%mDmp0 = 0.0_dp sup%kDmp0 = 0.0_dp sup%dmpScl = 1.0_dp + sup%massScl = 1.0_dp sup%stifScl = 1.0_dp sup%dmpSclIdx = 0 + sup%massSclIdx = 0 sup%stifSclIdx = 0 sup%savePos = .false. sup%saveVar = .false. @@ -1273,6 +1278,7 @@ module SupElNamelistModule integer :: recoveryFlag !< Flag for enabling stress recovery integer :: BC(maxGenDofs_p) !< Boundary conditions for component modes integer :: strDmpEngineId !< Base ID for structural damping function + integer :: massEngineId !< Base ID for mass scaling function integer :: stiffEngineId !< Base ID for stiffness scaling function integer :: savePos !< Flag indicating whether position should be saved @@ -1307,7 +1313,7 @@ module SupElNamelistModule & refTriad2Id, offset2, & & refTriad3Id, offset3, & & stressStiffFlag, massCorrFlag, recoveryFlag, & - & stiffEngineId, stiffScale, massScale, & + & stiffEngineId, massEngineId, stiffScale, massScale, & & dragParams, slamParams, & & strDmpEngineId, alpha1, alpha2, alpha3, alpha4, BC, & & savePos, saveVar @@ -1343,7 +1349,8 @@ subroutine read_SUP_EL (infp,stat) stressStiffFlag=-1; massCorrFlag=-1; recoveryFlag=-1 stiffScale=1.0_dp; massScale=1.0_dp; dragParams=0.0_dp; slamParams=0.0_dp alpha1=0.0_dp; alpha2=0.0_dp; alpha3=0.0_dp; alpha4=0.0_dp - stiffEngineId=0; strDmpEngineId=0; BC=1; savePos=0; saveVar=0 + massEngineId=0; stiffEngineId=0; strDmpEngineId=0 + BC=1; savePos=0; saveVar=0 read(infp,nml=SUP_EL,iostat=stat) diff --git a/src/vpmSolver/addInSysModule.f90 b/src/vpmSolver/addInSysModule.f90 index 6f932240..44332692 100644 --- a/src/vpmSolver/addInSysModule.f90 +++ b/src/vpmSolver/addInSysModule.f90 @@ -40,7 +40,7 @@ module AddInSysModule !> @param Rhs Right-hand-side force vector !> !> @details The Newton matrix is a linear combination of the mass-, damping- - !> and stiffness matrices, Nmat = scaleM*M + scaleC*C + scaleK*K. + !> and stiffness matrices, Nmat = scaleM*M + scaleC*C + scaleK*K. !> If @a Rhs is provided, contributions due to prescribed motions, !> if any, will be added. !> @@ -142,6 +142,7 @@ subroutine BuildNewtonMat (Nmat, scaleM, scaleC, scaleK, mech, sam, iter, & do i = 1, size(mech%sups) if ( reComputeSupMat .or. mech%sups(i)%stressStiffFlag(1) == 1 .or. & + abs(mech%sups(i)%massScl(1)-mech%sups(i)%massScl(2)) > eps_p .or. & abs(mech%sups(i)%stifScl(1)-mech%sups(i)%stifScl(2)) > eps_p ) then call BuildSupNewtonMat (newTanStiff, scaleM, scaleC, scaleK, & @@ -353,7 +354,8 @@ subroutine BuildMassMat (Mmat, mech, sam, ierr) do i = 1, size(mech%sups) if (associated(mech%sups(i)%Mmat)) then - call addInSupMat (mech%sups(i)%Mmat, Mmat, mech%sups(i), sam, ierr) + call addInSupMat (mech%sups(i)%Mmat, Mmat, mech%sups(i), sam, ierr, & + & scale=mech%sups(i)%massScl(1)) if (ierr < 0) goto 900 end if if (associated(mech%sups(i)%Mamat)) then diff --git a/src/vpmSolver/initiateSupElTypeModule.f90 b/src/vpmSolver/initiateSupElTypeModule.f90 index 54319797..7bc1154a 100644 --- a/src/vpmSolver/initiateSupElTypeModule.f90 +++ b/src/vpmSolver/initiateSupElTypeModule.f90 @@ -551,7 +551,7 @@ subroutine ReadSupEls (infp,env,triads,sups,err) end if end if - scaledStrDmp = scaledStructDamp .or. stiffEngineId > 0 + scaledStrDmp = scaledStructDamp .or. stiffEngineId>0 .or. massEngineId>0 if (scaledStrDmp .or. .not.haveDamping(0)) then !! Apply stiffness and/or mass scaling factors @@ -610,6 +610,19 @@ subroutine ReadSupEls (infp,env,triads,sups,err) sups(idIn)%genDOFs%alpha2 = alpha4(1:numGenDOFs) end if + if (massEngineId > 0) then + + !! Time-dependent mass scaling + if (haveDamping(2)) then + err = err - 1 + call ReportInputError('SUP_EL',idIn,sups(idIn)%id, & + 'Stiffness-proportional damping can not be combined with '// & + 'time-dependent mass scaling.') + else + sups(idIn)%massSclIdx = massEngineId + end if + + end if if (stiffEngineId > 0) then !! Time-dependent stiffness scaling @@ -903,6 +916,14 @@ subroutine InitiateSupEls2 (sam,sys,sups,masses,engines,err) & 'Non-existing engine referred') end if end if + if (sups(idIn)%massSclIdx > 0) then + if (.not. associated(GetPtrToId(engines,sups(idIn)%massSclIdx, & + & index=sups(idIn)%massSclIdx))) then + lerr = lerr - 1 + call ReportInputError('SUP_EL',idIn,sups(idIn)%id, & + & 'Non-existing engine referred') + end if + end if if (sups(idIn)%stifSclIdx > 0) then if (.not. associated(GetPtrToId(engines,sups(idIn)%stifSclIdx, & & index=sups(idIn)%stifSclIdx))) then diff --git a/src/vpmSolver/supElRoutinesModule.f90 b/src/vpmSolver/supElRoutinesModule.f90 index 119c7a8a..f6f22ed4 100644 --- a/src/vpmSolver/supElRoutinesModule.f90 +++ b/src/vpmSolver/supElRoutinesModule.f90 @@ -298,7 +298,8 @@ subroutine UpdateSupEls (sups,supLoads,env,beta,gamma,time,timeStep, & if (associated(sups(i)%Mmat)) then !! Calculate the inertia forces in the superelement, Fi = M*uldd - call scaledMatmul (m,n,1.0_dp,sups(i)%Mmat,sups(i)%uldd,sups(i)%Fi) + call scaledMatmul (m,n,sups(i)%massScl(1), & + & sups(i)%Mmat,sups(i)%uldd,sups(i)%Fi) else call DCOPY (m,0.0_dp,0,sups(i)%Fi(1),1) end if @@ -335,7 +336,7 @@ subroutine UpdateSupEls (sups,supLoads,env,beta,gamma,time,timeStep, & if (associated(sups(i)%fg)) then !! Calculate the gravitational forces which is an external force - call scaledMatmul (m,3,sups(i)%stifScl(1),sups(i)%fg, & + call scaledMatmul (m,3,sups(i)%massScl(1),sups(i)%fg, & & matmul(env%gravity,sups(i)%supTr(:,1:3)),sups(i)%Q) else call DCOPY (m,0.0_dp,0,sups(i)%Q(1),1) @@ -698,7 +699,8 @@ end subroutine UpdateSupElLoad !> @param[out] ierr Error flag !> !> @details Effective with time-dependent structural damping coefficients. - !> This subroutine also updates the time-dependent stiffness scaling factor. + !> This subroutine also updates the time-dependent stiffness- and mass + !> scaling factors. !> !> @callgraph @callergraph !> @@ -739,6 +741,10 @@ subroutine UpdateSupElDamping (sups,engs,alpha,ierr) else newDmp = .false. end if + if (sups(i)%massSclIdx > 0) then + sups(i)%massScl(2) = sups(i)%massScl(1) + sups(i)%massScl(1) = EngineValue(engs(sups(i)%massSclIdx),ierr) + end if if (sups(i)%stifSclIdx > 0) then sups(i)%stifScl(2) = sups(i)%stifScl(1) sups(i)%stifScl(1) = EngineValue(engs(sups(i)%stifSclIdx),ierr) @@ -1271,7 +1277,7 @@ subroutine BuildSupNewtonMat (newTangent,scaleM,scaleC,scaleK,sup,ierr) !! Structural mass matrix contribution call DCOPY (n,sup%Mmat(1,1),1,sup%Nmat(1,1),1) - call DSCAL (n,scaleM,sup%Nmat(1,1),1) + call DSCAL (n,scaleM*sup%massScl(1),sup%Nmat(1,1),1) else call DCOPY (n,0.0_dp,0,sup%Nmat(1,1),1)