Newer
Older
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
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! ||===\\ |
! || \\ |
! || || //==\\ || || //==|| ||/==\\ |
! || || || || || || || || || || |
! || // || || || || || || || |
! ||===// \\==// \\==\\ \\==\\ || |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! MAKE_PRESSURE Apr. 2007 |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
subroutine make_pressure (params,icon,xg,yg,zg,viscosity0,penalty,expon, &
unode,vnode,wnode,temp,pressure,strain,nnode, &
r0,s0,t0,rst,plasticity_type,plasticity_parameters,eviscosity)
use definitions
use invariants
!------------------------------------------------------------------------------|
!(((((((((((((((( Purpose of the routine ))))))))))))))))))))))))))))))))))))))
!------------------------------------------------------------------------------|
! this routine is a stripped down version of make_matrix but is used in the
! calculation of the pressure
! icon is connectivity array for the current element
! xg,yg,zg are the global coordinate arrays of length nnode
! viscosity0 is viscosity
! penalty is incompressibility penalty factor
! expon is nonlinear viscosity exponent
! unode,vnode,wnode is the velocity array (obtained from previous time step
! or at least containing the proper velocity at the fixed dofs)
! temp and pressure anre temperature and pressure
! strain is total accumulated strain as interpolated from the 3D cloud
! nnode is number of nodes
! r0,s0,t0 are the bottom left back corner of the part of the element
! we are now integrating (in local r,s,t coordinates)
! rst is the size of the part of the element we are integrating
!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine arguments ))))))))))))))))))))
!------------------------------------------------------------------------------|
implicit none
type(parameters) params
integer icon(params%mpe)
double precision xg(nnode),yg(nnode),zg(nnode)
double precision viscosity0
double precision penalty
double precision expon
double precision unode(nnode),vnode(nnode),wnode(nnode)
double precision temp(nnode)
double precision pressure
double precision strain(nnode)
integer nnode
double precision r0,s0,t0,rst
character (len=8) plasticity_type
double precision plasticity_parameters(9)
double precision eviscosity
!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
!------------------------------------------------------------------------------|
!integer nint
integer iproc,nproc,ierr,ic,k,iint,err
double precision viscosity,straintot
double precision exx,eyy,ezz,exy,eyz,ezx,E2d,E3d
double precision f1(3,3),defgrad(3,3)
double precision detf1,detf,e2dref,viscomean
double precision velox,veloy,veloz,tau,uvwnorm,xmin,xmax,ymin,ymax
double precision zmin,zmax,dx,dy,dz,alpha,div
double precision eps,r,s,t,w,volume
double precision jcb(3,3),jcbi(3,3),jcbp(3,3),jcbip(3,3)
!double precision,dimension(:),allocatable :: rr,ss,tt,ww,h
double precision,dimension(:),allocatable :: x,y,z
double precision,dimension(:),allocatable :: dhdx,dhdy,dhdz
logical is_plastic_dummy,flag_vrm_pb
double precision :: rr,ss,tt,ww
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
INCLUDE 'mpif.h'
call mpi_comm_size (mpi_comm_world,nproc,ierr)
call mpi_comm_rank (mpi_comm_world,iproc,ierr)
eps=tiny(eps)
allocate (x(params%mpe),stat=err) ; if (err.ne.0) call stop_run ('Error alloc x in make_pressure$')
allocate (y(params%mpe),stat=err) ; if (err.ne.0) call stop_run ('Error alloc y in make_pressure$')
allocate (z(params%mpe),stat=err) ; if (err.ne.0) call stop_run ('Error alloc z in make_pressure$')
allocate (dhdx(params%mpe),stat=err) ; if (err.ne.0) call stop_run ('Error alloc dhdx in make_pressure$')
allocate (dhdy(params%mpe),stat=err) ; if (err.ne.0) call stop_run ('Error alloc dhdy in make_pressure$')
allocate (dhdz(params%mpe),stat=err) ; if (err.ne.0) call stop_run ('Error alloc dhdz in make_pressure$')
rr=0.d0
ss=0.d0
tt=0.d0
ww=8.d0
do k=1,params%mpe
x(k)=xg(icon(k))
y(k)=yg(icon(k))
Dave Whipp
committed
z(k)=zg(icon(k))*params%vex
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
enddo
! integrate over nint integration points
! here nint has been chosen to 1 because incompresibility
! is imposed in the FE matrix with a single integration point
r=r0+rst*(rr+1.d0)/2.d0
s=s0+rst*(ss+1.d0)/2.d0
t=t0+rst*(tt+1.d0)/2.d0
w=ww
! added by Jean 13/2/08
r=rr
s=ss
t=tt
w=ww
!end addition
call compute_dhdx_dhdy_dhdz (params%mpe,r,s,t,x,y,z,dhdx,dhdy,dhdz,volume)
exx=0.d0
eyy=0.d0
ezz=0.d0
exy=0.d0
eyz=0.d0
ezx=0.d0
do k=1,params%mpe
ic=icon(k)
exx=exx+dhdx(k)*unode(ic)
eyy=eyy+dhdy(k)*vnode(ic)
ezz=ezz+dhdz(k)*wnode(ic)
!exy=exy+(dhdx(k)*vnode(ic)+dhdy(k)*unode(ic))/2.d0
!eyz=eyz+(dhdy(k)*wnode(ic)+dhdz(k)*vnode(ic))/2.d0
!ezx=ezx+(dhdz(k)*unode(ic)+dhdx(k)*wnode(ic))/2.d0
enddo
! Change by Jean 13/2/08
! divergence should not be scaled by volume or weight
! div=(exx+eyy+ezz)*volume*w
div=exx+eyy+ezz
! end of change
pressure=div*penalty*eviscosity
deallocate (x)
deallocate (y)
deallocate (z)
deallocate (dhdx)
deallocate (dhdy)
deallocate (dhdz)
return
end
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|