1 / 17

# Artificial Dissipation - PowerPoint PPT Presentation

Artificial Dissipation. The artificial dissipation term used in FEM is different than in FD,. Note: A is included to match terms in other equations. Artificial Dissipation (cont’d 2). Conservation of Mass:. Conservation of Momentum. Artificial Dissipation (cont’d 3).

I am the owner, or an agent authorized to act on behalf of the owner, of the copyrighted work described.

## PowerPoint Slideshow about ' Artificial Dissipation' - jered

Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author.While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server.

- - - - - - - - - - - - - - - - - - - - - - - - - - E N D - - - - - - - - - - - - - - - - - - - - - - - - - -
Presentation Transcript

The artificial dissipation term used in FEM is different than in FD,

Note: A is included to match terms in other equations

Conservation of Mass:

Conservation of Momentum

Applying method of weighted residuals and integration by

parts to the artificial dissipation term in the conservation of

mass equation,

Conservation of mass

Conservation of momentum

Implementation of boundary conditions will be done in the

same way as for the FD equations.

program hw6c integer nodes(131,2), lnode(2) real kglobal(258,6),klocal(4,4),x(131),svec(258), real area(131),coord(2),phi(2),dphi(2), real larea(2),x(3),w(3),lu(2),lrho(2),lmu(2),mulocal real fvec(258),flocal(4),uvec(131),rhovec(131) real kbound1(2,2),kbound2(2,2)cc initialize variablesc maxnode=131 maxmat=2*(maxnode-2) maxel=130 nlast=10000 dt=1.e-3 dx=0.1 temp=300. gamma=5./3. xi(1)=0.7745966692414830 xi(2)=0.0 xi(3)=-x(1) w(1)=0.5555555555555560 w(2)=0.8888888888888890 w(3)=w(1)

Declare variables and

Initialize certain variables

c

x(i)=(i-1)*dx

c

if(x(i).le.3.0)then

c

area(i)=4-2./3.*x(i)

c

else

c

area(i)=0.6*(x(i)-3.)+2.

c

endif

c

10 continue

c

do 20 i=1,maxel

c

nodes(i,1)=i

nodes(i,2)=i+1

c

20 continue

Set up nodal coordinates, the area

at each node, and the connectivity

between local and global node numbers

c

do 40 j=1,6

c

do 30 i=1,maxmat

c

kglobal(i,j)=0.0

c

30 continue

c

40 continue

c

do 200 ie=1,maxel

c

lnode(1)=nodes(ie,1)

lnode(2)=nodes(ie,2)

lrho(1)=rhovec(nodes(ie,1))

lrho(2)=rhovec(nodes(ie,2))

lu(1)=uvec(nodes(ie,1))

lu(2)=uvec(nodes(ie,2))

coord(1)=x(nodes(ie,1))

coord(2)=x(nodes(ie,2))

larea(1)=area(nodes(ie,1))

larea(2)=area(nodes(ie,2))

lmu(1)=0.125*lrho(1)*abs(lu(1))*dx

lmu(2)=0.125*lrho(2)*abs(lu(2))*dx

hlocal=coord(2)-coord(1)

hlocalinv=1./hlocal

Initialize sparse matrix.

Loop through all the elements, one at a time. Set up local system variables at the two local nodes for the element of interest

c

do 100 i=1,2

c

iloc=2*i-1

jloc=2*j-1

c

do 60 jinit=1,4

c

do 50 iinit=1,4

c

klocal(iinit,jinit)=0.0

c

50 continue

c

60 continue

c

do 80 nint=1,3

c

xlocal=0.5*(xi(nint)*hlocal+(coord(2)+coord(1)))

phi(1)=hlocalinv*(coord(2)-xlocal)

phi(2)=hlocalinv*(xlocal-coord(1))

dphi(1)=-hlocalinv

dphi(2)=hlocalinv

alocal=larea(1)*phi(1)+larea(2)*phi(2)

dalocal=larea(1)*dphi(1)+larea(2)*dphi(2)

ulocal=lu(1)*phi(1)+lu(2)*phi(2)

rholocal=lrho(1)*phi(1)+lrho(2)*phi(2)

