2
! CalculiX - A 3-dimensional finite element program
3
! Copyright (C) 1998-2014 Guido Dhondt
5
! This program is free software; you can redistribute it and/or
6
! modify it under the terms of the GNU General Public License as
7
! published by the Free Software Foundation(version 2);
10
! This program is distributed in the hope that it will be useful,
11
! but WITHOUT ANY WARRANTY; without even the implied warranty of
12
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13
! GNU General Public License for more details.
15
! You should have received a copy of the GNU General Public License
16
! along with this program; if not, write to the Free Software
17
! Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
19
subroutine dflux(flux,sol,kstep,kinc,time,noel,npt,coords,
20
& jltyp,temp,press,loadtype,area,vold,co,lakonl,konl,
21
& ipompc,nodempc,coefmpc,nmpc,ikmpc,ilmpc,iscale,mi)
23
! user subroutine dflux
28
! sol current temperature value
30
! kinc increment number
31
! time(1) current step time
32
! time(2) current total time
34
! npt integration point number
35
! coords(1..3) global coordinates of the integration point
36
! jltyp loading face kode:
44
! temp currently not used
45
! press currently not used
46
! loadtype load type label
47
! area for surface flux: area covered by the
49
! for body flux: volume covered by the
51
! vold(0..4,1..nk) solution field in all nodes
53
! 1: displacement in global x-direction
54
! 2: displacement in global y-direction
55
! 3: displacement in global z-direction
57
! co(3,1..nk) coordinates of all nodes
58
! 1: coordinate in global x-direction
59
! 2: coordinate in global y-direction
60
! 3: coordinate in global z-direction
61
! lakonl element label
62
! konl(1..20) nodes belonging to the element
63
! ipompc(1..nmpc)) ipompc(i) points to the first term of
64
! MPC i in field nodempc
65
! nodempc(1,*) node number of a MPC term
66
! nodempc(2,*) coordinate direction of a MPC term
67
! nodempc(3,*) if not 0: points towards the next term
68
! of the MPC in field nodempc
69
! if 0: MPC definition is finished
70
! coefmpc(*) coefficient of a MPC term
71
! nmpc number of MPC's
72
! ikmpc(1..nmpc) ordered global degrees of freedom of the MPC's
73
! the global degree of freedom is
74
! 8*(node-1)+direction of the dependent term of
75
! the MPC (direction = 0: temperature;
76
! 1-3: displacements; 4: static pressure;
78
! ilmpc(1..nmpc) ilmpc(i) is the MPC number corresponding
79
! to the reference number in ikmpc(i)
80
! mi(1) max # of integration points per element (max
82
! mi(2) max degree of freedomm per node (max over all
83
! nodes) in fields like v(0:mi(2))...
87
! flux(1) magnitude of the flux
88
! flux(2) not used; please do NOT assign any value
89
! iscale determines whether the flux has to be
90
! scaled for increments smaller than the
91
! step time in static calculations
93
! 1: scaling (default)
100
integer kstep,kinc,noel,npt,jltyp,konl(20),ipompc(*),
101
& nodempc(3,*),nmpc,ikmpc(*),ilmpc(*),node,idof,id,iscale,mi(*)
103
real*8 flux(2),time(2),coords(3),sol,temp,press,vold(0:mi(2),*),
104
& area,co(3,*),coefmpc(*)
106
! the code starting here up to the end of the file serves as
107
! an example for combined mechanical-lubrication problems.
108
! Please replace it by your own code for your concrete application.
110
integer ifaceq(9,6),ifacet(7,4),ifacew(8,5),ig,nelem,nopes,
113
real*8 xl21(3,9),xi,et,al,rho,um,h,pnode1(3),pnode2(3),
114
& ratio(9),dist,xl22(3,9),dpnode1(3,3),dpnode2(3,3),v1(3),
115
& v2(3),dh(3),xsj2(3),xs2(3,7),shp2(7,8)
117
data ifaceq /4,3,2,1,11,10,9,12,21,
118
& 5,6,7,8,13,14,15,16,22,
119
& 1,2,6,5,9,18,13,17,23,
120
& 2,3,7,6,10,19,14,18,24,
121
& 3,4,8,7,11,20,15,19,25,
122
& 4,1,5,8,12,17,16,20,26/
123
data ifacet /1,3,2,7,6,5,11,
127
data ifacew /1,3,2,9,8,7,0,0,
128
& 4,5,6,10,11,12,0,0,
129
& 1,2,5,4,7,14,10,13,
130
& 2,3,6,5,8,15,11,14,
131
& 4,6,3,1,12,15,9,13/
138
if(lakonl(4:4).eq.'2') then
141
elseif(lakonl(4:4).eq.'8') then
144
elseif(lakonl(4:5).eq.'10') then
147
elseif(lakonl(4:4).eq.'4') then
150
elseif(lakonl(4:5).eq.'15') then
152
elseif(lakonl(4:4).eq.'6') then
156
! treatment of wedge faces
158
if(lakonl(4:4).eq.'6') then
165
if(lakonl(4:5).eq.'15') then
173
! first side of the oil film
177
if((nope.eq.20).or.(nope.eq.8)) then
179
node=konl(ifaceq(i,ig))
181
call nident(ikmpc,idof,nmpc,id)
182
if((id.eq.0).or.(ikmpc(id).ne.idof)) then
183
write(*,*) '*ERROR in dflux: node ',node
184
write(*,*) ' is not connected to the structure'
187
node=nodempc(1,nodempc(3,ipompc(ilmpc(id))))
189
xl21(j,i)=co(j,node)+
193
elseif((nope.eq.10).or.(nope.eq.4)) then
194
write(*,*) '*ERROR in dload: tetrahedral elements'
195
write(*,*) ' are not allowed'
199
node=konl(ifacew(i,ig))
201
call nident(ikmpc,idof,nmpc,id)
202
if((id.eq.0).or.(ikmpc(id).ne.idof)) then
203
write(*,*) '*ERROR in dflux: node ',node
204
write(*,*) ' is not connected to the structure'
207
node=nodempc(1,nodempc(3,ipompc(ilmpc(id))))
209
xl21(j,i)=co(j,node)+
215
! projecting the integration point on the first side of the
222
call attach(xl21,pnode1,nopes,ratio,dist,xi,et)
224
! derivative of the shape functions in (xi,et)
227
call shape8q(xi,et,xl21,xsj2,xs2,shp2,iflag)
228
elseif(nopes.eq.4) then
229
call shape4q(xi,et,xl21,xsj2,xs2,shp2,iflag)
230
elseif(nopes.eq.6) then
231
call shape6tri(xi,et,xl21,xsj2,xs2,shp2,iflag)
233
call shape3tri(xi,et,xl21,xsj2,xs2,shp2,iflag)
236
! the gradient of pnode1
237
! dpnode1(j,k)=dpnode1(j)/dx(k)
243
dpnode1(i,j)=dpnode1(i,j)+shp2(j,k)*xl21(i,k)
248
! second side of the oil film
252
if((nope.eq.20).or.(nope.eq.8)) then
254
node=konl(ifaceq(i,ig))
256
call nident(ikmpc,idof,nmpc,id)
257
if((id.eq.0).or.(ikmpc(id).ne.idof)) then
258
write(*,*) '*ERROR in dflux: node ',node
259
write(*,*) ' is not connected to the structure'
262
node=nodempc(1,nodempc(3,ipompc(ilmpc(id))))
264
xl22(j,i)=co(j,node)+
268
elseif((nope.eq.10).or.(nope.eq.4)) then
269
write(*,*) '*ERROR in dload: tetrahedral elements'
270
write(*,*) ' are not allowed'
274
node=konl(ifacew(i,ig))
276
call nident(ikmpc,idof,nmpc,id)
277
if((id.eq.0).or.(ikmpc(id).ne.idof)) then
278
write(*,*) '*ERROR in dflux: node ',node
279
write(*,*) ' is not connected to the structure'
282
node=nodempc(1,nodempc(3,ipompc(ilmpc(id))))
284
xl22(j,i)=co(j,node)+
290
! projecting the integration point on the second side of the
297
call attach(xl22,pnode2,nopes,ratio,dist,xi,et)
299
! derivative of the shape functions in (xi,et)
302
call shape8q(xi,et,xl22,xsj2,xs2,shp2,iflag)
303
elseif(nopes.eq.4) then
304
call shape4q(xi,et,xl22,xsj2,xs2,shp2,iflag)
305
elseif(nopes.eq.6) then
306
call shape6tri(xi,et,xl22,xsj2,xs2,shp2,iflag)
308
call shape3tri(xi,et,xl22,xsj2,xs2,shp2,iflag)
311
! the gradient of pnode1
312
! dpnode2(j,k)=dpnode2(j)/dx(k)
318
dpnode2(i,j)=dpnode2(i,j)+shp2(j,k)*xl22(i,k)
323
! calculating the thickness of the oil film
325
h=dsqrt((pnode1(1)-pnode2(1))**2+
326
& (pnode1(2)-pnode2(2))**2+
327
& (pnode1(3)-pnode2(3))**2)
329
! calculating the gradient of the oil film thickness
332
dh(i)=((pnode1(1)-pnode2(1))*(dpnode1(1,i)-dpnode2(1,i))
333
& +(pnode1(2)-pnode2(2))*(dpnode1(2,i)-dpnode2(2,i))
334
& +(pnode1(3)-pnode2(3))*(dpnode1(3,i)-dpnode2(3,i)))/h
337
! velocity of the parts adjoining the film
338
! the axis or rotation is assumed to be the x-axis
344
v2(2)=-26000.d0*coords(3)/dsqrt(coords(2)**2+coords(3)**2)
345
v2(3)=26000.d0*coords(2)/dsqrt(coords(2)**2+coords(3)**2)
347
! density (oil, N-mm-s-K system)
353
flux(1)=-rho*((v1(1)+v2(1))*dh(1)+
354
& (v1(2)+v2(2))*dh(2)+
355
& (v1(3)+v2(3))*dh(3))/2.d0
added
removed
CalculiX/ccx_2.7/src/dflux.f
