Skip to content

Commit

Permalink
Merge pull request #2561 from luwang00/f/YawFrctUpdate
Browse files Browse the repository at this point in the history 
ED: Extended yaw-friction modeling
  • Loading branch information
@andrew-platt
andrew-platt authored Dec 18, 2024
2 parents b996397 + 09b1e42 commit 4377c13
Show file tree
Hide file tree
Showing 9 changed files with 249 additions and 72 deletions.
Binary file modified docs/source/user/elastodyn/figs/YawFrictionModel.jpg
View file
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
20 changes: 16 additions & 4 deletions docs/source/user/elastodyn/input.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -238,13 +238,25 @@ Rotor-Teeter
Yaw-Friction
~~~~~~~~~~~~

**YawFrctMod** - Yaw-friction model {0: none, 1: friction without Fz term at the yaw bearing, 2: friction includes Fz term at yaw bearing, 3: user defined model}
**YawFrctMod** - Yaw-friction model {0: none, 1: friction independent of yaw-bearing force and bending moment, 2: friction with Coulomb terms depending on yaw-bearing force and bending moment, 3: user defined model}

**M_CSmax** - Maximum Coulomb friction torque (N-m)[mu_s*D_eff when YawFrctMod=1 and Fz*mu_s*D_eff when YawFrctMod=2]
**M_CSmax** - Maximum static Coulomb friction torque (N-m) [M_CSmax when YawFrctMod=1; -min(0,Fz)*M_CSmax when YawFrctMod=2]

**M_CD** - Dynamic friction moment at null yaw rate (N-m) [mu_d*D_eff when YawFrctMod=1 and Fz*mu_d*D_eff when YawFrctMod=2]
**M_FCSmax** - Maximum static Coulomb friction torque proportional to yaw bearing shear force (N-m) [sqrt(Fx^2+Fy^2)*M_FCSmax; only used when YawFrctMod=2]

**sig_v** - Viscous friction coefficient (N-m/(rad/s))
**M_MCSmax** - Maximum static Coulomb friction torque proportional to yaw bearing bending moment (N-m) [sqrt(Mx^2+My^2)*M_MCSmax; only used when YawFrctMod=2]

**M_CD** - Dynamic Coulomb friction moment (N-m) [M_CD when YawFrctMod=1; -min(0,Fz)*M_CD when YawFrctMod=2]

**M_FCD** - Dynamic Coulomb friction moment proportional to yaw bearing shear force (N-m) [sqrt(Fx^2+Fy^2)*M_FCD; only used when YawFrctMod=2]

**M_MCD** - Dynamic Coulomb friction moment proportional to yaw bearing bending moment (N-m) [sqrt(Mx^2+My^2)*M_MCD; only used when YawFrctMod=2]

**sig_v** - Linear viscous friction coefficient (N-m/(rad/s))

**sig_v2** - Quadratic viscous friction coefficient (N-m/(rad/s)^2)

**OmgCut** - Yaw angular velocity cutoff below which viscous friction is linearized (rad/s)


Drivetrain
Expand Down
65 changes: 47 additions & 18 deletions docs/source/user/elastodyn/theory.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -153,16 +153,12 @@ The total moment on the given degree of freedom is:
.. math:: :label: TFTotTorque

Q = Q_\text{lin} + Q_\text{stop,spr} + Q_\text{stop,dmp}






.. _ed_yawfriction_theory:

Yaw-friction
------------
Yaw-friction model
------------------
A yaw-friction model is implemented in ElastoDyn based on a Coulomb-viscous approach.
The yaw-friction moment as a function of yaw rate (:math:`\omega`) is shown below in :numref:`figYawFriction`

Expand All @@ -172,29 +168,62 @@ The yaw-friction moment as a function of yaw rate (:math:`\omega`) is shown belo

Yaw-friction model

