Civil Engineering Reference
In-Depth Information
!-----------------------derivative boundary conditions--------------------
IF(iel==1)THEN
det=x_coords(2)-x_coords(1)
kc(2,2)=kc(2,2)+uy*det/d6; kc(2,3)=kc(2,3)+uy*det/d12
kc(3,2)=kc(3,2)+uy*det/d12; kc(3,3)=kc(3,3)+uy*det/d6
ELSE IF(iel==nels)THEN
det=x_coords(2)-x_coords(1)
kc(1,1)=kc(1,1)+uy*det/d6; kc(1,4)=kc(1,4)+uy*det/d12
kc(4,1)=kc(4,1)+uy*det/d12; kc(4,4)=kc(4,4)+uy*det/d6
END IF; CALL fsparv(kv,kc,num,kdiag); CALL fsparv(bp,mm,num,kdiag)
END DO elements_2;f1=uy*det/(two*theta); f2=f1; bp=bp+kv; kv=bp-kv/theta
!-----------------------factorise equations-------------------------------
CALL sparin(bp,kdiag); loads=zero
!-----------------------time stepping loop--------------------------------
WRITE(11,'(a,i3,a)')" Time Concentration (node",nres,")"
WRITE(11,'(2e12.4)')0.0,loads(nres)
timesteps: DO j=1,nstep
time=j*dtim; CALL linmul_sky(kv,loads,ans,kdiag)
ans(neq)=ans(neq)+f1; ans(neq-1)=ans(neq-1)+f2
CALL spabac(bp,ans,kdiag); loads=ans
IF(nod==4.AND.j==ntime)CALL contour(loads,g_coord,g_num,nci,13)
IF(j/npri*npri==j)WRITE(11,'(2e12.4)')
&
time,loads(nres)*exp(-ux*g_coord(1,nres)/two/kay(1,1))*
&
exp(-uy*g_coord(2,nres)/two/kay(2,2))
END DO timesteps
STOP
END PROGRAM p87
New scalar reals:
d6
set to 6
.
0
0
f1
used to fix derivative boundary condition
f2
used to fix derivative boundary condition
pt25
set to 0
set to 12
.
d12
.
25
velocity in
x
-direction
ux
velocity in
y
-direction
uy
New dynamic real arrays:
ans
rhs vector in equilibrium equations
When convection terms are retained in the simplified flow equations, (2.132) has to
be solved. Again many techniques could be employed but, in the present program, an
implicit algorithm based on equation (3.94) is used. Thus this program is an extension of
Program 8.2.
When the transformation of equation (2.134) is employed, the equation to be solved
becomes,
u
2
2
2
2
c
x
∂
h
∂x
2
+
c
y
∂
h
v
∂h
∂t
−
+
h
=
(8.3)
2
∂y
4
c
x
4
c
y
distinguishes the process from a simple diffusion one.
However, reference to Table 2.1 shows that the semi-discretised “stiffness” matrix for this
problem will still be symmetrical, the
thus the extra term involving
h
h
term involving an element matrix of the “mass
matrix” type, namely
N
i
N
j
d
x
d
y
.
