-
Notifications
You must be signed in to change notification settings - Fork 1
/
Asp_calc.f90
261 lines (237 loc) · 12.5 KB
/
Asp_calc.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
! Subroutine for calculating the shape of the surface roughness.
!The shape is either a anayltical equatio och measurements of a surface
SUBROUTINE Asp_calc(H00, t, SS, NYs, NN)
implicit none
! Input
real H00
integer t, SS, NYs, NN
! Calculation Parameters
integer I, J, JJ
real H23, PAI, asph1
real position, Ab, Wb, bump, bump2, bumpL, H0, asp
! Other
real aa, sign
DATA PAI/3.14159265/
include 'inc_Grid.h'
include 'inc_Itt.h'
include 'inc_Asp.h'
include 'inc_Ref.h'
include 'inc_CurrentH.h'
include 'inc_Current.h'
include 'inc_COMH.h'
include 'inc_Setup.h'
include 'inc_Init_H00.h'
! Output
save /currentH/
save /Init_H00/
! Calculate the dimensionless position of the surface defect.
position=X0+(XE-X0)*t/Ntime ! = Ua/Um*t/Ntime*Tend !D_IN: /Grid/ -> XE, X0; /Itt/ -> Ntime
! This part is used for simulatios where the load increases over the time. This part is not fully developed yet
H23=-2.0 ! This is a constant for how much the film thickness hsould change. Idealy whould be to have it as in input parameter
If( asp_shape .EQ. 122)then !D_IN: /Ref/ -> asp_shape
if( t .EQ. 0) H00_t0=H00 !D_OUT: H00_t0 -> /Init_H00/
if( t .GT. 0) then
H00=H00_t0+(H23-H00_t0)*t/Ntime
endif
endif
! Choice of asperity chape ------------------------------------------------------------------------------------------------
! The parameter asp_shape controuls which equations should be used
! No asperity
IF(asp_shape .EQ. 0 .OR. t .LT. -10)THEN
DO J=1,NN,SS
DO I=1,NX,SS !D_IN: /Grid/ -> NX
H0=RAD(I,J)+W(I,J) !D_IN: /COMH/ -> RAD; /Current/ -> W
H(I,J)=H0 +H00 !D_OUT: H(I,J) -> /CurrentH/
ENDDO
ENDDO
! Asperity shape based upone Venner CH and Lubrecht AA. Transient analysis of surfacefeatures in an EHL line contact in the case of sliding. Trans ASME J Tribol 1994; 116: 186–193.
ELSE IF( asp_shape .EQ. 1) THEN
DO J=1,NN,SS
DO I=1,NX,SS
asp=asph*(10**(-10*((X(I)-position)/aspw)**2)) !D_IN: /Setup/ -> X(I); /asp/ -> asph, aspw
bump=asp*cos(2*PAI*((X(I)-position)/aspw))
H0=RAD(I,J)+W(I,J)+bump
H(I,J)=H0+H00
ENDDO
ENDDO
! Input parameters acc Holmes et al. Transient EHL point contact analysis 2003. Also the same as Venner 1994 ball and cylinder
ELSE IF(asp_shape .EQ. 3)THEN
DO J=1,NN,SS
DO I=1,NX,SS
Ab=asph
Wb=aspw
bump=Ab*10**(-10*((X(I)-position)/Wb)**2)*cos(2*PAI*(X(I)-position)/Wb)
H0=RAD(I,J)+W(I,J)-bump
H(I,J)=H0 +H00
ENDDO
ENDDO
! The asperity shape used by D. Hannes in D. Hannes 2011 Rolling contact fatigue crack path prediction but formulates as a line defect
ELSE IF(asp_shape .EQ. 6) THEN
DO J=1,NN,SS
DO I=1,NX,SS
bump=0
IF ( abs(X(I)-position)/aspw .LT. 1.0)THEN
bump=asph/2*cos((X(I)-position)*PAI/(aspw))+asph/2;
ENDIF
H0=RAD(I,J)+W(I,J)-bump
H(I,J)=H0 +H00
ENDDO
ENDDO
! The asperity shape used by D. Hannes in Rolling contact fatigue crack path prediction 2011 but formulates as a line defect plus on as a point defect
ELSE IF(asp_shape .EQ. 7) THEN
DO J=1,NN,SS
DO I=1,NX,SS
! line
bumpL=0
IF ( (X(I)-position)**2/aspw2 .LT. 1.0)THEN
bumpL=asph2/2*cos((X(I)-position)*PAI/(aspw2))+asph2/2; !D_IN: /asp/ -> aspw2, asph2
ENDIF
! point
bump=0
IF ( ((X(I)-position)**2+(Y(J))**2)**(0.5)/aspw .LT. 1.0)THEN
bump=asph/2*cos(((X(I)-position)**2+(asp_ratio*Y(J))**2)**(0.5)*2*PAI/(aspw))+asph/2;
ENDIF
H0=RAD(I,J)+W(I,J)-bump-bumpL
H(I,J)=H0 +H00
ENDDO
ENDDO
! Two times the asperity shape used by D. Hannes in D. Hannes 2011 Rolling contact fatigue crack path prediction
ELSE IF(asp_shape .EQ. 8) THEN
DO J=1,NN
DO I=1,NX
! Point
bump=0
IF ( 2*((X(I)-position)**2+(Y(J))**2)**(0.5)/aspw .LT. 1.0)THEN
bump=asph/2*cos(((X(I)-position)**2+(Y(J))**2)**(0.5)*2*PAI/(aspw))+asph/2;
ENDIF
! Point
bump2=0
IF ( 2*((X(I)-position)**2+(Y(J))**2)**(0.5)/aspw2 .LT. 1.0)THEN
bump2=asph2/2*cos(((X(I)-position)**2+(Y(J))**2)**(0.5)*2*PAI/(aspw2))+asph2/2;
ENDIF
H0=RAD(I,J)+W(I,J)-bump-bump2
H(I,J)=H0 +H00
ENDDO
ENDDO
! Input parameters acc Holmes et al. Transient EHL point contact analysis 2003. Also the sam as Venner 1994 ball and cylinder But reshaped as a point defect
ELSE IF(asp_shape .EQ. 13)THEN
DO J=1,NN,SS
DO I=1,NX,SS
Ab=asph
Wb=aspw
bump=Ab*10**(-10*((X(I)-position)**2+Y(J)**2)/Wb**2)*cos(2*PAI*sqrt((X(I)-position)**2+Y(J)**2)/Wb) !D_IN: /Setup/ -> Y(J)
H0=RAD(I,J)+W(I,J)-bump
H(I,J)=H0 +H00
ENDDO
ENDDO
! The asperity shape used by D. Hannes in Rolling contact fatigue crack path prediction 2011 with variable length to width ratio
ELSE IF(asp_shape .EQ. 12 .OR. asp_shape .EQ. 122)THEN
DO J=1,NN,SS
DO I=1,NX,SS
bump=0
IF ( 2*((X(I)-position)**2+(asp_ratio*Y(J))**2)**(0.5)/aspw .LT. 1.0)THEN !D_IN: /asp/ -> asp_ratio
bump=asph/2*cos(((X(I)-position)**2+(asp_ratio*Y(J))**2)**(0.5)*2*PAI/(aspw))+asph/2;
ENDIF
H0=RAD(I,J)+W(I,J)-bump
H(I,J)=H0 +H00
ENDDO
ENDDO
! The asperity shape used by D. Hannes in Rolling contact fatigue crack path prediction 2011 with variable length to width ratio
ELSE IF(asp_shape .EQ. 112)THEN
DO J=1,NN,SS
DO I=1,NX,SS
bump=0
IF ( 2*((X(I)-position)**2+(asp_ratio*Y(J))**2)**(0.5)/aspw .LT. 1.0)THEN ! We're at the asperity
asph1=asph/2*cos(((X(I)-position)**2+(asp_ratio*Y(J))**2)**(0.5)*2*PAI/(aspw))+asph/2
asph1=asph1*(2.0)**(0.5) ! To make highest poing = asph and lowest = -asph
aa=X(I)-position
IF ( 2*abs(aa)/aspw .LT. 0.5)THEN ! Central part
bump=asph1*sin(aa*2*PAI/(aspw));
ELSE ! At the outer region, half the amplitude and the wavelength to get continiuos derivatives.
IF (X(I) .LT. Position) then
bump=asph1/2*cos(2*aa*2*PAI/(aspw))-asph1/2 ! rise or lower the mean value depending on if in front or behind the center of the asperity
else
bump=-asph1/2*cos(2*aa*2*PAI/(aspw))+asph1/2
endif
endif
ENDIF
H0=RAD(I,J)+W(I,J)-bump
H(I,J)=H0 +H00
ENDDO
ENDDO
! The asperity shape used by D. Hannes in Rolling contact fatigue crack path prediction 2011 with changes length to width ratio for multiple defects
ELSE IF(asp_shape .EQ. 32)THEN
DO J=1,NN,SS
DO I=1,NX,SS
bump=0
IF ( 2*((X(I)-position)**2+(asp_ratio*Y(J))**2)**(0.5)/aspw .LT. 1.0)THEN
bump=asph/2*cos(((X(I)-position)**2+(asp_ratio*Y(J))**2)**(0.5)*2*PAI/(aspw))+asph/2
ELSE IF ( 2*((X(I)-position+aspw)**2+(asp_ratio*Y(J))**2)**(0.5)/aspw .LT. 1.0)THEN
bump=asph/2*cos(((X(I)-position+1*aspw)**2+(asp_ratio*Y(J))**2)**(0.5)*2*PAI/(aspw))+asph/2
ELSE IF ( 2*((X(I)-position+2*aspw)**2+(asp_ratio*Y(J))**2)**(0.5)/aspw .LT. 1.0)THEN
bump=asph/2*cos(((X(I)-position+2*aspw)**2+(asp_ratio*Y(J))**2)**(0.5)*2*PAI/(aspw))+asph/2
ELSE IF ( 2*((X(I)-position+3*aspw)**2+(asp_ratio*Y(J))**2)**(0.5)/aspw .LT. 1.0)THEN
bump=asph/2*cos(((X(I)-position+3*aspw)**2+(asp_ratio*Y(J))**2)**(0.5)*2*PAI/(aspw))+asph/2
ELSE IF ( 2*((X(I)-position+4*aspw)**2+(asp_ratio*Y(J))**2)**(0.5)/aspw .LT. 1.0)THEN
bump=asph/2*cos(((X(I)-position+4*aspw)**2+(asp_ratio*Y(J))**2)**(0.5)*2*PAI/(aspw))+asph/2
ENDIF
H0=RAD(I,J)+W(I,J)-bump
H(I,J)=H0 +H00
ENDDO
ENDDO
! The asperity shape used by D. Hannes in Rolling contact fatigue crack path prediction 2011 with changes length to width ratio for multiple defects but with no offset. I.E half asperity half dent.
ELSE IF(asp_shape .EQ. 33)THEN
DO J=1,NN,SS
DO I=1,NX,SS
bump=0
IF ( 2*((X(I)-position)**2+(asp_ratio*Y(J))**2)**(0.5)/aspw .LT. 1.0)THEN
bump=asph/2*cos(((X(I)-position)**2+(asp_ratio*Y(J))**2)**(0.5)*2*PAI/(aspw))
ELSE IF ( 2*((X(I)-position+aspw)**2+(asp_ratio*Y(J))**2)**(0.5)/aspw .LT. 1.0)THEN
bump=asph/2*cos(((X(I)-position+1*aspw)**2+(asp_ratio*Y(J))**2)**(0.5)*2*PAI/(aspw))
ELSE IF ( 2*((X(I)-position+2*aspw)**2+(asp_ratio*Y(J))**2)**(0.5)/aspw .LT. 1.0)THEN
bump=asph/2*cos(((X(I)-position+2*aspw)**2+(asp_ratio*Y(J))**2)**(0.5)*2*PAI/(aspw))
ELSE IF ( 2*((X(I)-position+3*aspw)**2+(asp_ratio*Y(J))**2)**(0.5)/aspw .LT. 1.0)THEN
bump=asph/2*cos(((X(I)-position+3*aspw)**2+(asp_ratio*Y(J))**2)**(0.5)*2*PAI/(aspw))
ELSE IF ( 2*((X(I)-position+4*aspw)**2+(asp_ratio*Y(J))**2)**(0.5)/aspw .LT. 1.0)THEN
bump=asph/2*cos(((X(I)-position+4*aspw)**2+(asp_ratio*Y(J))**2)**(0.5)*2*PAI/(aspw))
ENDIF
H0=RAD(I,J)+W(I,J)-bump
H(I,J)=H0 +H00
ENDDO
ENDDO
! Using real surface measurements imported by the Asp_read routine
ELSEIF( asp_shape .ge. 200 .and. asp_shape .le. 209) then
if( t .LT. 1) then
DO J=1,NN,SS
JJ=NYs-J+1
DO I=1,NX,SS
H0=RAD(I,J)+W(I,J)
H(I,J)=H0 +H00
H(I,JJ)=H(I,J)
ENDDO
ENDDO
else
Call Asp_surf_calc(NX,NYs,NN,SS, t, H00) !CALL: Only if asperity surface is imported Nx, Nys, NN = (Nys+1)/2, SS, t, H00
endif
! The asperity shape used by D. Hannes in Rolling contact fatigue crack path prediction 2011 formulated as a indent
ELSE
DO J=1,NN,SS
DO I=1,NX,SS
bump=0
IF ( 2*((X(I)-position)**2+(Y(J))**2)**(0.5)/aspw .LT. 1.0)THEN
bump=asph/2*cos(((X(I)-position)**2+(Y(J))**2)**(0.5)*2*PAI/(aspw))+asph/2;
ENDIF
H0=RAD(I,J)+W(I,J)-bump
H(I,J)=H0 +H00
ENDDO
ENDDO
ENDIF
!Mirror the surface around the symmetry line
DO J=1,NN,SS
JJ=NYs-J+1
DO I=1,NX,SS
H(I,JJ)=H(I,J)
ENDDO
ENDDO
RETURN
END