The yaw-friction torque :math:`M_f` can be calcualated as follows.
When ``YawFrctMod`` = 1, the maximum static or dynamic Coulomb friction does not depend on the external load on the yaw bearing. The yaw-friction torque :math:`M_f` can be calculated as follows.
If :math:`\omega\neq0`, we have dynamic friction of the form

.. math::
M_f = -(\mu_d\bar{D})\cdot\textrm{sign}(\omega) - M_{f,vis},
where :math:`\bar{D}` is the effective yaw-bearing diameter and :math:`\mu_d` is the dynamic Coulomb friction coefficient. Their product, :math:`\mu_d\bar{D}`, is specified in the input file through ``M_CD``. The first term on the right-hand side is the dynamic Coulomb friction.
The viscous friction, :math:`M_{f,vis}`, is of the form

.. math::
M_{f,vis} = \sigma_v\omega + \sigma_{v2}\omega\left|\omega\right|\qquad\qquad\text{if}~\left|\omega\right|\ge\omega_c,
if :math:`\omega\neq0`:
or

.. math::
M_f = \mu_d\bar{D}\times\text{min}(0,F_z)\times\text{sign}(\omega) - \sigma_v\times\omega
if :math:`\omega=0` and :math:`\dot{\omega}\neq 0`:
M_{f,vis} = (\sigma_v + \sigma_{v2}\omega_c)\omega\qquad\qquad\text{if}~\left|\omega\right|\le\omega_c,
where :math:`\sigma_v` and :math:`\sigma_{v2}` are the linear and quadratic viscous friction coefficients and :math:`\omega_c` is the cutoff yaw rate below which viscous friction is linearized. Setting :math:`\omega_c=0` disables the linearization of viscous friction.

If :math:`\omega=0` and :math:`\dot{\omega}\neq 0`, we have a slightly modified dynamic Coulomb friction of the form

.. math::
M_f = \mu_d\bar{D}\times\text{min}(0,F_z)\times\text{sign}(\dot{\omega})
M_f = -\textrm{min}\!\left(\mu_d\bar{D},\left|M_z\right|\right)\cdot\textrm{sign}(M_z),
if :math:`\omega=0` and :math:`\dot{\omega}=0`:
where :math:`M_z` is the external yaw torque.
If :math:`\omega=0` and :math:`\dot{\omega}=0`, we have static Coulomb friction of the form

.. math::
M_f = -\text{min}(\mu_s\bar{D}\times|\text{min}(0,F_z)|,|M_z|)\times\text{sign}(M_z)
M_f = -\textrm{min}\!\left(\mu_s\bar{D},\left|M_z\right|\right)\cdot\textrm{sign}(M_z),
where :math:`\mu_s` is the static Coulomb friction coefficient. The product :math:`\mu_s\bar{D}` is specified in the input file through ``M_CSmax``.

where :math:`\bar{D}` is the effective yaw-bearing race diameter, :math:`\mu_d` is the dynamic friction coefficient, :math:`\mu_s` is the static friction coefficient, :math:`F_z` is effective axial load on yaw-bearing, :math:`M-z` is the external torque on yaw-bearing.
The static 'stiction' (where the static contribution exceeds the dynamic Coulomb friction) is only applied if both the yaw rotational velocity and acceleration at the current time-step are zero.
The static portion of the friction is omitted if the rotational acceleration is not null, (sign(0) is taken as 1).
This is to account for the fact that a 'warm' joint may not feel stiction when crossing through zero velocity in a dynamic sense :cite:`ed-hammam2023`

When ``YawFrctMod`` = 2, the maximum static or dynamic Coulomb friction depends on the external load on the yaw bearing, with proportional contributions from :math:`\left|F_z\right|`, the magnitude of the bearing axial load, if :math:`F_z<0`, from the bearing shear force magnitude, :math:`\sqrt{F_x^2+F_y^2}`, and from the bearing bending moment magnitude, :math:`\sqrt{M_x^2+M_y^2}`.
If :math:`\omega\neq0`, we have dynamic friction of the form

