CalculiX  2.13
A Free Software Three-Dimensional Structural Finite Element Program
Functions/Subroutines
mafillk.f File Reference

Go to the source code of this file.

Functions/Subroutines

subroutine mafillk (nef, ipnei, neifa, neiel, vfa, xxn, area, au, ad, jq, irow, nzs, b, vel, umfa, xlet, xle, gradkfa, xxi, body, volume, ielfa, lakonf, ifabou, nbody, neq, dtimef, velo, veloo, cvfa, hcfa, cvel, gradvel, xload, gamma, xrlfa, xxj, nactdohinv, a1, a2, a3, flux, nefa, nefb, iau6, xxni, xxnj, iturbulent)
 

Function/Subroutine Documentation

◆ mafillk()

subroutine mafillk ( integer, intent(in)  nef,
integer, dimension(*), intent(in)  ipnei,
integer, dimension(*), intent(in)  neifa,
integer, dimension(*), intent(in)  neiel,
real*8, dimension(0:7,*), intent(in)  vfa,
real*8, dimension(3,*), intent(in)  xxn,
real*8, dimension(*), intent(in)  area,
real*8, dimension(*), intent(inout)  au,
real*8, dimension(*), intent(inout)  ad,
integer, dimension(*), intent(in)  jq,
integer, dimension(*), intent(in)  irow,
integer, intent(in)  nzs,
real*8, dimension(neq), intent(inout)  b,
real*8, dimension(nef,0:7), intent(in)  vel,
real*8, dimension(*), intent(in)  umfa,
real*8, dimension(*), intent(in)  xlet,
real*8, dimension(*), intent(in)  xle,
real*8, dimension(3,*), intent(in)  gradkfa,
real*8, dimension(3,*), intent(in)  xxi,
real*8, dimension(0:3,*), intent(in)  body,
real*8, dimension(*), intent(in)  volume,
integer, dimension(4,*), intent(in)  ielfa,
character*8, dimension(*), intent(in)  lakonf,
integer, dimension(*), intent(in)  ifabou,
integer, intent(in)  nbody,
integer, intent(in)  neq,
real*8, intent(in)  dtimef,
real*8, dimension(nef,0:7), intent(in)  velo,
real*8, dimension(nef,0:7), intent(in)  veloo,
real*8, dimension(*), intent(in)  cvfa,
real*8, dimension(*), intent(in)  hcfa,
real*8, dimension(*), intent(in)  cvel,
real*8, dimension(3,3,*), intent(in)  gradvel,
real*8, dimension(2,*), intent(in)  xload,
real*8, dimension(*), intent(in)  gamma,
real*8, dimension(3,*), intent(in)  xrlfa,
real*8, dimension(3,*), intent(in)  xxj,
integer, dimension(*), intent(in)  nactdohinv,
real*8, intent(in)  a1,
real*8, intent(in)  a2,
real*8, intent(in)  a3,
real*8, dimension(*), intent(in)  flux,
integer, intent(in)  nefa,
integer, intent(in)  nefb,
integer, dimension(6,*)  iau6,
real*8, dimension(3,*)  xxni,
real*8, dimension(3,*)  xxnj,
integer, intent(in)  iturbulent 
)
25 !
26 ! filling the matrix for the conservation of energy
27 !
28  implicit none
29 !
30  logical knownflux
31 !
32  character*8 lakonf(*)
33 !
34  integer i,nef,indexf,ipnei(*),j,ifa,iel,neifa(*),
35  & neiel(*),jq(*),irow(*),nzs,ielfa(4,*),nefa,nefb,
36  & ipointer,ifabou(*),nbody,neq,indexb,nactdohinv(*),
37  & iau6(6,*),iturbulent
38 !
39  real*8 xflux,vfa(0:7,*),xxn(3,*),area(*),au(*),ad(*),b(neq),
40  & vel(nef,0:7),umfa(*),xlet(*),xle(*),coef,gradkfa(3,*),
41  & xxi(3,*),body(0:3,*),volume(*),dtimef,velo(nef,0:7),
42  & veloo(nef,0:7),rhovol,cvel(*),gradvel(3,3,*),sk,
43  & cvfa(*),hcfa(*),div,xload(2,*),gamma(*),xrlfa(3,*),
44  & xxj(3,*),a1,a2,a3,flux(*),xxnj(3,*),xxni(3,*),difcoef
45 !
46  intent(in) nef,ipnei,neifa,neiel,vfa,xxn,area,
47  & jq,irow,nzs,vel,umfa,xlet,xle,gradkfa,xxi,
48  & body,volume,ielfa,lakonf,ifabou,nbody,neq,
49  & dtimef,velo,veloo,cvfa,hcfa,cvel,gradvel,xload,gamma,xrlfa,
50  & xxj,nactdohinv,a1,a2,a3,flux,nefa,nefb,iturbulent
51 !
52  intent(inout) au,ad,b
53 !
54  if(iturbulent.eq.1) then
55 !
56 ! k-epsilon
57 !
58  sk=1.d0
59  else
60 !
61 ! k-omega
62 !
63  sk=0.5d0
64  endif
65 !
66  do i=nefa,nefb
67  do indexf=ipnei(i)+1,ipnei(i+1)
68 !
69 ! convection
70 !
71  ifa=neifa(indexf)
72  iel=neiel(indexf)
73  xflux=flux(indexf)
74 !
75  if(xflux.ge.0.d0) then
76 !
77 ! outflowing flux
78 !
79  ad(i)=ad(i)+xflux
80 !
81  b(i)=b(i)-gamma(ifa)*(vfa(6,ifa)-vel(i,6))*xflux
82 !
83  else
84  if(iel.gt.0) then
85 !
86 ! incoming flux from neighboring element
87 !
88  au(indexf)=au(indexf)+xflux
89 !
90  b(i)=b(i)-gamma(ifa)*(vfa(6,ifa)-vel(iel,6))*xflux
91 !
92  else
93 !
94 ! incoming flux through boundary
95 !
96  if(ielfa(2,ifa).lt.0) then
97  indexb=-ielfa(2,ifa)
98  if(((ifabou(indexb+1).ne.0).and.
99  & (ifabou(indexb+2).ne.0).and.
100  & (ifabou(indexb+3).ne.0)).or.
101  & (dabs(xflux).lt.1.d-10)) then
102  b(i)=b(i)-vfa(6,ifa)*xflux
103  else
104  write(*,*) '*ERROR in mafillk: the turbulent'
105  write(*,*) ' kinetic energy'
106  write(*,*) ' of an incoming flux'
107  write(*,*) ' through face ',
108  & indexf-ipnei(i),'of'
109  write(*,*)' element ',nactdohinv(i),
110  & ' is not given'
111  endif
112  else
113  write(*,*) '*ERROR in mafillk: the turbulent'
114  write(*,*) ' kinetic energy'
115  write(*,*) ' of an incoming flux'
116  write(*,*) ' through face ',
117  & indexf-ipnei(i),'of'
118  write(*,*)' element ',nactdohinv(i),
119  & ' is not given'
120  endif
121  endif
122  endif
123 !
124 ! diffusion
125 !
126  difcoef=umfa(ifa)+sk*vfa(5,ifa)*vfa(6,ifa)/vfa(7,ifa)
127 !
128  if(iel.ne.0) then
129 !
130 ! neighboring element
131 !
132  coef=difcoef*area(ifa)/xlet(indexf)
133  ad(i)=ad(i)+coef
134  au(indexf)=au(indexf)-coef
135 !
136 ! correction for non-orthogonal grid
137 !
138  b(i)=b(i)+difcoef*area(ifa)*
139  & (gradkfa(1,ifa)*xxnj(1,indexf)+
140  & gradkfa(2,ifa)*xxnj(2,indexf)+
141  & gradkfa(3,ifa)*xxnj(3,indexf))
142  else
143 !
144 ! boundary; either temperature given or adiabatic
145 ! or outlet
146 !
147  knownflux=.false.
148  ipointer=abs(ielfa(2,ifa))
149  if(ipointer.gt.0) then
150  if((ifabou(ipointer+5).ne.0).or.
151  & (ifabou(ipointer+1).ne.0).or.
152  & (ifabou(ipointer+2).ne.0).or.
153  & (ifabou(ipointer+3).ne.0)) then
154 !
155 ! no outlet:
156 ! (i.e. no wall || no sliding || at least one velocity given)
157 ! turbulent variable is assumed fixed
158 !
159  coef=difcoef*area(ifa)/xle(indexf)
160  ad(i)=ad(i)+coef
161  b(i)=b(i)+coef*vfa(6,ifa)
162  else
163 !
164 ! outlet: no diffusion
165 !
166  endif
167  endif
168 !
169 ! correction for non-orthogonal grid
170 !
171  if(.not.knownflux) then
172  b(i)=b(i)+difcoef*area(ifa)*
173  & (gradkfa(1,ifa)*xxni(1,indexf)+
174  & gradkfa(2,ifa)*xxni(2,indexf)+
175  & gradkfa(3,ifa)*xxni(3,indexf))
176  endif
177  endif
178  enddo
179 !
180 ! viscous dissipation
181 !
182  rhovol=vel(i,5)*volume(i)
183 !
184 ! sink terms are treated implicitly (lhs)
185 !
186  ad(i)=ad(i)+rhovol*0.09d0*vel(i,7)
187 !
188 ! source terms are treated explicitly (rhs)
189 !
190  b(i)=b(i)+rhovol*vel(i,6)*((
191  & (2.d0*(gradvel(1,1,i)**2+gradvel(2,2,i)**2+
192  & gradvel(3,3,i)**2)+
193  & (gradvel(1,2,i)+gradvel(2,1,i))**2+
194  & (gradvel(1,3,i)+gradvel(3,1,i))**2+
195  & (gradvel(2,3,i)+gradvel(3,2,i))**2))/vel(i,7)
196 c & -0.09d0*vel(i,7)
197  & )
198 !
199 ! transient term
200 !
201  b(i)=b(i)-(a2*velo(i,6)+a3*veloo(i,6))*rhovol
202  rhovol=a1*rhovol
203  ad(i)=ad(i)+rhovol
204 !
205  enddo
206 !
207  return
flux
subroutine flux(node1, node2, nodem, nelem, lakon, kon, ipkon, nactdog, identity, ielprop, prop, kflag, v, xflow, f, nodef, idirf, df, cp, R, rho, physcon, g, co, dvi, numf, vold, set, shcon, nshcon, rhcon, nrhcon, ntmat_, mi, ider, ttime, time, iaxial)
Definition: flux.f:24
Hosted by OpenAircraft.com, (Michigan UAV, LLC)