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

Go to the source code of this file.

Functions/Subroutines

subroutine mafillo (nef, ipnei, neifa, neiel, vfa, xxn, area, au, ad, jq, irow, nzs, b, vel, umfa, xlet, xle, gradofa, 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, gradkel, gradoel)
 

Function/Subroutine Documentation

◆ mafillo()

subroutine mafillo ( 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)  gradofa,
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,
real*8, dimension(3,*), intent(in)  gradkel,
real*8, dimension(3,*), intent(in)  gradoel 
)
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,gradofa(3,*),
41  & xxi(3,*),body(0:3,*),volume(*),dtimef,velo(nef,0:7),
42  & veloo(nef,0:7),rhovol,cvel(*),gradvel(3,3,*),sw,be,ga,
43  & cvfa(*),hcfa(*),div,xload(2,*),gamma(*),xrlfa(3,*),
44  & xxj(3,*),a1,a2,a3,flux(*),xxnj(3,*),xxni(3,*),difcoef,
45  & gradkel(3,*),gradoel(3,*)
46 !
47  intent(in) nef,ipnei,neifa,neiel,vfa,xxn,area,gradkel,
48  & jq,irow,nzs,vel,umfa,xlet,xle,gradofa,xxi,gradoel,
49  & body,volume,ielfa,lakonf,ifabou,nbody,neq,
50  & dtimef,velo,veloo,cvfa,hcfa,cvel,gradvel,xload,gamma,xrlfa,
51  & xxj,nactdohinv,a1,a2,a3,flux,nefa,nefb,iturbulent
52 !
53  intent(inout) au,ad,b
54 !
55  if(iturbulent.eq.1) then
56 !
57 ! k-epsilon
58 !
59  sw=0.856d0
60  be=0.0828d0
61  ga=0.4404d0
62  else
63 !
64 ! k-omega
65 !
66  sw=0.5d0
67  be=0.075d0
68  ga=0.5532d0
69  endif
70 !
71  do i=nefa,nefb
72  do indexf=ipnei(i)+1,ipnei(i+1)
73 !
74 ! convection
75 !
76  ifa=neifa(indexf)
77  iel=neiel(indexf)
78  xflux=flux(indexf)
79 !
80  if(xflux.ge.0.d0) then
81 !
82 ! outflowing flux
83 !
84  ad(i)=ad(i)+xflux
85 !
86  b(i)=b(i)-gamma(ifa)*(vfa(7,ifa)-vel(i,7))*xflux
87 !
88  else
89  if(iel.gt.0) then
90 !
91 ! incoming flux from neighboring element
92 !
93  au(indexf)=au(indexf)+xflux
94 !
95  b(i)=b(i)-gamma(ifa)*(vfa(7,ifa)-vel(iel,7))*xflux
96 !
97  else
98 !
99 ! incoming flux through boundary
100 !
101  if(ielfa(2,ifa).lt.0) then
102  indexb=-ielfa(2,ifa)
103  if(((ifabou(indexb+1).ne.0).and.
104  & (ifabou(indexb+2).ne.0).and.
105  & (ifabou(indexb+3).ne.0)).or.
106  & (dabs(xflux).lt.1.d-10)) then
107  b(i)=b(i)-vfa(7,ifa)*xflux
108  else
109  write(*,*) '*ERROR in mafillo: the turbulence'
110  write(*,*) ' frequency'
111  write(*,*) ' of an incoming flux'
112  write(*,*) ' through face ',
113  & indexf-ipnei(i),'of'
114  write(*,*)' element ',nactdohinv(i),
115  & ' is not given'
116  endif
117  else
118  write(*,*) '*ERROR in mafillk: the turbulence'
119  write(*,*) ' frequency'
120  write(*,*) ' of an incoming flux'
121  write(*,*) ' through face ',
122  & indexf-ipnei(i),'of'
123  write(*,*)' element ',nactdohinv(i),
124  & ' is not given'
125  endif
126  endif
127  endif
128 !
129 ! diffusion
130 !
131  difcoef=umfa(ifa)+sw*vfa(5,ifa)*vfa(6,ifa)/vfa(7,ifa)
132 !
133  if(iel.ne.0) then
134 !
135 ! neighboring element
136 !
137  coef=difcoef*area(ifa)/xlet(indexf)
138  ad(i)=ad(i)+coef
139  au(indexf)=au(indexf)-coef
140 !
141 ! correction for non-orthogonal grid
142 !
143  b(i)=b(i)+difcoef*area(ifa)*
144  & (gradofa(1,ifa)*xxnj(1,indexf)+
145  & gradofa(2,ifa)*xxnj(2,indexf)+
146  & gradofa(3,ifa)*xxnj(3,indexf))
147  else
148 !
149 ! boundary; either temperature given or adiabatic
150 ! or outlet
151 !
152  knownflux=.false.
153  ipointer=abs(ielfa(2,ifa))
154  if(ipointer.gt.0) then
155  if((ifabou(ipointer+5).ne.0).or.
156  & (ifabou(ipointer+1).ne.0).or.
157  & (ifabou(ipointer+2).ne.0).or.
158  & (ifabou(ipointer+3).ne.0)) then
159 !
160 ! no outlet:
161 ! (i.e. no wall || no sliding || at least one velocity given)
162 ! turbulent variable is assumed fixed
163 !
164  coef=difcoef*area(ifa)/xle(indexf)
165  ad(i)=ad(i)+coef
166  b(i)=b(i)+coef*vfa(7,ifa)
167  else
168 !
169 ! outlet: no diffusion
170 !
171  endif
172  endif
173 !
174 ! correction for non-orthogonal grid
175 !
176  if(.not.knownflux) then
177  b(i)=b(i)+difcoef*area(ifa)*
178  & (gradofa(1,ifa)*xxni(1,indexf)+
179  & gradofa(2,ifa)*xxni(2,indexf)+
180  & gradofa(3,ifa)*xxni(3,indexf))
181  endif
182  endif
183  enddo
184 !
185 ! viscous dissipation
186 !
187  rhovol=vel(i,5)*volume(i)
188 !
189 ! sink terms are treated implicitly (lhs)
190 !
191  ad(i)=ad(i)+rhovol*2.d0*be*vel(i,7)
192 !
193 ! source terms are treated explicitly (rhs)
194 !
195  b(i)=b(i)+rhovol*((ga*
196  & (2.d0*(gradvel(1,1,i)**2+gradvel(2,2,i)**2+
197  & gradvel(3,3,i)**2)+
198  & (gradvel(1,2,i)+gradvel(2,1,i))**2+
199  & (gradvel(1,3,i)+gradvel(3,1,i))**2+
200  & (gradvel(2,3,i)+gradvel(3,2,i))**2))
201 c & -be*vel(i,7)*vel(i,7)
202  & +be*vel(i,7)*vel(i,7)
203  & )
204 !
205  if(iturbulent.eq.1) then
206  b(i)=b(i)+2.d0*rhovol*sw*
207  & (gradkel(1,i)*gradoel(1,i)+
208  & gradkel(2,i)*gradoel(2,i)+
209  & gradkel(3,i)*gradoel(3,i))/vel(i,7)
210  endif
211 !
212 ! transient term
213 !
214  b(i)=b(i)-(a2*velo(i,7)+a3*veloo(i,7))*rhovol
215  rhovol=a1*rhovol
216  ad(i)=ad(i)+rhovol
217 !
218  enddo
219 !
220  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)