.. math::
M_f = \left(\mu_d\bar{D}\cdot\textrm{min}\!\left(0,F_z\right)-\mu_{df}\bar{D}\sqrt{F_x^2+F_y^2}-\mu_{dm}\sqrt{M_x^2+M_y^2}\right)\cdot\textrm{sign}(\omega) - M_{f,vis},
where :math:`M_{f,vis}` is defined in the same way as when ``YawFrctMod`` = 1. The product :math:`\mu_{df}\bar{D}` and :math:`\mu_{dm}` are specified in the input file through ``M_FCD`` and ``M_MCD``, respectively.
If :math:`\omega=0` and :math:`\dot{\omega}\neq 0`, we have a modified dynamic Coulomb friction of the form

.. math::
M_f = -\textrm{min}\!\left(\mu_d\bar{D}\left|\textrm{min}(0,F_z)\right| + \mu_{df}\bar{D}\sqrt{F_x^2+F_y^2} + \mu_{dm}\sqrt{M_x^2+M_y^2},\left|M_z\right|\right)\cdot\textrm{sign}(M_z).
If :math:`\omega=0` and :math:`\dot{\omega}=0`, we have static Coulomb friction of the form

.. math::
M_f = -\textrm{min}\!\left(\mu_s\bar{D}\left|\textrm{min}(0,F_z)\right| + \mu_{sf}\bar{D}\sqrt{F_x^2+F_y^2} + \mu_{sm}\sqrt{M_x^2+M_y^2},\left|M_z\right|\right)\cdot\textrm{sign}(M_z),
where the product :math:`\mu_{sf}\bar{D}` and :math:`\mu_{sm}` are specified in the input file through ``M_FCSmax`` and ``M_MCSmax``, respectively.

The static 'stiction' (where the static contribution exceeds the dynamic Coulomb friction) is only applied if both the yaw rotational velocity and acceleration at the current time-step are zero.
The static portion of the friction is omitted if the rotational acceleration is not null.
This is to account for the fact that a 'warm' joint may not feel stiction when crossing through zero velocity in a dynamic sense :cite:`ed-hammam2023`.
When :math:`\omega=0`, the yaw-bearing static or dynamic friction is formulated such that the frictional resistance opposes the external applied moment, :math:`M_z`, without overcoming it.

.. _ed_dev_notes:

Expand Down
5 changes: 3 additions & 2 deletions modules/elastodyn/src/ED_UserSubs.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -101,7 +101,7 @@ SUBROUTINE UserTeet ( TeetDef, TeetRate, ZTime, DirRoot, TeetMom )
RETURN
END SUBROUTINE UserTeet
!=======================================================================
SUBROUTINE UserYawFrict ( ZTime, Fz, Mzz, Omg, OmgDot, DirRoot, YawFriMf )
SUBROUTINE UserYawFrict ( ZTime, F, M, Mzz, Omg, OmgDot, DirRoot, YawFriMf )

! This is a dummy routine for holding the place of a user-specified
! Yaw Friction. Modify this code to create your own device.
Expand All @@ -115,7 +115,8 @@ SUBROUTINE UserYawFrict ( ZTime, Fz, Mzz, Omg, OmgDot, DirRoot, YawFriMf )

! Passed Variables:
REAL(DbKi), INTENT(IN ) :: ZTime ! Current simulation time, sec.
REAL(R8Ki), INTENT(IN ) :: Fz, Mzz ! Yaw Bering normal force (positive if upward) and torque, N and N*m
REAL(ReKi), INTENT(IN ) :: F(3),M(3) ! Yaw bearing force and moment N and N*m
REAL(R8Ki), INTENT(IN ) :: Mzz ! External axial yaw bearing torque N*m
REAL(R8Ki), INTENT(IN ) :: Omg ! Yaw rotational speed, rad/s.
REAL(R8Ki), INTENT(IN ) :: OmgDot ! Yaw rotational acceleration, rad/s^2.

