2
! CalculiX - A 3-dimensional finite element program
3
! Copyright (C) 1998-2015 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(8,6),ifacet(6,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
intent(in) sol,kstep,kinc,time,noel,npt,coords,
118
& jltyp,temp,press,loadtype,area,vold,co,lakonl,konl,
119
& ipompc,nodempc,coefmpc,nmpc,ikmpc,ilmpc,mi
121
intent(out) flux,iscale
125
ifaceq=reshape((/4,3,2,1,11,10,9,12,
126
& 5,6,7,8,13,14,15,16,
127
& 1,2,6,5,9,18,13,17,
128
& 2,3,7,6,10,19,14,18,
129
& 3,4,8,7,11,20,15,19,
130
& 4,1,5,8,12,17,16,20/),(/8,6/))
131
ifacet=reshape((/1,3,2,7,6,5,
134
& 1,4,3,8,10,7/),(/6,4/))
135
ifacew=reshape((/1,3,2,9,8,7,0,0,
136
& 4,5,6,10,11,12,0,0,
137
& 1,2,5,4,7,14,10,13,
138
& 2,3,6,5,8,15,11,14,
139
& 4,6,3,1,12,15,9,13/),(/8,5/))
144
if(lakonl(4:4).eq.'2') then
147
elseif(lakonl(4:4).eq.'8') then
150
elseif(lakonl(4:5).eq.'10') then
153
elseif(lakonl(4:4).eq.'4') then
156
elseif(lakonl(4:5).eq.'15') then
158
elseif(lakonl(4:4).eq.'6') then
162
! treatment of wedge faces
164
if(lakonl(4:4).eq.'6') then
171
if(lakonl(4:5).eq.'15') then
179
! first side of the oil film
183
if((nope.eq.20).or.(nope.eq.8)) then
185
node=konl(ifaceq(i,ig))
187
call nident(ikmpc,idof,nmpc,id)
188
if((id.eq.0).or.(ikmpc(id).ne.idof)) then
189
write(*,*) '*ERROR in dflux: node ',node
190
write(*,*) ' is not connected to the structure'
193
node=nodempc(1,nodempc(3,ipompc(ilmpc(id))))
195
xl21(j,i)=co(j,node)+
199
elseif((nope.eq.10).or.(nope.eq.4)) then
200
write(*,*) '*ERROR in dload: tetrahedral elements'
201
write(*,*) ' are not allowed'
205
node=konl(ifacew(i,ig))
207
call nident(ikmpc,idof,nmpc,id)
208
if((id.eq.0).or.(ikmpc(id).ne.idof)) then
209
write(*,*) '*ERROR in dflux: node ',node
210
write(*,*) ' is not connected to the structure'
213
node=nodempc(1,nodempc(3,ipompc(ilmpc(id))))
215
xl21(j,i)=co(j,node)+
221
! projecting the integration point on the first side of the
228
call attach(xl21,pnode1,nopes,ratio,dist,xi,et)
230
! derivative of the shape functions in (xi,et)
233
call shape8q(xi,et,xl21,xsj2,xs2,shp2,iflag)
234
elseif(nopes.eq.4) then
235
call shape4q(xi,et,xl21,xsj2,xs2,shp2,iflag)
236
elseif(nopes.eq.6) then
237
call shape6tri(xi,et,xl21,xsj2,xs2,shp2,iflag)
239
call shape3tri(xi,et,xl21,xsj2,xs2,shp2,iflag)
242
! the gradient of pnode1
243
! dpnode1(j,k)=dpnode1(j)/dx(k)
249
dpnode1(i,j)=dpnode1(i,j)+shp2(j,k)*xl21(i,k)
254
! second side of the oil film
258
if((nope.eq.20).or.(nope.eq.8)) then
260
node=konl(ifaceq(i,ig))
262
call nident(ikmpc,idof,nmpc,id)
263
if((id.eq.0).or.(ikmpc(id).ne.idof)) then
264
write(*,*) '*ERROR in dflux: node ',node
265
write(*,*) ' is not connected to the structure'
268
node=nodempc(1,nodempc(3,ipompc(ilmpc(id))))
270
xl22(j,i)=co(j,node)+
274
elseif((nope.eq.10).or.(nope.eq.4)) then
275
write(*,*) '*ERROR in dload: tetrahedral elements'
276
write(*,*) ' are not allowed'
280
node=konl(ifacew(i,ig))
282
call nident(ikmpc,idof,nmpc,id)
283
if((id.eq.0).or.(ikmpc(id).ne.idof)) then
284
write(*,*) '*ERROR in dflux: node ',node
285
write(*,*) ' is not connected to the structure'
288
node=nodempc(1,nodempc(3,ipompc(ilmpc(id))))
290
xl22(j,i)=co(j,node)+
296
! projecting the integration point on the second side of the
303
call attach(xl22,pnode2,nopes,ratio,dist,xi,et)
305
! derivative of the shape functions in (xi,et)
308
call shape8q(xi,et,xl22,xsj2,xs2,shp2,iflag)
309
elseif(nopes.eq.4) then
310
call shape4q(xi,et,xl22,xsj2,xs2,shp2,iflag)
311
elseif(nopes.eq.6) then
312
call shape6tri(xi,et,xl22,xsj2,xs2,shp2,iflag)
314
call shape3tri(xi,et,xl22,xsj2,xs2,shp2,iflag)
317
! the gradient of pnode1
318
! dpnode2(j,k)=dpnode2(j)/dx(k)
324
dpnode2(i,j)=dpnode2(i,j)+shp2(j,k)*xl22(i,k)
329
! calculating the thickness of the oil film
331
h=dsqrt((pnode1(1)-pnode2(1))**2+
332
& (pnode1(2)-pnode2(2))**2+
333
& (pnode1(3)-pnode2(3))**2)
335
! calculating the gradient of the oil film thickness
338
dh(i)=((pnode1(1)-pnode2(1))*(dpnode1(1,i)-dpnode2(1,i))
339
& +(pnode1(2)-pnode2(2))*(dpnode1(2,i)-dpnode2(2,i))
340
& +(pnode1(3)-pnode2(3))*(dpnode1(3,i)-dpnode2(3,i)))/h
343
! velocity of the parts adjoining the film
344
! the axis or rotation is assumed to be the x-axis
350
v2(2)=-26000.d0*coords(3)/dsqrt(coords(2)**2+coords(3)**2)
351
v2(3)=26000.d0*coords(2)/dsqrt(coords(2)**2+coords(3)**2)
353
! density (oil, N-mm-s-K system)
359
flux(1)=-rho*((v1(1)+v2(1))*dh(1)+
360
& (v1(2)+v2(2))*dh(2)+
361
& (v1(3)+v2(3))*dh(3))/2.d0
added
removed
CalculiX/ccx_2.11/src/dflux.f