mulocal=lmu(1)*phi(1)+lmu(2)*phi(2)

Loop through i,j, and also

Set up indices of p,q.

Zero out local matrix

c

term=2.*alocal/dt*phi(i)*phi(j)-ulocal*

! alocal*dphi(i)*phi(j)+0.5*mulocal*(alocal*

! dphi(i)+phi(i)*dalocal)*dphi(j)

klocal(iloc,jloc)=klocal(iloc,jloc)+term*

! w(nint)*hlocal/2.

c

c k_ij^(12)

c

term=-rholocal*alocal*dphi(i)*phi(j)

klocal(iloc,jloc+1)=klocal(iloc,jloc+1)+

! term*w(nint)*hlocal/2.

c

c k_ij^(21)

c

term=alocal/dt*ulocal*phi(i)*phi(j)-0.5*

! ulocal**2*alocal*dphi(i)*phi(j)-temp/(2*

! gamma)*(alocal*dphi(i)+phi(i)*

! dalocal)*phi(j)

klocal(iloc+1,jloc)=klocal(iloc+1,jloc)+

! term*w(nint)*hlocal/2.

c

term=alocal/dt*rholocal*phi(i)*phi(j)-

! rholocal*ulocal*alocal*dphi(i)*phi(j)+0.5*

! mulocal*(alocal*dphi(i)+phi(i)*dalocal)*dphi(j)

klocal(iloc+1,jloc+1)=klocal(iloc+1,jloc+1)+

! term*w(nint)*hlocal/2.

c

c f_i^(1)

c

term=2.*rholocal*ulocal*alocal*dphi(i)-mulocal*

! rholocal*(alocal*dphi(i)+phi(i)*dalocal)

flocal(iloc)=flocal(iloc)+term*w(nint)*hlocal/2.

c

c f_i^(2)

c

term=rholocal*ulocal**2*alocal*dphi(i)+

! (temp/gamma*rholocal-mulocal*ulocal)*

! (alocal*dphi(i)+phi(i)*dalocal)

flocal(iloc+1)=flocal(iloc+1)+term*

! w(nint)*hlocal/2.

c

80 continue

c

100 continue

c

110 continue

c

if(ie.eq.1)then

c

do 130 jj=1,2

c

do 120 ii=1,2

c

kglobal(idfree+2+ii,jj+2)=kglobal(idfree+ii,jj+2)

! +klocal(ii+2,jj+2)

kbound1(ii,jj)=klocal(ii+2,jj)

c

120 continue

c

130 continue

c

else if(ie.eq.maxel)then

c

do 150 jj=1,2

c

do 140 ii=1,2

c

kglobal(idfree+ii,2+jj+2)=kglobal(idfree+ii,2+jj+2)+

! klocal(ii,jj+2)

kbound2(ii,jj)=klocal(ii,jj+2)

c

140 continue

150 continue

c

do 180 jj=1,4

c

do 170 ii=1,2

c

kglobal(idfree+ii,2+jj)=kglobal(idfree+ii,2+jj)+klocal(ii,jj)

kglobal(idfree+2+ii,jj)=kglobal(idfree+ii,2+jj)+klocal(ii+2,jj)

c

170 continue

c

180 continue

c

endif

c

200 continue

c

c Boundary Conditions

c

c Solution of Matrix

c

If not in first or

last element,

transfer entire

local matrix to

global matrix.

Boundary condition

and matrix solution

c

do 250 i=1,maxnode-2

c

i1=i*2-1

rhovec(i+1)=svec(i1)+rhovec(i+1)

uvec(i+1)=svec(i1+1)+uvec(i+1)

c

250 continue

c

uvec(1)=2.*uvec(2)-uvec(3)

uvec(maxnode)=2.*uvec(maxnode-1)-uvec(maxnode-2)

rhovec(1)=

rhovec(maxnode)=2*rhovec(maxnod-1)-rhovec(maxnode-2)

c

500 continue

c

stop

end

Note: index shift to

account for reduction

ff 1st and last set of

rows of matrix equation

Solution at first and last

node is given, but may

not be necessary.