OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
i25therm.F
Go to the documentation of this file.
1Copyright> OpenRadioss
2Copyright> Copyright (C) 1986-2025 Altair Engineering Inc.
3Copyright>
4Copyright> This program is free software: you can redistribute it and/or modify
5Copyright> it under the terms of the GNU Affero General Public License as published by
6Copyright> the Free Software Foundation, either version 3 of the License, or
7Copyright> (at your option) any later version.
8Copyright>
9Copyright> This program is distributed in the hope that it will be useful,
10Copyright> but WITHOUT ANY WARRANTY; without even the implied warranty of
11Copyright> MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12Copyright> GNU Affero General Public License for more details.
13Copyright>
14Copyright> You should have received a copy of the GNU Affero General Public License
15Copyright> along with this program. If not, see <https://www.gnu.org/licenses/>.
16Copyright>
17Copyright>
18Copyright> Commercial Alternative: Altair Radioss Software
19Copyright>
20Copyright> As an alternative to this open-source version, Altair also offers Altair Radioss
21Copyright> software under a commercial license. Contact Altair to discuss further if the
22Copyright> commercial version may interest you: https://www.altair.com/radioss/.
23!||====================================================================
24!|| i25therm ../engine/source/interfaces/int25/i25therm.F
25!||--- called by ------------------------------------------------------
26!|| i25mainf ../engine/source/interfaces/int25/i25mainf.F
27!||--- calls -----------------------------------------------------
28!|| finter ../engine/source/tools/curve/finter.F
29!||====================================================================
30 SUBROUTINE i25therm(JLT ,KTHE ,TEMPI ,AREAS ,IELESI ,
31 2 IELEMI ,GAPV ,IFUNCTK,XTHE ,FNI ,
32 3 NPC ,TF ,FRAD ,DRAD ,EFRICT ,
33 4 FHEATS ,FHEATM,CONDINT,IFORM ,TEMP ,
34 5 H1 ,H2 ,H3 ,H4 ,FCOND ,
35 6 DCOND ,TINT ,XI ,YI ,ZI ,
36 7 X1 ,Y1 ,Z1 ,X2 ,Y2 ,
37 8 Z2 ,X3 ,Y3 ,Z3 ,X4 ,
38 9 Y4 ,Z4 ,IX1 ,IX2 ,IX3 ,
39 A IX4 ,PHI ,PHI1 ,PHI2 ,PHI3 ,
40 B PHI4 ,PM ,NSV ,ITAB ,THEACCFACT)
41C-----------------------------------------------
42C I m p l i c i t T y p e s
43C-----------------------------------------------
44#include "implicit_f.inc"
45C-----------------------------------------------
46C G l o b a l P a r a m e t e r s
47C-----------------------------------------------
48#include "mvsiz_p.inc"
49C-----------------------------------------------
50C C o m m o n B l o c k s
51C-----------------------------------------------
52#include "param_c.inc"
53#include "com08_c.inc"
54C-----------------------------------------------
55C D u m m y A r g u m e n t s
56C-----------------------------------------------
57 INTEGER JLT, IFUNCTK, FCOND,ITAB(*) ,NSV(*),
58 . NPC(*),IFORM,IELESI(MVSIZ) ,IELEMI(MVSIZ),
59 . IX1(MVSIZ), IX2(MVSIZ), IX3(MVSIZ), IX4(MVSIZ)
60 my_real, intent(in) :: theaccfact
61 my_real DCOND, TINT,
62 . PM(NPROPM,*)
63 my_real
64 . KTHE, XTHE, FRAD, DRAD, FHEATS, FHEATM,
65 . TEMPI(MVSIZ), XI(MVSIZ), YI(MVSIZ),TEMP(*),
66 . ZI(MVSIZ), PHI(MVSIZ), AREAS(MVSIZ), ASI(MVSIZ),
67 . bsi(mvsiz), gapv(mvsiz), condint(mvsiz),
68 . fni(mvsiz), tf(*), efrict(mvsiz),
69 . h1(mvsiz), h2(mvsiz), h3(mvsiz), h4(mvsiz),
70 . phi1(mvsiz),phi2(mvsiz),phi3(mvsiz),phi4(mvsiz),
71 . x1(mvsiz), x2(mvsiz), x3(mvsiz), x4(mvsiz),
72 . y1(mvsiz), y2(mvsiz), y3(mvsiz), y4(mvsiz),
73 . z1(mvsiz), z2(mvsiz), z3(mvsiz), z4(mvsiz)
74 my_real
75 . finter
76 EXTERNAL finter
77C-----------------------------------------------
78C L o c a l V a r i a b l e s
79C-----------------------------------------------
80 INTEGER I, MATS , MATM
81C REAL
82 my_real
83 . SX1 , SY1 , SZ1 , SX2 , SY2 , SZ2,
84 . ts, rstif, tstifm, tstift, cond, p, dydx,
85 . tm, ddcond, dd, hcond, dist, penrad, norm, areac,
86 . phim,area,conds,condm
87C-----------------------------------------------
88 phi(1:jlt) = zero
89 IF(ifunctk==0)THEN
90 rstif = one/max(em30,kthe)
91 DO i=1,jlt
92 ts = tempi(i)
93 IF(iform ==0) THEN
94 tm = tint
95 ELSE
96 tm = h1(i)*temp(ix1(i)) + h2(i)*temp(ix2(i))
97 . + h3(i)*temp(ix3(i)) + h4(i)*temp(ix4(i))
98 ENDIF
99 condint(i) = zero
100 ddcond = max(dcond-gapv(i),em20)
101C
102C------------------------------------------
103C CALCUL DE LA SURFACE VECTORIELLE (*2.)
104C------------------------------------------
105 sx1=(y1(i)-y3(i))*(z2(i)-z4(i)) - (z1(i)-z3(i))*(y2(i)-y4(i))
106 sy1=(z1(i)-z3(i))*(x2(i)-x4(i)) - (x1(i)-x3(i))*(z2(i)-z4(i))
107 sz1=(x1(i)-x3(i))*(y2(i)-y4(i)) - (y1(i)-y3(i))*(x2(i)-x4(i))
108C
109 norm = sqrt(sx1**2 + sy1**2 + sz1**2)
110C--------+---------+---------+---------+---------+---------+---------+--
111C CALCUL DE LA DISTANCE ENTRE CENTRE ET SURFACE
112C-------------------------------------------------------------
113 IF(ix3(i)/=ix4(i))THEN
114 sx2 = fourth*(x1(i) + x2(i) + x3(i) + x4(i)) - xi(i)
115 sy2 = fourth*(y1(i) + y2(i) + y3(i) + y4(i)) - yi(i)
116 sz2 = fourth*(z1(i) + z2(i) + z3(i) + z4(i)) - zi(i)
117 ELSE
118 sx2 = third*(x1(i) + x2(i) + x3(i)) - xi(i)
119 sy2 = third*(y1(i) + y2(i) + y3(i)) - yi(i)
120 sz2 = third*(z1(i) + z2(i) + z3(i)) - zi(i)
121 END IF
122C
123C-----------------------------------------------
124C Distance between secnd node
125C and main segment
126C-----------------------------------------------
127 dist = -(sx2*sx1+sy2*sy1+sz2*sz1) / max(em15,norm)
128
129 penrad = dist - gapv(i)
130C
131 IF(areas(i) == zero )THEN
132C
133 areac =half*norm
134C
135 ELSE
136C
137 areac = areas(i)
138C
139 ENDIF
140 tstifm = zero
141 cond = zero
142C---------------------------------
143C DISTANCE BETWEEN SECOND NODE
144C AND SURFACE MAIN
145C---------------------------------
146C-------------------------------------
147C Conduction : close distance
148C-------------------------------------
149
150 IF(penrad <= zero)THEN
151C---------------------------------
152C CONDUCTANCE COMPUTATION
153C---------------------------------
154
155 mats = ielesi(i)
156 IF(mats > 0) THEN
157 conds=pm(75,mats)+pm(76,mats)*ts
158 ELSE
159 conds =zero
160 ENDIF
161 IF(iform == 1 ) THEN
162 matm = ielemi(i)
163 IF(matm > 0) THEN
164 condm=pm(75,matm)+pm(76,matm)*tm
165 ELSE
166 condm = zero
167 ENDIF
168 IF(condm == zero) THEN
169 cond = conds
170 ELSEIF(conds == zero) THEN
171 cond = condm
172 ELSE
173 cond = two*conds*condm/(condm + conds)! harmonic mean
174 ENDIF
175 ELSE
176 cond = conds
177 ENDIF
178 IF(cond /= zero) tstifm = abs(dist) / cond
179 tstift = tstifm + rstif
180 condint(i) = areac * theaccfact /tstift
181C---------------------------------
182 phi(i) = areac*(tm - ts)*dt1*theaccfact / tstift
183C----------------------------------------------------------------------
184C Conduction + Radiation : Heat exchange depending on distance
185C----------------------------------------------------------------------
186 ELSEIF(penrad <= ddcond)THEN
187C---------------------------------
188 dist = gapv(i)
189
190 mats = ielesi(i)
191 IF(mats > 0) THEN
192 conds=pm(75,mats)+pm(76,mats)*ts
193 ELSE
194 conds =zero
195 ENDIF
196
197 IF(iform == 1 ) THEN
198 matm = ielemi(i)
199 IF(matm > 0) THEN
200 condm=pm(75,matm)+pm(76,matm)*tm
201 ELSE
202 condm = zero
203 ENDIF
204 IF(condm == zero) THEN
205 cond = conds
206 ELSEIF(conds == zero) THEN
207 cond = condm
208 ELSE
209 cond = two*conds*condm/(condm + conds)! harmonic mean
210 ENDIF
211 ELSE
212 cond = conds
213 ENDIF
214
215 IF(cond /= zero) tstifm = max(dist,zero) / cond
216 tstift = tstifm + rstif
217 dd = penrad /ddcond
218 hcond = finter(fcond,dd,npc,tf,dydx) / tstift
219 condint(i) = areac*hcond*theaccfact
220
221 phi(i) = areac * (tm - ts)*dt1* hcond*theaccfact
222
223 phi(i) = phi(i) + frad * areac * (tm*tm+ts*ts)
224 . * (tm + ts) * (tm - ts) * dt1 * theaccfact
225C-------------------------------------
226C Radiation :
227C-------------------------------------
228 ELSEIF(penrad <= drad)THEN
229C---------------------------------
230 phi(i) = frad * areac * (tm*tm+ts*ts)
231 . * (tm + ts) * (tm - ts) * dt1 * theaccfact
232 END IF
233C
234 IF(iform == 1 )THEN
235C
236 phi1(i) = -phi(i) *h1(i)
237 phi2(i) = -phi(i) *h2(i)
238 phi3(i) = -phi(i) *h3(i)
239 phi4(i) = -phi(i) *h4(i)
240C---------------------------------
241C HEAT GENERATION DUE TO FRICTION
242C---------------------------------
243 phim = fheatm * efrict(i)
244 phi1(i) = phi1(i) + phim*h1(i) ! main HEATING
245 phi2(i) = phi2(i) + phim*h2(i)
246 phi3(i) = phi3(i) + phim*h3(i)
247 phi4(i) = phi4(i) + phim*h4(i)
248
249 ENDIF
250C---------------------------------
251C HEAT GENERATION DUE TO FRICTION
252C---------------------------------
253 phi(i) = phi(i) + fheats * efrict(i) * theaccfact
254 ENDDO
255C
256 ELSE
257 DO i=1,jlt
258 ts = tempi(i)
259 IF(iform ==0) THEN
260 tm = tint
261 ELSE
262 tm = h1(i)*temp(ix1(i)) + h2(i)*temp(ix2(i))
263 . + h3(i)*temp(ix3(i)) + h4(i)*temp(ix4(i))
264 ENDIF
265 condint(i) = zero
266 ddcond = max(dcond-gapv(i),em20)
267C
268C------------------------------------------
269C CALCUL DE LA SURFACE VECTORIELLE (*2.)
270C------------------------------------------
271 sx1=(y1(i)-y3(i))*(z2(i)-z4(i)) - (z1(i)-z3(i))*(y2(i)-y4(i))
272 sy1=(z1(i)-z3(i))*(x2(i)-x4(i)) - (x1(i)-x3(i))*(z2(i)-z4(i))
273 sz1=(x1(i)-x3(i))*(y2(i)-y4(i)) - (y1(i)-y3(i))*(x2(i)-x4(i))
274C
275 norm = sqrt(sx1**2 + sy1**2 + sz1**2)
276C--------+---------+---------+---------+---------+---------+---------+--
277C CALCUL DE LA DISTANCE ENTRE CENTRE ET SURFACE
278C-------------------------------------------------------------
279 IF(ix3(i)/=ix4(i))THEN
280 sx2 = fourth*(x1(i) + x2(i) + x3(i) + x4(i)) - xi(i)
281 sy2 = fourth*(y1(i) + y2(i) + y3(i) + y4(i)) - yi(i)
282 sz2 = fourth*(z1(i) + z2(i) + z3(i) + z4(i)) - zi(i)
283 ELSE
284 sx2 = third*(x1(i) + x2(i) + x3(i)) - xi(i)
285 sy2 = third*(y1(i) + y2(i) + y3(i)) - yi(i)
286 sz2 = third*(z1(i) + z2(i) + z3(i)) - zi(i)
287 END IF
288C
289C-----------------------------------------------
290C Distance between secnd node
291C and main segment
292C-----------------------------------------------
293 dist = -(sx2*sx1+sy2*sy1+sz2*sz1) / max(em15,norm)
294
295 penrad = dist - gapv(i)
296C
297 IF(areas(i) == zero )THEN
298C
299 areac =half*norm
300C
301 ELSE
302C
303 areac = areas(i)
304C
305 ENDIF
306 tstifm = zero
307 cond = zero
308C-------------------------------------
309C Conduction : close distance
310C-------------------------------------
311 IF(penrad <= zero)THEN
312C---------------------------------
313C CALCUL DE LA CONDUCTIBILITE
314C---------------------------------
315 p = xthe * abs(fni(i)) / areas(i)
316 rstif = one / max(em30,kthe * finter(ifunctk,p,npc,tf,dydx))
317 mats = ielesi(i)
318 IF(mats > 0) THEN
319 conds=pm(75,mats)+pm(76,mats)*ts
320 ELSE
321 conds =zero
322 ENDIF
323
324 IF(iform == 1 ) THEN
325 matm = ielemi(i)
326 IF(matm > 0) THEN
327 condm=pm(75,matm)+pm(76,matm)*tm
328 ELSE
329 condm = zero
330 ENDIF
331 IF(condm == zero) THEN
332 cond = conds
333 ELSEIF(conds == zero) THEN
334 cond = condm
335 ELSE
336 cond = two*conds*condm/(condm + conds)! harmonic mean
337 ENDIF
338 ELSE
339 cond = conds
340 ENDIF
341
342 IF(cond /= zero) tstifm = abs(dist) / cond
343 tstift = tstifm + rstif
344 condint(i) = areac*theaccfact /tstift
345C---------------------------------
346 phi(i) = areac * (tm- ts)*dt1 * theaccfact / tstift
347C---------------------------------
348C----------------------------------------------------------------------
349C Conduction + Radiation : Heat exchange depending on distance
350C----------------------------------------------------------------------
351 ELSEIF(penrad <= ddcond)THEN
352C---------------------------------
353 dist = gapv(i)
354 mats = ielesi(i)
355 IF(mats > 0) THEN
356 conds=pm(75,mats)+pm(76,mats)*ts
357 ELSE
358 conds =zero
359 ENDIF
360
361 IF(iform == 1 ) THEN
362 matm = ielemi(i)
363 IF(matm > 0) THEN
364 condm=pm(75,matm)+pm(76,matm)*tm
365 ELSE
366 condm = zero
367 ENDIF
368 IF(condm == zero) THEN
369 cond = conds
370 ELSEIF(conds == zero) THEN
371 cond = condm
372 ELSE
373 cond = two*conds*condm/(condm + conds)! harmonic mean
374 ENDIF
375 ELSE
376 cond = conds
377 ENDIF
378 p = zero
379 rstif = one / max(em30,kthe * finter(ifunctk,p,npc,tf,dydx))
380
381 IF(cond /= zero)tstifm = max(dist,zero) / cond
382 tstift = tstifm + rstif
383
384 dd = penrad /ddcond
385 hcond = finter(fcond,dd,npc,tf,dydx) / tstift
386 condint(i) = areac*hcond*theaccfact
387
388 phi(i) = areac * (tm - ts)*dt1 * hcond * theaccfact
389
390 phi(i) = phi(i) + frad * areac * (tm*tm+ts*ts)
391 . * (tm + ts) * (tm - ts) * dt1 * theaccfact
392C-------------------------------------
393C Radiation :
394C-------------------------------------
395C---------------------------------
396 ELSEIF(penrad <= drad)THEN
397C---------------------------------
398 phi(i) = frad * areac * (tm*tm+ts*ts)
399 . * (tm + ts) * (tm - ts) * dt1 * theaccfact
400 END IF
401C
402 IF(iform == 1 )THEN
403C
404 phi1(i) = -phi(i) *h1(i)
405 phi2(i) = -phi(i) *h2(i)
406 phi3(i) = -phi(i) *h3(i)
407 phi4(i) = -phi(i) *h4(i)
408C---------------------------------
409C HEAT GENERATION DUE TO FRICTION
410C---------------------------------
411 phim = fheatm * efrict(i)
412 phi1(i) = phi1(i) + phim*h1(i) ! main HEATING
413 phi2(i) = phi2(i) + phim*h2(i)
414 phi3(i) = phi3(i) + phim*h3(i)
415 phi4(i) = phi4(i) + phim*h4(i)
416 ENDIF
417C---------------------------------
418C HEAT GENERATION DUE TO FRICTION
419C---------------------------------
420 phi(i) = phi(i) + fheats * efrict(i) * theaccfact
421 ENDDO
422 END IF
423C
424 RETURN
425 END
norm
norm(diag(diag(diag(inv(mat))) -id.SOL), 2) % destroy mumps instance id.JOB
area
subroutine area(d1, x, x2, y, y2, eint, stif0)
Definition hm_read_prop32.F:567
i25therm
subroutine i25therm(jlt, kthe, tempi, areas, ielesi, ielemi, gapv, ifunctk, xthe, fni, npc, tf, frad, drad, efrict, fheats, fheatm, condint, iform, temp, h1, h2, h3, h4, fcond, dcond, tint, xi, yi, zi, x1, y1, z1, x2, y2, z2, x3, y3, z3, x4, y4, z4, ix1, ix2, ix3, ix4, phi, phi1, phi2, phi3, phi4, pm, nsv, itab, theaccfact)
Definition i25therm.F:41
max
#define max(a, b)
Definition macros.h:21