Expand Down
82 changes: 47 additions & 35 deletions modules/elastodyn/src/ElastoDyn.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1247,11 +1247,12 @@ SUBROUTINE ED_CalcOutput( t, u, p, x, xd, z, OtherState, y, m, ErrStat, ErrMsg )
m%AllOuts( YawBrMxp) = DOT_PRODUCT( MomBNcRt, m%CoordSys%b1 )
m%AllOuts( YawBrMyp) = -DOT_PRODUCT( MomBNcRt, m%CoordSys%b3 )
m%AllOuts(YawFriMom) = OtherState%Mfhat*0.001_ReKi !KBF add YawFricMom as an output based on HSSBrTq (kN-m)
m%AllOuts(YawFriMfp) = OtherState%YawFriMfp*0.001_ReKi
m%AllOuts(YawFriMz) = m%YawFriMz*0.001_ReKi
m%FrcONcRt = m%AllOuts( YawBrFzn)*1000_ReKi
m%AllOuts(OmegaYF) = OtherState%OmegaTn*R2D
m%AllOuts(dOmegaYF) = OtherState%OmegaDotTn*R2D
m%AllOuts(YawFriMfp) = OtherState%YawFriMfp*0.001_ReKi
m%AllOuts(YawFriMz) = m%YawFriMz*0.001_ReKi
m%FrcONcRt = (/m%AllOuts( YawBrFxn),m%AllOuts( YawBrFyn),m%AllOuts( YawBrFzn)/) * 1000_ReKi
m%MomONcRt = (/m%AllOuts( YawBrMxn),m%AllOuts( YawBrMyn),m%AllOuts( YawBrMzn)/) * 1000_ReKi
m%AllOuts(OmegaYF) = OtherState%OmegaTn*R2D
m%AllOuts(dOmegaYF) = OtherState%OmegaDotTn*R2D

! Tower Base Loads:

Expand Down Expand Up @@ -1874,8 +1875,7 @@ SUBROUTINE ED_CalcContStateDeriv( t, u, p, x, xd, z, OtherState, m, dxdt, ErrSta
INTEGER(IntKi) :: ErrStat2 ! The error status code
CHARACTER(ErrMsgLen) :: ErrMsg2 ! The error message, if an error occurred
CHARACTER(*), PARAMETER :: RoutineName = 'ED_CalcContStateDeriv'
Real(R8Ki) :: YawFriMz ! Loops through some or all of the DOFs.
Real(R8Ki) :: Fz ! Loops through some or all of the DOFs.
Real(R8Ki) :: YawFriMz ! External loading on yaw bearing not including inertial contributions


! Initialize ErrStat
Expand Down Expand Up @@ -1917,11 +1917,10 @@ SUBROUTINE ED_CalcContStateDeriv( t, u, p, x, xd, z, OtherState, m, dxdt, ErrSta
CALL RFurling( t, p, x%QT(DOF_RFrl), x%QDT(DOF_RFrl), m%RtHS%RFrlMom ) ! Compute moment from rotor-furl springs and dampers, RFrlMom
CALL TFurling( t, p, x%QT(DOF_TFrl), x%QDT(DOF_TFrl), m%RtHS%TFrlMom ) ! Compute moment from tail-furl springs and dampers, TFrlMom
! Compute the yaw friction torque
Fz= m%FrcONcRt !YawBrFzn force from CalcOutput
YawFriMz=DOT_PRODUCT( m%RtHS%MomBNcRtt, m%CoordSys%d2 ) + u%YawMom
m%YawFriMz = YawFriMz

CALL YawFriction( t, p, Fz, YawFriMz, OtherState%OmegaTn, OtherState%OmegaDotTn, m%RtHS%YawFriMom ) !Compute yaw Friction #RRD
CALL YawFriction( t, p, m%FrcONcRt, m%MomONcRt, YawFriMz, OtherState%OmegaTn, OtherState%OmegaDotTn, m%RtHS%YawFriMom ) !Compute yaw Friction #RRD


!bjj: note m%RtHS%GBoxEffFac needed in OtherState only to fix HSSBrTrq (and used in FillAugMat)
Expand Down Expand Up @@ -3334,10 +3333,17 @@ SUBROUTINE SetPrimaryParameters( InitInp, p, InputFileData, ErrStat, ErrMsg )
p%TeetHSSp = 0.0
END IF

! Yaw friction model inputs
p%YawFrctMod = InputFileData%YawFrctMod
p%M_CD = InputFileData%M_CD
p%M_CSmax = InputFileData%M_CSmax
p%sig_v = InputFileData%sig_v
p%M_CD = InputFileData%M_CD
p%M_FCD = InputFileData%M_FCD
p%M_MCD = InputFileData%M_MCD
p%M_CSmax = InputFileData%M_CSmax
p%M_FCSmax = InputFileData%M_FCSmax
p%M_MCSmax = InputFileData%M_MCSmax
p%sig_v = InputFileData%sig_v
p%sig_v2 = InputFileData%sig_v2
p%OmgCut = InputFileData%OmgCut


CALL AllocAry( p%TipMass, p%NumBl, 'TipMass', ErrStat, ErrMsg )
Expand Down Expand Up @@ -6403,58 +6409,64 @@ SUBROUTINE Teeter( t, p, TeetDef, TeetRate, TeetMom )
END SUBROUTINE Teeter
!----------------------------------------------------------------------------------------------------------------------------------
!> This routine computes the Yaw Friction Torque due to yaw rate and acceleration.
SUBROUTINE YawFriction( t, p, Fz, Mzz, Omg, OmgDot, YawFriMf )
SUBROUTINE YawFriction( t, p, F, M, Mzz, Omg, OmgDot, YawFriMf )
!..................................................................................................................................

! Passed Variables:
REAL(DbKi), INTENT(IN) :: t !< simulation time
TYPE(ED_ParameterType), INTENT(IN) :: p !< parameters from the structural dynamics module
REAL(R8Ki), INTENT(IN ) :: Fz, Mzz !< Effective yaw bearing force and external yaw bearing torque
REAL(R8Ki), INTENT(IN ) :: Omg !< The yaw rate (rotational speed), x%QDT(DOF_Yaw).
REAL(R8Ki), INTENT(IN ) :: OmgDot !< The yaw acceleration (derivative of rotational speed), x%QD2T(DOF_Yaw).

REAL(ReKi), INTENT(OUT) :: YawFriMf !< The total friction torque (Coulomb + viscous).
REAL(ReKi), INTENT(IN ) :: F(3), M(3) !< Effective yaw bearing force and moment
REAL(R8Ki), INTENT(IN ) :: Mzz !< External yaw bearing torque
REAL(R8Ki), INTENT(IN ) :: Omg !< The yaw rate (rotational speed), x%QDT(DOF_Yaw).
REAL(R8Ki), INTENT(IN ) :: OmgDot !< The yaw acceleration (derivative of rotational speed), x%QD2T(DOF_Yaw).
REAL(ReKi), INTENT(OUT) :: YawFriMf !< The total friction torque (Coulomb + viscous).

! Local variables:
REAL(ReKi) :: temp ! It takes teh value of Fz or -1.
REAL(ReKi) :: temp, Fs, Mb, Mf_vis ! temp takes the value of Fz or -1.


SELECT CASE ( p%YawFrctMod )
! Yaw-friction model {0: none, 1: does not use Fz at yaw bearing, 2: does, 3: user defined model} (switch)
! Yaw-friction model {0: none, 1: does not use F and M at yaw bearing, 2: does, 3: user defined model} (switch)

CASE ( 0_IntKi ) ! None!


YawFriMf = 0.0_ReKi

CASE ( 1_IntKi, 2_IntKi ) ! 1 = F and M not used. 2 = F and M used

CASE ( 1_IntKi, 2_IntKi ) ! 1= no Fz use. 2=Fz used
temp = -1.0_ReKi ! In the case of YawFrctMod=1
Fs = 0.0_ReKi
Mb = 0.0_ReKi

temp = -1.0_ReKi !In the case of YawFrctMod=1

IF (p%YawFrctMod .EQ. 2) THEN
temp = MIN(0.0_R8Ki, Fz) !In the case of YawFrctMod=2
temp = MIN(0.0_ReKi, F(3)) ! In the case of YawFrctMod=2
Fs = SQRT(F(1)**2+F(2)**2) ! Effective shear force on yaw bearing
Mb = SQRT(M(1)**2+M(2)**2) ! Effective bending moment on yaw bearing
ENDIF

IF (EqualRealNos( Omg, 0.0_R8Ki ) )THEN
YawFriMf = -MIN(real(p%M_CD,ReKi) * ABS(temp), ABS(real(Mzz,ReKi))) * SIGN(1.0_ReKi, real(Mzz,ReKi))
IF (EqualRealNos( Omg, 0.0_R8Ki )) THEN
IF (EqualRealNos( OmgDot, 0.0_R8Ki )) THEN
YawFriMf = -MIN(real(p%M_CSmax,ReKi) * ABS(temp), ABS(real(Mzz,ReKi))) * SIGN(1.0_ReKi, real(Mzz,ReKi))
YawFriMf = -MIN( real(p%M_CSmax,ReKi) * ABS(temp) + real(p%M_FCSmax,ReKi) * Fs + real(p%M_MCSmax,ReKi) * Mb, ABS(real(Mzz,ReKi)) ) * SIGN(1.0_ReKi, real(Mzz,ReKi))
ELSE
YawFriMf = -MIN( real(p%M_CD, ReKi) * ABS(temp) + real(p%M_FCD, ReKi) * Fs + real(p%M_MCD, ReKi) * Mb, ABS(real(Mzz,ReKi)) ) * SIGN(1.0_ReKi, real(Mzz,ReKi))
ENDIF
ELSE
YawFriMf = real(p%M_CD,ReKi) * temp * sign(1.0_ReKi, real(Omg,ReKi)) - real(p%sig_v,ReKi) * real(Omg,ReKi)
! Viscous friction
IF ( ABS(Omg) > p%OmgCut ) THEN ! Full quadratic viscous friction
Mf_vis = - real(p%sig_v,ReKi) * real(Omg,ReKi) - real(p%sig_v2,ReKi) * real(Omg,ReKi) * ABS(real(Omg,ReKi))
ELSE ! Linearized viscous friction
Mf_vis = - ( real(p%sig_v,ReKi) + real(p%sig_v2,ReKi) * real(p%OmgCut,ReKi) ) * real(Omg,ReKi)
ENDIF
YawFriMf = ( real(p%M_CD,ReKi) * temp - real(p%M_FCD,ReKi) * Fs - real(p%M_MCD,ReKi) * Mb ) * sign(1.0_ReKi, real(Omg,ReKi)) & ! Coulomb friction
+ Mf_vis
ENDIF

CASE ( 3_IntKi ) ! User-defined YawFriMf model. >>>> NOT IMPLEMENTED YET

CASE ( 3_IntKi ) ! User-defined YawFriMf model. >>>> NOT IMPLEMENTED YET


CALL UserYawFrict ( t, Fz, Mzz, Omg, OmgDot, p%RootName, YawFriMf )

CALL UserYawFrict ( t, F, M, Mzz, Omg, OmgDot, p%RootName, YawFriMf )

END SELECT


RETURN
END SUBROUTINE YawFriction
!----------------------------------------------------------------------------------------------------------------------------------
Expand Down
Loading

0 comments on commit 4377c13

Please sign in to comment.