Numeriska metoder för tjälanalys

V Lfnotat

Nummer: V99 Datum:

-1989-08

Titel: _{Numeriska metoder för tjälana lys}

Författare: Jan Ziobrowsk i

Avdelning: _{Vägavde lningen (M-sekt)}

Projektnummer: 4203901-6 _

Pråjektnamn: Tjäle och tjälskydd Uppdragsgivare: VTI

Distribution: fri / begränsad /

Statens väg- och trafikinstitut

W Väg- 06/1 TF.äfI/f& Pa: 58101 Linköping. Tel. 013-204000. Telex 50125 VTISGIS. Telefax 013-14 1436 I StitUtet Besök: Olaus Magnus väg 37 Linköping

NUMERISKA mmmm FÖR TJÃLANALYS av Jan Ziobrowski

VTI, Vägavdo, Linköping, Juni 1989

Arbetet med denna rapport har bedrivits under ledning av Lars Stenberg

:

I

L

FÖRORD

Sendan en lägre tid har bristen på ett generellt program.för

be-räkning av tjälförlopp i mark upplevts starkt. De vanligen före-kommande metoderna idag används till beräkning av tjäldjup eller isotermernas lägen i mark kring fjärrvärme och VA-ledningar.

Betydligt svårare är det att beräkna tjällyfning och isanrikningar

på olika nivåer i frusen jord. Tjällyftningen är ju beroende av

frysprocessens förlopp.

Det nu presenterade arbetet utgör en begränsad genomgång av presen-terade metoder. Ett arbete av det här slaget är mycket komplicerat och programmen är svårforcerade pga sin ofullständighet.

"God-bitarna" vi programmakarna behålla för sig själva. Det får vi

acceptera då fungerande program kan nyttjas gentemot konkurrenter

och dessutom är av stor ekonomisk betydelse.

Detta arbete har utförts under hård tidspress. Resultatet är mot den bakgrunden mycket gott. VTI har nu en god grund att bygga

vidare pâo

Målsättningen är att ta fram ett generellt och beskrivande program, där flertalet av de i tjälförloppet ingående fysikaliska proces-serna har kunnat modelleras matematiskt. Ett tvådimensionellt pro-gram. kan då nyttjas vid dimensionering av vägöverbyggnader,

ut-spetsningar kring trummor och kulvertar, skärningar och bankar,

jordarts- och jord/bergövergångar mm. Kopplingen till VTI:s frys-testutrustning är också en viktig ingrediens.

Projektarbetet har bedrivits inom ramen för VTI egen FoU-verksam-het, tjäle och tjälskydd.

8.

BILAGA 1 (XYFREZ programmet)

BILAGA 2 (FEHEAT programmet)

NNEHÃLLSFÖRTECKNING INLEDNING Metoder för problemlösning 91. Analytiska metoder .20 Approximativa metoder 63° Empiriska metoder v 94° Numeriska metoder TABELL l

Testkörning av program FEHEAT

Kort översikt över existerande algoritmer i tema tjällyftning

Tvådimensionell tjälmodell från CRREL

Indataparametrarna 01. Geometriska indata 02, Tidsparametrar 03. Värmeparametrar 04 Fuktighetsparametrar Slutledning LITTERATUR REFERENSER Sidan 10 10 10 11 11 13 från Bl/l tom 31/31 från BZ/l tom.B2/18

Motto:

"If the mathematical modell can be used to estimate the

likely range of frost heave

of candidate pavement profile, we have made a significant improvement in the procedure for designing roadway and air-field pavements in seasonal frost areas"

(CRREL, Report 80-10) 1. INLEDNING

Denna rapport är delresultat av ett utvecklingsarbete omfattande

probleminventering och analys av numeriska metoder, tillämpliga för problemlösning inom. ämnesområdet tjällyftning. Arbetet är baserat på sex utplockade referenser av tillgänglig litteratur med detta

tema. (Se kapitel LITTERATUR ref. från /1/ tom /6/).

Taber (1930) har skrivit att frysning och upptining resulterar i

stora skador på vägytan i områden med kallt klimat, men de styrande

fysikaliska processerna är inte tillräckligt klarlagda. Johnson och

Lovell (1953) räknade fram.13 yttre faktorer och 22 inre faktorer

sam påverkar tjällyftning. Det innebär att matematiska formler som

fullständigt beskriver tjällyftning måste se ut som.(1).

y = f (xl,x2,x3,...,x35) (1)

där

y = fysikalisk funktion som ska räknas ut

(till ex. djupet av nollisotermen), f = funktionssymbol,

Xi för i = 1,2,3...,35 = fysikaliska variabler som ska ingå i beräkningar (t.ex. en av klimatdata).

Flertalet av ingående faktorer kan inte heller ges en enkel karak-tär, som t.ex. struktur, porositet, snö osv. Därför skriver alla

författare i inledningen till sina rapporter att ämnet är mycket,

komplicerat och otillräckligt känt trots att tjällyftningsskador kostar mycket att reparera.

2 . rumpan FÖR Pmsmmösnme

Det finns olika metoder för analys av fysikaliska effekter som tjällyftning. Några av dem är baserade på att beräkna energibalans och värmeflöde. Alla metoder kan bli uppdelade i följande fyra grupper: 2,1. analytiska metoder, 2,2. approximativa metoder, 2,3. empiriska metoder, 204. numeriska metoder. 2.1. Anlytiska metoder

Analytiska metoder är av värde till lösning av enbegränsad klass

av problem för vilka noggranna ekvationer efter teoretisk analys, kan bli formulerade. Antalet problemställningar, för vilka ana»

lytiska metoder är tillämpliga, är få beroende på att bara enkla geometrier med ideala gränstillstånd kan beskrivas i formler. Tjälw lyftningen, som är berörd av 35 möjliga parametrar hör inte till de enkla problem,som kan bli lösta på analytiskt sätta Men analytiska

metoder kan bli tänkbara för att finna approximativa lösningar° 2920 Approximativa.metoder

Approximativa metoder kan förlängas till analytiska metoder. Om nogranna lösningar av problem.inte kan bli möjliga, då ärdet tänkbart att utföra approximativa beräkningar. Men ett ungefärligt värde kan användas om det inte behövs att veta ett riktigt noga

grannt värde,

2939 Empiriska metoder

Empiriska metoder är grundade på erfarenhet. På experimentell väg kan man finna empiriska formler som beskriver den fysikaliska verkligheten i probelemet. De empiriska formlerna ska istället lösas av de teoretiska ekvationerna samicke är tillgängligan Ett

representerade i analogimaskin av elektriska parametrar, likartade

dem i verkligheten.

Empiriska metoder är begränsade till situationer där empiriska formler kan bli framtagna. Empiriska ekvationer kan bli tillämpliga till mycket specifika problem.och kan inte bli extrapolerade för

att beskriva liknande problem" T.ex. ekvation 12) (se sidan 10 från

referens /5/).

U II 0.328 + 0.0578 F (2)

där

frostdjup / meter /,

köldmängd / grader-dagar Celsius/.

För att ta fram ekvation (2), samlades nödvändiga data från olika

platser i provinsen Ontario, Canada under 5 år. Regressiv analys av samlade data gav ekvation (2). Den är icke användbar till att be-räkna frostdjupet i någon annan provins än Ontario och således inte

i Norrland heller. För att skapa en sådan formel för t.ex. Norrland behövs det för en regressiv analys att använda mätdata, samlade i

toex. VTIs tjäldataarkiv, sam nu omfattar uppgifter från cirka 200 tjälgränsmätare. Observera också att en dylik empirisk ekvation

inte tar hänsyn till lokala förhållanden såsom jordlagerföljd och

grundvattenförhâllanden.

2ø4 Nümeriska.metoder

Numeriska metoder har breda applikationer i varje område av veten-skap och vårt liv. Problemlösning ska i korthet tillgå så att ut-gående från den fysikaliska verkligheten måste nödvändiga

ekva-tioner skapas. Dessa ska- sedan överföras till lämplig form och prepareras till numerisk lösning. De stora fördelarna med numeriska

metoder är deras möjlighet att modellera breda problemområden med behövlig noggranhetsgrad, också sådana med komplicerade gräns-tillstånd. Nackdelen med numeriska metoder är att de i sig själva är mycket komplicerade. En numerisk metods kriteria angående

innan metoden ska bli godkänd som.lämplig. Lyckligtvis kan datam

program användas till samma problemområden och användaren behöver endast veta en liten del om programmet och problemen för att utföra

beräkningen° I praktiken behövs det att utföra ett stort arbete ("tedious job ) för att samla nödvändiga indata för programmen° Som hjälpmedel för detta kan ibland dataugenerativa program.användas, sam i sig själva är dyra att skapa och svåra att använda (se sid. 19, 20 av referens /3/). För numeriska metoder ska också tas hänsyn

till följande:

- tillgänglighet av nödvändiga indata,

- samlingsbehov av indata,

- behov och tillgänglighet av maskintid.

Nackdelen med. numeriska metoder är att också vid enkla problem

måste dator användas. Sist men inte minst behöver man för numerisk

lösning ett kodat program.

2.4.1 Tillgång till program

För att spara tid och kostnader är det möjligt att söka lämpliga

program i publikationero Om det inte finns sådana program, då

kräver situationen att man grundligt studerar problemet och samlar

information om, existerande metoder för problemlösning. Samlade metoder ska bedömas med hänsyn till alla möjliga synvinklar för att plocka fram de lämpligaste och lovande lösningarna. Bearbetning av sådana metoder kan leda till en numerisk lösning. I ämnesområdet

tjällyftning finns det i samdade publikationer (se LiTTERATUR ref.

/2/, /3/, /5/) inte bara metodbeskrivning utan också text till

tillämpande programo Fullkomlig information om detta program är samlad i Tabell 10 Bara hälften av alla publikationer innehåller

programtexterna dvs. ref. /2/, /3/ och /5/. Två av dem (ref. /2/ och /5/) är av okänd anledning icke kompletta. Det enda som finns i sin helhet är FEHEAT programmet från Hanover, USA. Angående (ref. /2/) programmet EFIRING kan konstateras att i tillagd text fattas 7

subrutiner till EFTRING, trots att (se sid. 4 av LITTERATUR /2/)

T a b e l l 1: Sa mm an st äl ln in g av in fo rm at io n om ti ll gån g ti ll pr og ra m ba se ra d på sa ml ad e li tt er at ur st ud ie r. R e f e r e n s -n um m e r i L I T T E -R A T U R Or t, L a n d /1/ HHHHHHHH Ir vi ne , US A HHI-lHI-lI-GHl-l

I

-

--I

-I

--I

---I

--I St oc kh ol m, I S V E R I G E

I

-

-I

-I

-I

=

-I Ha no ve r, I U S A I Ha no ve r, I U S A

I-

---I

---I

-1

-I Ca rl et on , I C a n a d a I SV ER IG E I I Be näm ni ng av I sys te me t/ I pr og ra mm et I 19 82 I I I I 19 85 I I I I 19 83 I I I I 19 80 I I I I 19 81 I I I I 19 89 I I I

Är

pr

og

ra

mm

et

s

te xt bi la gd ? I

I

I-i

Är

pr

og

ra

m

I ko mple tt ? l

FR

OS

TZ

X

x=

A,

B,

c,

D

EF

TR

IN

G

FE

HE

AT

,

SS CO ND UC T F r o s t P r e -d i c t i o n I I I I I I I I I I I I I I NE J I I JA I HHHHHHHH JA I JA

"Huvudprogrammet och de subrutiner som detta är uppbyggt av redon

visas i sin helhet i bilaga 1." Ingen helhet, tyvärr, det fattas subrutiner: COORD, INBZFR, INM, INPUT2, PRINT3, SOIZFR och UTDZFR. Men kännedom om dessa fakta kan dyka upp vid försök att köra pro-gramförbindning med de subrutiner som redovisas i verkligheten. Då

är det försent att ta reda på. Bästa sättet att bedöma programmets användbarhet är testkörningo

3 .,

TESTKÖRNING av PROGRAM FEHEAT.

Programmet FEHEAT är sammansatt av följande rutiner:

FEHEAT (huvudrutin), FORMK, FRHS, ISOTHM, MCHB, TSM; Alla program»

texterna i sin helhet är tagna från publikation /3/. För att

an-passa program» FEHEAT till VAX/VMS omprogrammerades all

diskfil-hantering. Alla nya texterna är i sin helhet bifogade i Bilaga 2 till denna rapport och inskrivna i diskminnet i VAX/VMS. som.be= räkningsexempel för testkörning har tagits den som är i publikation /3/. Svårighetsgraden av dataförberedningen till FEHEAT (eller liknande program i tema tjällyfting) är stor. Trots att testexempw let utan förändringar överfördes till diskfilerna, har det uppstått

några datafelo De datafelen är så svåra att upptäcka att speciella

program, som hjälpmedel för detta har fått skrivas. En sådan utväg föreskriver också bruksanvisningen för FEHEAT (se side 39 av /3/). De är PRGl.SAS OCH PRG2.SAS, Erogramtexterna och ritningarna av geometriska indata till FEHEAT och beräknade av FEHEAT isotermer som PRG1.SAS och PRG2.SAS har ritat finns i diskfilerna på VTI:s dator VAX/VMB. På grund av överenstämmelse mellan resultatet (koor= dinater för isotermer) som redovisas i /3/ och resultat som har räknats av FEHEAT kan konstateras att VAX/VMS kopia av FEHEAT är den rätta, Den kan tas i drift vid behov. Användbarheten av FEHEAT är liten med hänsyn till begränsningarna. FEHEAT är förberedd att utföra beräkningar för enbart stabila tillstånd (dvs alla värme-flöden i beräknat område ska bli lika med noll). Sådan situation är

endast tänkbar t.eX. i norra Norrland med värmeledningar under väg

i ett område med permafrost. Den andra, SSCONDUCT program.i puba likation /3/ är sannerligen också. duglig för överföring till VAX/VMS miljö. Men det är ointressant på grund av följande:

- SSCONDUCT program. har lika funktioner som.FEHEAT. Enda skill-naden dem. emellan är att lösningsalgoritmen i SSCONDUCT är bam serad på finita differens metod och i FEHEAT program är den

algoritm av finita element metod,

- användbarhet hos de båda är lika litet.

Rekapitulerar man kapitel 2.491 och kapitel 3. kan det konstateras att det inte finns tillräcklig tillgång till lämpliga program i de studerade publikationerna inom.ämnésområdet tjällyftning. Då är det nödvändigt att ta i beaktande alla möjliga lösningsalgoritmer (inom

ämnesområdet) för att förbereda vägen till en egen lösning av

prob-lemet.

4 . Kom' övmzsnm: AV EXISTERANDE ALGORI'JMER I :rum TJÃLLLYFTNING Flödesprocessen för samtidig värme- och fukttransport med

fas-övergångar började Stefan att modellera för nästan 100 år sedan (1890). Stefans modell antar att:

s vatten fryser vid konstant temperatur (0°C), - det existerar en rörlig gränsyta mellan faserna.

Neuman har vidareutvecklat Stefans analys med partiella

diffe-rentialekvationer som. beskriver temperaturfördelningen på båda

sidor av den rörliga gränsytan i modellen. Neumans analys til-lämpades av Berggren (1943), Aldrich (1956), samt utvecklades

matematiskt av Goodman (1964), Sikarski och Boley (1965).

Om det finns ofruset vatten i den frusna delen då förekommer också

fukttransport. Williams (1972) har upptäckt att försumbarhet av

effekter till följd av fuktighetsflöde leder till grova fel.

De sista två decenierna har vittnat om att parallellt med utveck-lingen av datorer utvecklades också analysen av tjälproblem.

Nakano och Brown har modellerat (1971) värmeflödesprocessen. Harlan

(1973) har skapat en modell för kopplad värme- och fuktighetsflöde med hjälp av en finita differensbaserad numerisk algoritm. Jame

(1978) har testat Harlands modell med experimentella data. Alla de

modellerna är begränsde till rymdkoordinatssteg 3 om och tidsteg en halvtimme.

En metod för breddning av dessa modeller till stor skala i rum och tid för tvådimensionella tjälproblem.fanns inte är 1987.

En endimensionell modell är tillräcklig för ett stort antal av tillämpningar, men tvådimensionella modeller skulle behövas för att

toex. analysera fjärrvärmeledningar, vägbankar och jordbank på permafrost. Värmeflödesmodeller (dvs fuktflödeseffekten är för=

summad) med isotermisk fasövergång är rapporterade t.ex. av: Bafus

och Guymon (1976), Guymon och Hromadka (1977). Kontinuerligt de-formerad rutsystemsmodell presenterades av O'Niell och Lynch

(1981)° Den modellen är endimensionell och baserad på en finita element algoritmr Samma typ men mer avancerad tvådimensionell modell publicerades av Albert (1984). De båda modellerna kan inte tillämpas till tjälproblemanalys med fuktighetsflöde och icke homou gen jordart. Det behövs en tvådimensionell modell för tjälproblem-analys angående temperatur- och fuktfördelningen i marken som kan klara parametrarna från områden med kallt klimat till en fysikalisk

verklighet. Modellen ska bli lagom.kostsam att köra, dvs utan behov

av dataparametrar som är dyra att skapa. Sist men inte minst ska ett kodat program av denna modell också bli körbar på PC.

5 .

TVÅDDENSIONELL TJÃLMODELL FRÅN CRREL

Arbetet med denna modell påbörjades i Cold Regions Research and

Engineering Laboratory (CRREL, Hanover, New Hampshire, USA) är 1979. Delresultatet publiserades av Guymon (1980), Berg (1980),

Denna modell antar följande:

» Vattnets flöde sker i den ofrusna delen och drivs av hydraliska gradienter, som följer Darcys lag.

- Vattenflödet i den frusna delen är försumbar.

- Värmeflödet beskrivs av formel (1), sid. 45 av ref /1/).

- Vattenflöde beskrivs matematiskt av (3), (se sidan 46 av ref.

/1/)9

w Vattnets fasövergång kan approximeras som en isotermisk process. w Ofrusen del deformeras ej och vid tjälgränsen eller den frusna

delen är deformationen en följd av isbildning eller smältning.

m Portrycket vid tjälgränsen är beroende av kvarvarande vatten-mängd som ska bestämmas genom.frystest.

m Hystereseffekt är försumbar.

- Konstanta parametrar (t.ex. porositet) är konstanta med avseende

på tiden (dvs att smältning, frysning eller stelning förändrar

inte parametrarna).

a Frysning och smältning i tjälprocessen föregår på så sätt, att det inte uppstår några spänningar mellan de olika zonerna.

Tjälmodell med fasövergång, som inräknar fuktighets- och värme-flöde, är baserat på en enkel volumetrisk kontrollberäknings-algoritm. (Guymon och Hromadka, 1977) för isotermisk tjälprocess analys. Algoritmen sysslar med beräkningar av kontrollvolym.(litet ytelement), som också kallas nodular domän, skapad av en knut=

punktig integrationsmetod. Enligt den metoden antages att den frusna delen behåller temperaturen noll grader Celsius tilldess

att allt latent smältvärme har inräknats (dvs bortförts). Den

modellen är kallad Nodular Domän Integration (NDI) modell. NDI

algoritmen är beskriven i litteratur /1/ med matematisk bakgrund och konvergens, konsistens, nogrannhetsgrad faställning för NDI

metod. En stor del av detta är publicerat i tidigare publikationer

av samma autorer och står intemed i referens /1/. I /1/ finns

också PROTOO och FROSTZB programbeskrivning men programtexterna

10

5 . INDATAPARMTRARNA

601 Geometriska indata

Geometriska indata definierar beräkningsområdet. Det behövs hjälpa program, för att underlätta data-förberedningen och kontroll. Vid

finita differens metod används en rektangel som randdel (cell),

finita elementmetoden nyttjar triangeln som.litet yteelement.

No-dulär Domän .Integration (NDI) metod utför beräkningar för modulär domän (också kallad kontroll volym) som.är summan av trianglarnas

tredjedelar som har anknytning till den nodulara punkten (knutu punkten). Tre mediander delar triangelytan i tredjedelar. (Se sid.

8 av ref. /l/)o

602 Tidparametrar

De beskriver alla tider som ingår i beräkningarna:

starttid, stopptid, tidsteg, intervall på utskrifterna, tidsför= skjutning för indatavärdena som varierar sinusformigt. Tidspara-metrar behandlas inte i FEHEAT programmet (se litteratur /3/) som är baserat på "steady state" metod (dvs alla värmeflöden i modellen

är lika med noll)a

603 Vårmeparametrar

m volymetrisk värmekapacitet, vatten (J/m**3 K)

- volymetrisk värmekapacitet, is (J/m**3 K)

m värmekonduktivitet i vatten (W/m K)

m värmekonduktivitet i is (W/m K)

m latent smältvärme (W/m**3)

w undre frysgränsen i den använda frysmodellen

(oftast = O°C) (°C)

m värmekonduktivitet i ofrusen jord (W/m K)

= volymetrisk värmekapacitet i ofrusen jord (J/m**3 K)

(Obso NDI metoden behöver inte .värmekonduktiviteten och värmen kapaciteten i frusen jord, men EFTRING-PROGRAMMET (ref. /2/) be-höver dessa sam indata)°

11

694 Fuktighetsparametrar

= fuktighetskonduktivitet i vattenmättad jord w exponent i Gardners fuktighetskonduktivitetstal

n koefficient i Gardners fuktighetskonduktivitetstal

w exponent i fuktighetskonduktivitets korrektiv tal (sjunkningfaktor med hänsyn till ishalt).

w jordporositet

m exponent i Gardners vattenhaltstal w koefficient i Gardners vattenhaltstal = kvarstående vattenhalt i frusen del a porvattentryckhöjd

w vattenhalt i vattenmättad jord

7 o SLUWNING

Mot bakgrund baserad på information från den samlade litteraturen

och med stöd av genomförda studier av andra publicerade metoder kan

-man dra följande slutsatser. Det behövs en tvådimensionell modell

för tjälproblem.anlays angående temperatur- och fuktfördelningen i

marken. Det kan betraktas som lämpligt att börja med en

endimen-sionell modell baserad på analys av temperatur- och fuktfördelning

eller alternativt att börja med en tvådimensionell modell som.syss-lar med analys av enbart temperaturfördelning. Den teoretiska lösning som är mest lovande (av dem som.är kända) är Nodular Domän Integration modell från CRREL. NDI metoden inkluderar i sig tre olika metoder som är anropade i algoritmen på ett helt automatiskt sätt för att minimera beräkningsfelen i hela programmet, baserat på

NDI metoden. Det kan vidare behövas att bearbeta den teoretiska

bakgrunden sam beskrives i litteraturen /1/ för att nå en

beräk-ningsalgoritm, lämplig till programskrivning. För att möta löpande behov kan man emellertid använda XYFREZ (från CRREL) som.autori-serades och har tagits i drift på Vägavdelningen under april/maj 1989. XYEREZprogramtexterna, indataexempel och resultatutskrifterna

till dem. är i sin helhet bifogade i bilaga 1 till denna rapport. Det programmet kan lösa tvådimensionella värmeekvationer,

inklu-derande effekten av fasövergång, genom att beräkna värmekapacitet.

12

behöver utföras flera testkörningar för en fullständig verfikation av XYFREZprogrammet. Det fattas också nödvändiga hjälppgram.till XYFREZ för indataförberedning till hela XYFREZ och ritning av iso

termer baserade på resultatet av beräkningar. De hjälpprogrammen är

nödvändiga och ska skrivas. Parallellt ska utvecklingsarbete

ut-föras med NDI modellen för att nå möjlighet till fullständig

nume-risk tjälprohlemanalys angående temperatur och fuktfördelning i

marken. Om det framtida programmet ska omprövas och testköras med bra 'resultat, då kan man säga att det mål som stod i mottot till

8.

/l/

/2/

/3/

/4/

/5/

/6/

13 LITTERATURREFERENSER Ted Hromadka, June 1987,

CRREL, US Army Corps of Engineers, Special report 87-9

A nodal domain integration model of two-dbmensional heat and

soil-water flow coupled by soil-water change". Lars-Eric Jansson, 1985,

VBB Stockhobm,

Teoretisk simulering

"Tjälproblam vid fjärrvärmeledningar i gator".

MOR.° Albert and G. Phetteplace, April 1983,

CRREL, US Army Corps of Engineers, Report 83-10

Computer models

ductions". for two-dimensional steady stateheat con-ROL. Berg, GoL. Guymon and T.C. Johnson,

February 1980,

CRREL, US Army Corps of Engineers, Report 80010

"Mathematical model to correlate frost heave of pavements with

laboratory predictions".

R.A. Chisholm, W.A. Phang, September 1978,

Ministry of Transportation and Communication,

Ontario, Canada, Report

"Aspects of prolonged exposure of pavements to sub-zero temperatures".

D. Sheng, S. Knutsson, K. Axelsson,

Department of Civil Engineering,

Luleå University, Mars 1989,

Report

"Verification and application of a numerical model for frost front penetration".

O O O O O O O O O O O O 701 772 10 773 7777 11 0 0 0 0 7001 110 7002 120 7003 130 1 1 2 3 V Bilaga 1 Program XYFREZ.4.Public3

This program solves the ZnD heat equation, including phase change effects through the heat capacity.

*fakir* _{Main Program} intimt*

This routine obtains input for determining array dimensions.

The "real" main program follows in SUBROUTINE MAIN

PROGRAM XYFREZ_4

PARAMETER (ID=2000,KD=1000,ID2=4000,ID10=20000, ID100=200000,ID1000=2000000) COMMON NDBC(ID),XG(ID),YG(ID),R(ID),T(ID),TOLD(ID),TLAST(ID), NGLNBC(ID),ANS(ID),MTYPE(ID2),NBCFLG(KD),NBCNGL(KD), NODEL(ID10),DT(ID),A(ID100),TBC(IDlOOO), HETU(10),HETF(10),CNDU(10),CNDF(10),TF(10),EL(10) CHARACTER*6 IDOCl ' CHARACTER ANSWER WRITE(6,701)

FONMAT(///////,5X,' WELCOME TO XYFREZ,

VERSION 4 ',///)

WRITE(6,772)

FORMAT(///,' File for input is DATFIL .',/,

° ' Is it convinient for you ? Answer (Y/N):')

READ(5,7777)ANSWER IDOCl = 'DATFIL' IF(ANSWER.EQ.'Y') GO T0 11 IF(ANSWER.EQ.'N') GO T0 10 GO TO 9

WRITE(6,773)

FORMAT(///,' Enter name of input data file:',/)

READ(5,7777) IDOCl FORMAT (A)

OPEN(UNIT=1,FILE=IDOC1,STATUS='OLD') READ(1,*) NBND, NBC1IN, NELS, NTS, MEDIA WRITE(6,*) NBND, NBClIN, NELS, NTS, MEDIA IF(NBND.LE.ID) GO TO 110

WRITE46,7001)NBND,ID

FORMAT(///,' NBND value of',I8,' is greater than ID =',I8) GO TO 1000

CONTINUE

IF(NBClIN.LE.KD) GO TO 120 WRITE(6,7002)NBC1IN,KD

FORMAT(///,' NBNBClIN value of',I8,' is greater than KD =',I8)

GO TO 1000 CONTINUE

IF(NELS.LE.ID2) GO TO 130 WRITE(6,7003)NELS,ID2

FORMAT(///,' NELS value of',18,' is greater than IDZ =',I8)

GO TO 1000 CONTINUE

IF(NTS.LE.ID) GO TO 140 WRITE(6,7004)NTS,ID

7004 140 7005 1000 7999 U1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -3 0 0 0 0 0 0 772 Mi l-5 FORMAT(///,' GO TO 1000 CONTINUE

NTS value of',I8,' is greater than ID =',I8)

IF(MEDIA°LE910) GO TO 150

WRITE (6,7005)MEDIA

FORMAT(///,' MEDIA value of',I8, is greater than 10')

WRITE(6,7999)

FORMAT(////,' Execution HALTED !')

CLOSE(UNIT=1) STOP

CONTINUE

NTSP = NTS + 1

The following is a dummy value of NEW, for use in the dimension statement in subroutine MAIN

n

NBW = l

CALL MAIN(A, NBND, NEW, NDBC, XG, YG, R, T, TOLD, TLAST, NGLNBC, ANS, MTYPE, NBCFLG, NBClIN, NBCNGL, DT, NTSP, TBC, NODEL, NELS, MEDIA)

STOP E N D

SUBROUTINE MAIN(A, NBND, NBW, NDBC, XG, YG, R, T, TOLD, TLAST, NGLNBC, ANS, MTYPE, NBCFLG; NBClIN, NBCNGL, DT, NTSP, TBC, NODEL, NELS, MEDIA)

DIMENSION A(NBND,NBW), NDBC(NBND), XG(NBND), YG(NBND) DIMENSION R(NBND), T(NBND), TOLD(NBND), TLAST(NBND)

DIMENSION NGLNBC(NBND), ANS(NBND), MTYPE(NELS), NBCFLG(NBC1IN) DIMENSION NBCNGL(NBC1IN), DT(NTSP),TBC(NTSP,NBC1IN)

DIMENSION NODEL(NELS,3), TITLE(20), TF(10), EL(10) DIMENSION HETU(10), HETF(10), CNDU(10), CNDF(lO) CHARACTER*6 IDOCZ

CHARACTER*6 OUTFIL, PLTFIL CHARACTER ANSWER I n p u t NTS = NTSP - 1 READ(1,*) IO IF(IO.EQ.6) GO T0 12 WRITE(6,772)

FORMAT(///,' File for general output is RESULT',/, ' is it convinient for you ? Answer (Y/N):')

READ(5,7777)ANSWER

10 773 7777 11 12 1791 791 707 719 C C C 19 672 20 673 21 112 8001 09 O O O O O Q 95 OUTFIL='RESULT' IF(ANSWER.EQ.'Y') GO T0 11 IF(ANSWER.EQ.'N') GO TO 10 GO TO 9 WRITE(6,773)

FORMAT(///,' Enter file name for general output',///) READ(5,7777) OUTFIL

FORMAT(A)

OPEN(UNIT=2,FILE=OUTFIL,STATUS='NEW')

CONTINUE

READ(1,1791) (TITLE(K), K=l,20)

FORMAT(20A4)

WRITE(IO,791)(TITLE(K), K=l,20) FORMAT(////20A4///)

READ(1,*) ITRLM, KPRTVL

WRITE(IO,707) NBND, NBClIN, NELS, NTS, ITRLM, KPRTVL

FORMAT(/,10X,'NBND',5X,'NMBC1',6X,'NELS',7X,'NTS',

5x,'ITRLM',4x,'KPRTVL',/,4x;6r10)

READ(1,*)

IZDPLT, MSHPRT, KBVPRT

WRITE(IO,719) IZDPLT, MSHPRT, KBVPRT

FORMAT(//,21X,'IZDPLT',6X,'MSHPRT',6X,

'KBVPRT',/,15x,4I12)

IF(IZDPLT.EQ.O) GO TO 112

WRITE(6,672)

FORMAT(///,' File for plotting output is DATPLT.',/,

' Is it convinient for you ? Answer (Y/N):')

READ(5,7777)ANSWER IDOC2='DATPLT' IF(ANSWER°EQ°'Y') GO T0 21 IF(ANSWEROEQ.'N') GO T0 20 GO T0 19 WRITE(6,673)

FORMAT(///,' Enter file name for plotting output',///)

READ(5,7777) IDOCZ

OPEN(UNIT=33,FILE=IDOC2,STATUS='NEW') CONTINUE

WRITE(6,8001)

FORMAT(//,20X,' Preparing Finite Element Mesh...',//) READ(1,*) KPAUSE

WRITE(IO,709) KPAUSE

FORMAT(//,3OX,'KPAUSE =',IZ)

The next loop reads each node, the type bc at the node and the ceordinates of the node.. Also sets bc flags. NMBCI O KBCZ 0 DO 100 NODE = 1,NBND READ(1,*) ND, NDBCV, XGV, YGV IF( NDBCVOEQ.2 ) KBC2 = 1 NGLNBC( ND ) = 0 IF ( NDBCV.NE01 ) GO T0 95 NMBCI = NMBCl + 1 NGLNBC ( ND ) = NMBCl NBCNGL ( NMBCl ) = ND NDBC( ND ) = NDBCV

81/ 3

XG( ND ) = XGV YG( ND ) = YGV 100 CONTINUE

IF(NMBC10EQ.NBC1IN) GO TO 111

WRITE(6,7770) NBClIN, NMBCl

7770 FORMAT(///,10X,' NMBCl value of',I6,'input is different from' 1 /,10X,' actual value of',IlO,

2

///////,20X,' EXECUTION HALTED!')

CLOSE(UNIT=1) CLOSE(UNIT=33) CLOSE(UNIT=IO) STOP lll CONTINUE C C

IF(MSHPRT.EQ.1)

1 WRITE(IO,701) ( K, NDBC(K), XG(K), YG(K), K = 1,NBND)

701

FORMAT(//,(12Xy'Node',4X,' BC ',12x,'x',

1 10X,'Y'),/,(9X,ZI7,3X,1P2E1162) )

The following loop reads the connectivity of the nodes and detemmines the material type of the elements

O O O O O DO 110 NLMI =1,NELS READ(1,*) NL, LTYPE, Nl, N2, N3 MTYPE(NLMT) LTYPE NODEL(NL,1) Nl NODEL(NL,2) NZ NODEL(NL,3) N3 lO CONTINUE

Check for flipped element : is X1*X2 - Yl*Y2 .LT. O ? Also detenmine NDG and NBW

O O O O O O H NNDG = 0 DO 250 L KNTERR = 249 Nl NODEL(L,1) N2 NODEL(L,2) N3 - NODEL(L,3) NLBCZ = NDBC(Nl) + NDBC(NZ) + NDBC(N3) IF(NLBC2.NE°6) GO TO 2249 WRITE(IÖ,7249) L = 1,NELS 0

7249 FORMAT(////,' Element no',15,' has three type 2 BCs',//, 1 ' FIX ITE')

GO TO 1111 c

2249 NDIF = IABS(N2 - Nl)

IF( NDIF.GT°NNDG) NNDG = NDIF NDIF : IABS(N3 « NZ) ' IF( NDIFOGTONNDG) NNDG = NDIF NDIF = IABS(N1 « N3)

IF( NDIF.GT.NNDG) NNDG = NDIF

X1 = XG(N2) - XG(N1) X2 = XG(N3) - XG(N1) Yl = YG(N2) n YG(N1) YZ = YG(N3) m YG(N1) PROD = X1*Y2 - Y1*X2 IF(PRODQGT.O.O) GO TO 250 KNTERR : KNTERR + 1

IF(KNTERR.GT.1) THEN

WRITE(IO,7111) L 81/ 5 7111 FORMAT(Å///,' Incurable mesh problem.in element no',I6)

GO TO 1111 ELSE ENDIF NHOLD = NODEL(L,3) NODEL(L,3) = NODEL(L,2) NODEL(LYZ) = NHOLD GO TO 249 250 CONTINUE C

C Arguments to be passed for later calculations involving the matrix A,

C eøgo NBW is the bandwidth of the handed matrix.

C NBW = 2 * NNDG + 1 NDG = NNDG + 1 O Q O IF ( MSHPRT.NE.1 ) GO TO 1250 WRITE(IO,702) (K,MIYPE(K),(NOQEL(K,J),J=1,3),K=1,NELS) 702 FORMAT(//,2(6x,'Element',3X,'Material',6X,' Nodes',2X),/, 1 2(9X,I4,5X,I2,6X,314) ) 1250 CONTINUE C

C Write mesh plotting info on disk

C

IF(IZDPLT.EQ.O) GO TO 1108

WRITE(33,7010) NELS

7010 FORMAT(I7)

DO 107 NL = 1,NELS

WRITE(33,7011) NL, MIYPE(NL), (NODEL(NL,J),J=1,3)

7011

FORMAT(5I7)

107

CONTINUE

WRITE(33,7010) NBND

DO 108 ND = 1,NBND

WRITE(33,7012) ND, NDBC(ND), XG(ND), YG(ND)

7012

FORMAT(217,1P2E14.5)

108 CONTINUE

1108 CONTINUE

C Input thermal properties DO 1109 I = l,MEDIA

READ(1,*) HETF(I), HETU(I), CNDF(I), CNDU(I), EL(I), TF(I)

WRITE(IO,705) I,HETF(I), HETU(I), CNDF(I), CNDU(I), EL(I), TF(I)

705 FORMAT(///,15X,'Thermal Properties of Material',13,//,9X,

1

'HETF',6X,'HETU',6X,'CNDF',6X,'CNDU',8X,'EL',8X,'TF',

2 3X,/,(3X,1P6E10.2),2X) 1109 CONTINUE

C

C THETA and AMP control implicitness in time stepping and iteration

C

READ(1,*) TIMLIM, AMP, THETA WRITE(IO,706) TIMLIM, AMP, THETA

706 FORMAT(////,15X,'Solution Control Parameters',//,21X, l 'TIMLIM',9X,'AMP',7X,'THETA',/,15" 'PiEL:.3) C C Determine IC s C READ(l,*) KIC IF(KIC.EQ.1) GO T0 60 READ(1,*) TEMPIC C C Initialize TOLD C DO 50 ND = 1,NBND TOLD(ND) = TEMPIC

50 60 70 99 733 0 0 0 0 0 0 0 0 0 0 0 0 0 0 301 0 0 0 305 310 320 w0 0 0 0 0 0 l CONTINUE GO T0 99 DO 70 ND = 1,NBND

READ(1,*) NODE, TEMPIC TOLD(NODE) = TEMPIC CONTINUE

CONTINUE

WRITE(IO,733) (K,TOLD(K), K = 1,NBND)

FORMAT(/////,15X,'Initial Temperatures at Each Node',//, 3(IlO,lP1E12.3) )

Set BC's and time step values

READ(1,*) KBCT

For KBCT = 1, bc values are constant over time, with

KDT 1 for uniform time steps, or

KDT 2 for uniform steps in SQRT(TIME)

For KBCT = 2 bc values over time and time steps are input IF(KBCT.EQ°2) GO TO 360

Constant BC's, initialize the RHS of the nonlinear system at bc's

DO 301 ND = 1,NMBC1

READ(1,*) NODE, TBCVAL

R(NODE) = TBCVAL

CONTINUE

READ(1,*) KDT

IF(KDT.EQ.2) GO TO 310

Unifonm time steps, calculate DT(nts) DTCNST = TIMLIM / NTS

DO 305 IDT = l,NTS DT(IDT) = DTCNST

GO TO 399

Unifomm steps in SQRT(TIME), calculate DT(nts) DTSQT SQRT(TIMLIM) / NTS DT(1) DTSQT**2 TIM2 = DT(1) DO 320 IDT = 2,NTS TIMl = TIM2 TSQT = SQRT(TIM2) TNWSQT = TSQT + DTSQT TIM2 = TNWSQT**2 DT(IDT) = TIM2 " TIMl CONTINUE

GO TO 399

Read in bc values and associated times, with number of intervening time steps. Set flags and interpolate

between nodes.

CONTINUE

IT = 0

TIMOLD = 0.0

DO 390 IBV = 1,10000

READ(1,*) TIMBD, NUMCHG, NTSBD

365 366 370 371 374 377 378 380 390 C 399 741 742 401 440 1743 743

DTIM = ( TIMBD = TIMOLD ) / NTSBD TIMOLD = TIMBD DO 365 ND = 1,NMBC1 NBCFLG(ND) = -1 CONTINUE Nl = 1 DO 380 ND = 1,NMBC1 IF (ND.GT.NUMCHG) GO TO 366 READ(1,*) NODE,TBCVAL TBCVLP = TOLD(NODE) NM = NGLNBC(NODE) IF(IT.NE.O) TBCVLP = TBC(IT,NM) DTMP = (TBCVAL - TBCVLP) / NTSBD GO TO 377 CONTINUE DO 370 NM = Nl, NMBCl IF(NBCFLG(NM).LT.0) GO TO 371 CONTINUE Nl = NM + 1 NODE = NBCNGL(NM) IF(IT9EQ°O) GO TO 374 TBCVLP = TBC(IT,NM) ITM = IT m 1 IF(IT0EQ.1) TBCVPP TOLD(NODE) IF(IT.GT.1) TBCVPP - TBC(ITM,NM) DTMP = (TBCVLP * TBCVPP) * DTIM / DT(IT)

GO TO 377

TBCVLP = TOLD(NODE)

DTMP = 0.0

NBCFLG(NM) = 1

DO 378 :IE = 1,NTSBD

ITC = IT + ITB DT(ITC) = DTIM TBC(ITC,NM) = TBCVLP + ITB * DTMP CONTINUE CONTINUE IT = IT + NTSBD IF(IT.GE.NTS) GO TO 399 CONTINUE CONTINUE IF(KBVPRT,EQ02) GO TO 440 IF(KBVPRT.NE.1) GO TO 499

WRITE(IO,741)

FORMAT(////,15X,'Time and DT values at Each Time Step',//, 11X,'TIME STEP',8X,'TIME',10X,'DT')

TIMCK = 000

DO 401 IT = 1,NTS

TIMCK = TIMCK + DT(IT)

WRITE(IO,742) IT, TIMCK, DT(IT)

FORMAT(I20,1P2E12.3)

CONTINUE

GO TO 499 CONTINUE

WRITE(IO,1743)

FORMAT(////,15X,'Time and DT at Each Time Step, with ',/, 20X,'Boundary Values at Each Node')

TIMCK = 090

DO 460 IT = 1,NTS

TIMCK = TIMCK + DT(IT)

WRITE(IO,743) IT, TIMCK, DT(IT)

FORMAT(/////,10X,'TIME STEP =',I4,5X,'TIME =',lPlE12.3,5X, 'DT =',1P1E12.3)

745 744 460 499 0 0 0 0 7744 44 997 7998 8000 999 O O O O W O O O O O O Q O Q O 1 WRITE(IO,745)

FORMAT(/,20X,'Nodes and their temperatures') IF(KBCT°EQ°2)

WRITE(IO,744) (NBCNGL(J), TBC(IT,J), J=1,NMBC1)

FORMAT(3(I10,1?1E12°3)) IF(KBCT.EQ°1) WRITE(IO,744) CONTINUE (NBCNGL(J), R(NBCNGL(J)),J = 1,NMBC1) CONTINUE

Input convective BC stuff

IF(KBC2°EQ°0) GO TO 44

READ(1,*) HH, TA

WRITE(IO,7744) HH, TA

F0RMAT(//,20x,'HH =',1P1E12°3,4X,'TA =',1PlE1203)

CONTINUE

IF(KBAUSE:EQ.0) GO TO 998 WRITE(6,7998)

FORMAT(///////,4X,'OK to proceed ( Y = YES

READ(5,7777) ANSWER WRITE(6,8000)

FORMAT(////////,10X, Computing Finite Element Solution...',////)

N = NO )')

IF(ANSWER°EQ.'Y') Go TO 998

IF(ANSWER6NE. N5) GO TO 997 WRITE(6,7999)

FORMAT(//////////,10X,'User mandated STOP E') Close all files and stop if user has so commanded CLOSE(UNIT=1)

CLOSE(UNIT=33) CLOSE(UNIT=IO) STOP

CONTINUE

Initialize time loop parameters TIME g 000 THETAM = THETA m 100 AMPM = 1.0 AMP INEW = 0 ITRLMP = ITRLM + 1 KTPRNT = 0 KTPLOT = 1 Adjust dt and tbc NTSP = NTS + 1 DT(NTSP) = DT(NTS) DO 999 NM = 1, NMBCl TBC(NTSP,NM) = TBC(NTS,NM) CONTINUE

IF(IZDPLT.NE.O) WRITE(33,7099) IZDPLT

IF(IZDPLT.EQ.0) IZDPLT = -1 '

FORMAT(I5)

Initialize ANS, TLAST

989 7001 0 O O O O 920 930 W O O O O O O O 0 0 0 0 0 0 716 713 990 800 C C

DO 989 ND = 1,NBND

ANS(ND) : TOLD(ND)

TLAST(ND) = TOLD(ND)

CONTINUE

WRITE(IO,7001)

FORMAT(5(/////),20X,'Calculated Temperature Solution',//////)

Time Loop

D0 1000 IT = 1, NTSP

TIMOL = TIME

TIME TIME + DT(IT) ITER KTPRNT KTPLOT II I I O KTPRNT + 1 KTPLOT + 1 IF(KBCT.EQ.1) GO TO 930 DO 920 NM = 1, NMBCl NODE = NBCNGL(NM) R(NODE) = TBC(IT,NM) CONTINUE CONTINUE

Call the routine which performs element by element integrations and returns coefficient matrix A and r.h.s. vector R

CALL TRINT(A,NBND,NBW,NDBC,NODEL,NELS,XG,YG,IO,NDG,DT(IT),R,

CNDF,CNDU,HETF,HETU,TF,THETA,THETAM,T,TOLD,TLAST,EL,ITER,

ITRLMP,TIMDL,AMP,AMPM,IZDPLT,KTPLOT,ANS,IT,MEDIA,MTYPE,HH,TA) .

Solve algebraic system of equns

CALL BANSAL(A,NBW,NBND,NBND,R,ANS,INEW,KERR,IO) IF(ITER.EQ.ITRLM) GO TO 950 ITER = ITER + 1 G0 TO 901 CONTINUE IF(KTPRNT°NE.KPRTVL) KTPRNT = 0 GO TO 990 Print results WRITE(IO,716) TIME, IT

FORMAT(//,25X,'TIME =',1P1E12.3,4X,'TIME STEP =',I6) WRITE(IO,713) (K, ANS(K), K = 1, NBND)

FORMAT(//,3(6X,'NODE',8X,'TEMP')r/,3(5X,IS,1P1E12.3)) CONTINUE

Update values of TLAST and TOLD DO 800 ND = 1, NBND TLAST(ND) = TOLD(ND) TOLD(ND) = ANS(ND) 1000 CONTINUE C FIL = 1.0 VRYNEG = _1.1E3O

IF(IZDPLT.GT.O) WRITE(33,7007) VRYNEG,FIL,FIL,FIL

7007 1111 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 100 FORMAT(1P4E11.2) 81/ 10

Close all files and stop CONTINUE CLOSE(UNIT=1) CLOSE(UNIT=33) CLOSE(UNIT=IO) RETURN E N D SUBROUTINE TRINT(A,NBND,NBW,NDBC,NODEL,NELS,XG,YG,IO,NDG,DT,R, CNDF,CNDU,HETF,HETU,TF,THETA,THETAM,T,TOLD,TLAST,EL,ITER, ITRLMP,TIMOL,AMP,AMPM,IZDPLT,KTPLOT,ANS,IT,MEDIA,MTYPE,HH,TA)

Subroutine to perform.finite element integrations and assemble

handed global matrix for linear triangles

COMMON/SMALL/ ZERO

DIMENSION A(NBND,NBW), NDBC(NBND), NODEL(NELS,3), XG(NBND)

DIMENSION YG(NBND), X(3), Y(3), R(NBND), T(NBND), TOLD(NBND) DIMENSION TLAST(NBND), ELA(3), ELB(3), TE(3), ANS(NBND), EL(10) DIMENSION TF(10), HETU(10), HETF(10)

DIMENSION CNDF(10), CNDU(10), MEYPE(NDLS)

Identifies a first pass IF(IT * (ITER + 1) .NEO 1) G0 TO 101

The next computations eliminate the problems that could start if the maximum initial temperature was zero

TMAX TOLD(1) IMIN TOLD(1) DO 100 ND = 2,NBND

TM = TOLD(ND)

IF(TM .GTO TMAX) TMAX = TM IF(TM .LTo DMIN) TMIN = TM CONIINUE

DELTF = 000025 * ( TMAX m TMIN ) TMAX = ABS( TMAX)

IF(ABS(TMIN).GT°TMAX) TMAX = ABS(TMIN)

IF(DELTFOLT.050025 * TMAX) DELTF = 0.0025 * TMAX CONTINUE

Update the T(ND) used to calculate th m * 1: I

DO 105 ND = 1, NBND

T(ND) = AMP * ANS(ND) + Amen * TOLD(ND)

CONTINUE

DO 10 ND = 1, NBND

O O O O O l -* k o O O O O O O O O O O O O O O 0 0 0 C C C11009 C O O O O O O OQ O D0 5 J = 1, NBW 81/ 11

A(ND,J) = 0.0

CONTINUE

IF(NDBC(ND).NE.1) GO TO 9

A(ND,NDG) = 1.0

GO T0 10

R(ND) = 0.0

CONTINUE

Element loop. This is the one that does the quadrature,

and assembles A and R. DO 1000 NLMT = 1,NELS

Initialize flags

KPHS = "3 : Element frozen KEHS = +3 : Element unfrozen KFRZ = 1 : Phase change element KFRZ = 3

KPHS = 0 RF = 100

Identify material type of this element

Il = MTYPE(NLMT)

'

The next loop works with each element to determine if it is

a phase change element or not and then does necessery computations to enable calculation of phase change effect.

DO 400 K = 1,3 NK = NODEL(NLMT,K) WRITE(IO,11009)NLMT,NELS,NK,NLMT,K,NODEL(NLMT,K) FORMAT(/,'NLMT,NELS,NK,NLMT,K,NODEL(NLMT,K)',6I6) X (K) XG (NK) Y (K) YG (NK) TE(K) = T(NK)

This section accounts for the possibility that the node

temperature equals the phase change temperature. If

it does, it is shifted in the direction of its last

recorded value.

TFMZ = TF(Il) n DELTF TFPZ = TF(I1) + DELTF

IF ((TE(K).GT.TFMZ)°AND.(TE(K).LT.TFPZ)) THEN IF (TLAST(NK).GT.TF(Il)) THEN

TE(K) = TF(Il) + DELTF ELSE

TE(K) = TF(Il) w DELTF

ENDIF '

ENDIF

Record frozen or unfrozen status of node

KPHS - l KPHS + 1 IF(TE(K).LT.TF(Il)) KPHS

IF(TE(K).GT.TF(I1)) KPHS

Identify next node around the element KP = K + 1

O Q O O Q O 360 0 0 ₀ O O O G O ÖO O O O O O O O O O O O D ) O O O G O O Q O O O O 390 11011 Bypass NE = NODEL(NLMT,KP) _{31/ 12}

the rest of the loop if no phase change

IF((

(T(NP).GT5TF(Il)).AND.(T(NK).GT°TF(I1))).O

((T(

NP)°LT°TF(II))OAND°(T(NK).LT.TF(II)

))) G

R.O TO 400 IF(KFRZ.EQ.3) NPL 2 NP IF(KFRZQEQ°3) GO TO 380 NDF = NPL NK IF(NDF.NE.0) GO TO 360

( XC,YC ) is vertex of triangular subwelement

NDl = NK XC = X(K) YC = Y(K) GO TO 380 NDl 2 NP XC = X(KP) "YC = Y(KP) CONTINUE KFRZ = KFRZ m 1

The next computations are for the phase change element, ccmpnting the cocrdinates ef the phase change isetherm,

change

and the value of the basis/waight function at the

iso-thenm ( ELA and ELB ).

The next few lines interpolate to find the phase isothetm.

WRITE(IO,11011)NP,NK,II,T(NP),T(NK),TF(Il)

FORMAT(/,'NP =',I5,'

NK =',I5,' 11 =',I5,/,

'

T(NP) =',191312.2,'

T(NK) =',1P1E12.2,

'

TF(Il) =',1PlE12.3)

ROB = ABS (T(NP) T(NK))

IF(ROB°LT.OolE°30) CLOSE(UNIT=IÖ) ROB = ABS (T(NP)°T(NK))

ROBl = ABS (TF(Il)-T(NK)) IF(ROBl.LTBOQlE=30) THEN ROBl = 000 DLTF = 0.0 ELSE DLTF = ABS((TF(Il) T(NK)) / (T(NP) - T(NK))) ENDIF DLTF = ABS((TF(Il) == T(NK)) / (T(NP) -° T(NK))) DELX = DLTF * (XG(NP) XG(NK))

DELY = DLTF * (YG(NP) * YG(UV * IF(KFRZ°EQ01) GO TO 390 KB = XG(NK) + DELX YB = YG(NK) + DELY ELB(KP) = DLTF ELB(K) = 100 ° DLTF KPP = KP + 1 IF( KPPoEQ°4 ) KPP = 1 ELB(KPP) = 0.0 GO TO 400 XA = XG(NK) + DELX YA = YG(NK) + DELY

400 102 0 0 0 410 498 499 0 0 0 0 0 50 O O O O O O h EEA(KP) = DLTF EEA(K) = 1.0 - DLTF KPP = KP + 1 IF(KPP.EQ°4) KPP = 1 ELA(KPP) = 00 CONTINUE

81/ 13

Update tlast. The The DO 102 ND = 1, NBND TLAST(ND) = T(ND) CONTINUE IF(KPHS.NE03) GO TO 410

next assignments are for an unfrozen element. HET = HETU(I1)

COND = CNDU(I1) GO TO 499

IF(KPHS.NE°°3) GO TO 420

next assignments are°for an frozen element. HET = HETF(I1)

COND = CNDF(Il) GO TO 499

Calculate area of triangular sub-element

XB*YC ° XC*YB XA*YC + XC*YA + XA*YB - XB*YA ABS(A1/2.0)

A1 A1

Assigne frozen and unfroszen conductivities to appropriate parts of element.

COND = 0.0 VKl = CNDU(Il) VKZ = CNDF(Il) IF(T(ND1).GT°TF(I1)) GO TO 498 VKl = CNDF(I1) VK2 = CNDU(Il) CONTINUE CONTINUE

Calculate the area of the element using cross product.

AE = X(2)*Y(3) » X(3)*Y(2) - X(1)*Y(3) + X(3)*Y(1) + X(1)*Y(2) - X(2)*Y(1) AE = AE/2.0

Assigne conductivity factor for non-phase change element.

C4A = 0.25*COND / AE IF(KFRZ°EQ.1) GO TO 450

BCF = HET * AE / ( 12.0 * DT \ GO TO 4450

CONTINUE

This next assignment accounts for the different sensible heat in the unfrozen and frozen portion of the phase change element.

IF( T(ND1).GT°TF(I1) ) C1 = HETU(Il) C2 = HETF(I1)

ELSE C1 C2 END IF

BCF = ( C1*A1 + C2 * (AE Al)) / (1200 * DT)

HETF(I1) 81/ 14

HETU(Il)

Calculaticns on element with phase change to account for the

latent heat effect discontinuity

Calculate length of tau

TAU = SQRT( (XB w XA)**2 + (YB w YA)**2)

The magnitude cf the gradient cf T is computed using the X and Y

derivatives of the equation of the plane T(X,Y) = A + BX + CY, ice° , expressing A and B in terms of X(I), Y(I), & T(I) using

CramerÅs rule. The 2*AE is the determinant used in Cramer's rule

0 0 0 0 0 0 0 0 0 0 0 0 0 0

SUMDX

TE(1)*(Y(2)WY(3))

TE(3)*(Y(l)-Y(2))

TE(3)*(X(2)uX(l))

TE(1)*(X(3)GX(2))

SQRT( SUMDX**2 + SUMDY**2) / (2.0 * AE)

= EL(I1) / ( DTDS * DT )

+ TE(2)*(Y(3)°Y(1))

SUMDY + TE(2)*(X(1)X(3)) H + I l + l l DTDS = ELDSDT

Stering zodo isothenm locations for plctting

0

0

0

IF(KTPLOT.EQ.IZDPLT.AND.ITER.EQ°O) THEN WRITE(33,7001) XA, YA, XB, YB

7001 FORMAT(1P4E13.4) END IF

4450 CONTINUE

Element node row loop. DO 600 I = 1,3 0 0 0 0

IG = N0DDL<NLMD,I)

IF(NDBC(IG)0EQ°1) GO TO 600

Treat convective BC'sc 0 0 0

IF(NDBC(IG)°NE02) G0 TO 599

DO 510 J = 1,3

IF(J°EQ°I) GO TO 510

JG = NODEL(NLMT,J)

IF(NDBC(JG)0NE°2) GO TO 510

DDX w X(I) = X(J)

DDY

Y(I) m Y(J)

DD = SQRT( DDX**2 + DDY**2 )

R(IG) : R(IG) + HH * TA * DD / 200

1

+ THETAM *HH*DD * (TOLD(IG)/3.0 + TOLD(JG)/6.0)

A(IG,NDG) = A(IG,NDG) + THETA*HH*DD/3.0

JGBD = JG - IG + NDG

A(IG,JGBD) = A(IG,JGBD) + IHETA*HH+DD / 6.0

510

CONTINUE

599

CONTINUE

IJ = I + 1 IF(IJ°EQ.4) IJ IK = IJ + 1 IF(IK0EQ°4) IK II H ll H

m O O O O O Q O 0 0 0 0 äO O O O O O O h O O O K O _U1 500 600 1000 7070 7701 1010 0 0 0 0 0 1 2 3 DO 500 J = 1,3 JI = J - 1 31/ 15 IF(JI.EQOO) JI = 3 JK = J + 1 IF(JK.EQ.4) JK = 1 JG = NODEL(NLMT,J) JGBD = JG ' IG + NDG

This term ( DRV ) is part of the element quadrature *

*

DRV (Y(IJ)-Y(IK)) +

(X(IK)-X(IJ)) BIJ = BCE

IF( IOEQDJ) BIJ = 2. IF(KFRZ.NE.1) GO TO

(Y(JK)-Y(JI)) (X(JI)-X(JK)) O * BCF

480 For the phase change element

BIJ = BIJ .

+ ELDSDT*TAU*(2.0*(ELB(I)-ELA(I))*(ELB(J)-ELA(J)) + 3.0 * ELB(I)'* ELA(J)

+ 300 * ELB(J) * ELA(I) ) / 6.0

The next term distributes the appropriate conductivity over

the correct proportions in the phase change element.

CIJ = VK2 * DRV / ( 460 * AE)

+ (VKl-VK2)*DRV*A1 / ( 4.0 * AE * AE) GO TO 495

For the other elements

CIJ = C4A * DRV

Load the matrix A and the vector R for the system.solver

The THETA parameter controls numerical implicitness

A(IG;JGBD) = A(IG,JGBD) + BIJ + THETA * CIJ R(IG) = R(IG) + ( BIJ + THETAM * CIJ ) * TOLD(JG) CONTINUE

CONTINUE CONTINUE

IF(KTPLOT.NE.IZDPLT .OR. ITER.NE.O) GO TO 1010 WRITE(IO,7070) TIMOL

FORMAT(//,12X,'Plotting data stored for TIME =',1P1E12.3)

KTPLOT = 0

VRYLGE = 1.1E30

VITM = IT " 100

WRITE(33,7701) VRYLGE, VITM, TIMOL, TIMOL

FORMAT(1P4E11.2)

CONTINUE

Determine "ZERO", i.e. a magnitude below which significance

is lost. IF( IT*(ITER+1) AMAX 0.0 DO 1100 I = 1,NBND DO 1100 J = 1,NBW ABSVAL = ABS(A(I,J))

IF(ABSVAL.GT.AMAX) AMAX = ABSVAL .NE. 1 ) GO TO 1101

1100 CONTINUE

ZERO = AMAX * 000001

31/ 16

C C 1101 CONTINUE RETURN E N D C C C

C The following routine solves the system AX = B C

C

SUBROUTINE BANSAL (A,NBW,NEQ,NMAX,B,X,INEW,KERR,IO)

C

C

COMMON /SMALL/ ZERO

C

DIMENSION A(NMAX,NBW), B(NMAX), X(NMAX)

C C * C Initialize counters C KERR = 0 NBWT = NBW N = NEQ NBZ = NBWT / 2 + 1 N31 3 NBZ m 1 Ni = N m 1 X(I) 10 CONTINUE IF(INEW.NE,O) GO TO 350

Forward elimination of matrix

01 0 0 0 0 Kl = NBl IMl = 1 DO 300 I = 2,N Kl = MINO( Kl+1,N) K = NBl AII = A(IM1,NBZ) IF( (AIIOLTOZERO).AND.(AII.GT."ZERO)) GO TO 1850 D0 200 J = I,K1 CX A(J,K) / AII Ll NB2 + 1 KPB K + NBl KPl K + 1 DO 100 L = KP1,KPB

A(J,L) = A(J,L) = CX * A(IM1,L1)

100 Ll = Ll + 1 C A(J,K) = CX X(J) = X(J) m X(IM1) * CX _ K = K 1 200 CONTINUE IMl = I 300 CONTINUE GO TO 550 C

C Operate on the new RHS a is done after the first pass C

350 IN = NBl IMl 2 1

DO 500 I = 2,N

600 700 1850 C 1 Jl = NBl DO ÄOO J = 1,IN X(J) = X(J) m X(IMl) * A(J,K1) Kl = Kl m 1 CONTINUE IMl = I

Back ward substitution

X(N) = X(N) / A(N,NB2)

DO 700 I = 1,N1

SUM.= 0.0

L : N32 + 1

K = N m I

Kl = MINO(N1, K + NBlMl)

DO 600 J = K,Kl

SUM = SUM.+ A(K,L) * X(J+1)

L = L + 1

_

CONTINUE

X(K) = ( X(K) - SUM) / A(K,NBZ)

CONTINUE

'

RETURN

WRITE(IO,9000) IMl

9000 FORMAT(1HO,20(4H*!?*)/10X,'Error from BANSAL -',IS, E

/'Diag0nal element reduced to zero')

N D

24 1 H U ' I O k O O O Q m U ' I 0 5 wN H 28 0.5083 0.7559 7.2E-3 4.2Ewå 36.0 0.5083 0.7559 7.2E-3 4.2E-3 36.0 10 2 H o m m q m m øwN ø-â H H P J F H A a wa r e h wd o wn m o n a h ua r 4 h ác n o c n -J N N N . m m m 6 28 4O XYFREZ t vb o h a un vb o h a un vb o h o h ua k a k a h wd r d h wa + 4 k a h ua r a k a h ua â -Ål -l H O O O O O O O O O , O O O O O O O O O H H H H o m m m q m m m m øwwww t a r d k áh ud øwøb o uo k a k n a r A h J N J O N U H J h a h a h ud t d F H D K Q O J G D 21 xl k D H M G Q M M L D N wøb N N W H H N H H l -' Q K D H 2

Sample Problem: Freezing of a Wedge..

°OOOE+02 .239E+Ol .071E+Ol 0500E+02 0386E+02 °061E+02 OOOOE+02 0848E+02 ,414E+02 .OOOE+02 .772E+02 .121E+02 .000E+02 .696E+02 .828E+02 .500E+02 °O81E+02 .889E+02 .500E+02 .929E+02 .303E+02 9000E+03 .239E+02 .O71E+02 ( D xl o m U ' l a h -J U ' I N H o b 1.0E9 0.65 0.65 G 4.0 2 0.2500E+O7 6 2 l 20501E+OO 2 20501E+00 Q W O M N O W N O N H O N H O H Q O H M O Q W O.OOOEEOl .827E+01 0071E+01 DOOOEa01 074OE+01 0061E+02 OOOOEEOI 0654E+01 e414E+02 .OOOEWOI .148E+02 0121E+02 QOOOEDOI 0531E+02 0828E+02 OOOOEGOI 0105E+02 .889E+02 .OOOE=01 0870E+02 0303E+02 °OOOE==01 .827E+02 .O71E+02 =0.1E=15 -O.1E-15

81/ 18

3 29501E+00 22 4°OOOE+OO 23 4.000E+OO 24 4.000E+00 0.1000E+08 6 19060E+00 1.06OE+00 1906OE+00 3.998E+00 3.998E+00 24 3.998E+OO O.2250E+08 6 1 1.605E-01 2 19605E-01 3 1.605E-01 22 30944E+OO 23 3.944E+00 24 3.944E+00 N N ( J O N U U N H O.4000E+08 6 2 l -4.7S4E-01 2 -40754E*01 3 -40754E°01 22 30796E+00 23 35796E+00 24 30796E+00 O.625E+08 6 2 l -9.669E°01 2 -99669EUOI 3 -9.669E-01 22 30588E+00 23 39588E+00 24 3.588E+00 0.9000E+08 6 1 -1.369E+00 2 -1.369E+00 3 -10369E+00 22 30359E+OO 23 3°359E+00 24 3.359E+00 O.1225E+09 6 l -1.709E+OO 2 -16709E+00 3 -l.709E+00 22 3.130E+00 23 3013OE+00 24 30130E+00 O.1600E+09 6 1 '2.004E+OO 2 -2.004E+OO 3 -2.004E+OO 22 20908E+00 23 20908E+OO 24 2.908E+OO O.2025E+O9 6 1 -2.264E+00 2 -2.264E+OO 3 -2.264E+OO 22 2.698E+00 23 2.698E+OO 24 20698E+OO O.2500E+09 6 1 -2.497E+OO 2 -2.497E+00 3 -2.497E+OO 22 2.501E+OO 23 2.501E+00 24 2.501E+OO 0.3025E+09 6 Bl/ 19

1 2,707E+OO 2 -2.707E+00 3 =2.707E+00 22 20316E+00 23 20316E+OO 24 2.316E+00 O.3600E+O9 6 1 20899E+00 2 20899E+00 3 20899E+OO 22 2°142E+00 23 2°l42E+OO 24 29142E+00 0.4225E+09 6 1 30076E+00 2 =30076E+00 3 °30076E+OO 22 19979E+OO 23 16979E+00 24 10979E+00 004900E+09 6 l =30240E+00 2 ©30240E+00 3 30240E+OO

22

1.825E+00

23

10825E+00

24

19825E+OO

O.5625E+O9 6

1 39393E+00 2 3.393E+00 3 30393E+OO 22 10680E+00 23 10680E+00 24 10680E+00 0064OOE+09 6 1 ==30535E+00 2 ©30535E+00 3 30535E+OO 22 10543E+00 23 1.543E+00 24 10543E+00 007225E+09 6 1 =3°669E+00 l =3°669E+OO 1 =3°669E+00 22 10413E+00 23 1,413E+OO 24 10413E+OO Oo81OOE+O9 6 1 30795E+00 2 3.795E+00 3 3.795E+OO 22 1029OE+00 23 10290E+00 24 1029OE+OO 009025E+09 6 1 30915E+00 2 m3°915E+OO 3 39915E+OO 22 10172E+00 23 19172E+OO 24 10172E+OO 0°1000E+1O 6 l ©4.028E+OO 2 G40028E+00 3 40028E+00 22 10060E+00 23 1°O6OE+OO

81/ 20

24

_1.06QE+00

XYFREZ Sample Problem: Freezing of a Wedge. NBND 24 Element NMBCl 6 IZDPLT 5 BC l 1 1 0 0 0 0 0 O 0 O 0 0 0 0 0 O 0 0 0 O 1 1 1 Material 1 1 1 1 1 1 1 l 1 2 2 2 2 2 Q K D H U I O N Q W W W N M Q N N W H H N H H H Q K D H Hü-åo l-*D -* I-*H l-I m m m m wm o xøq m p ww NELS 28 MSHPRT X °OOE+02 024E+01 007E+01 050E+02 039E+02 006E+02 .OOE+02 .85E+02 .41E+02 900E+02 °77E+02 012E+02 OOOE+02 970E+02 083E+02 050E+02 008E+02 .89E+02 °50E+02 093E+02 030E+02 OOOE+03 024E+02 007E+02 Nodes 1 KPAUSE = 1 d M O M N O M N O N H O N H O H Q O B -* W O Q M O NTS 40 KBVPRT Y .OOE+OO 083E+01 007E+01 000E+00 074E+01 006E+02 °OOE+00 065E+01 °4lE+02 500E+OO 015E+02 012E+02 OOOE+00 053E+02 .83E+02 .00E+00 .11E+02 089E+02 OOOE+00 087E+02 030E+02 OOOE+00 .83E+02 .O7E+02 ITRLM 2 Element

Thermal Properties of Material l

0 Mate t h J h 3 k J A D N H A §4 1 4 k 4 k 4 k 4 k 4 k 4 KPRTVL 10 rial R O C D O X W O O N 12 14 15 17 18 20 21

81/ 22

k o m m m 7 11 12 14 15 17 18 20 23 24

HETF 5.08E-01 HETF 5.08E-01 1 4 4 4 7 4 10 4 13 4 16 4 19 4 22 4 TIME 1 3 22 4 TIME 1 2 22 4 TIME 1 1 22 3 HETU CNDF CNDU EL TF

7.56E"01 7.20E°03 4.20E-O3 3.60E+01 -1.00E-16

Thermal Properties of Material 2

HETU CNDF CNDU EL TF

7.56E-01 7.20E-03 4.20E-03 3.60E+Ol -1.00E-16

Solution Control Parameters

TIMLIM

1.000E+09 6.500E-01AMP

Initial Temperatures at Each Node o000E+00 2 4.000E+00 0000E+OO 5 4.000E+OO 0000E+00 8 4.000E+00 0000E+00 ll 4.000E+OO GOOOE+OO 14 4.000E+OO GOOOE+OO 17 40000E+00 OOOOE+OO 20 4.000E+00 .000E+00 23 4.000E+00

Time and DT at Each Time Step, with

Boundary Values at Each Node

STEP = 1 TIME = 1.250E+06

Nodes and their temperatures

0250E+00 2 3.250E+00 .OOOE+OO 23 4.000E+00

STEP = 2 TIME = 2.500E+06 Nodes and their temperatures .501E+OO 2 2.501E+On .000E+OO 23 4.000E+OO

STEP = 3 TIME = 6.250E+06

Nodes and their temperatures

.780E+OO 2 1.780E+OO .999E+00 23 3.999E+OO THETA 6.500E-01 3 12 15 18 21 24 DT DT

DT,

24 k øb k t âo ân âøb h däc k) .000E+OO .000E+OO .OOOE+OO .OOOE+OO .000E+00 .000E+00 .000E+00 .000E+00 1.250E+06 .250E+00 .OOOE+OO 1.250E+06 .501E+00 .OOOE+OO 3.750E+06 .780E+00 .999E+OO

31/ 231

m m ( M I -l STEP = 4 Nodes ,060E+00 9998E+00 STEP = 5 Nodes .102E-Ol .971E+00 STEP = 6 Nodes .605E-01 .944E+00 STEP = 7 Nodes .575E°01 .87OE+OO STEP = 8 Nodes 0754E°01 .796E+OO STEP = 9 Nodes 0212E=Ol .692E+00 STEP = 10 Nodes TIME = 1.000E+O7 their temperatures 2 1.060E+00 23 39998E+00 TIME = 10625E+07 their temperatures . 2 69102E 01 23 36971E+00 TIME = 20250E+07 their temperatures 2 10605Ew01 23 30944E+00 TIME = 3°125E+07 their temperatures 2 -10575E-01 23 30870E+00 TIME = 4.000E+O7 their temperatures 2 =4.754E-01 23 3.796E+00

TIME = .5.1255+07

their temperatures 2 76212E-01 23 30692E+00 TIME : 6.250E+07 their temperatures wo n wH 3.750E+06 .060E+00 .998E+00 6.250E+06 .102E-01 .971E+00 6.250E+06 .605E-Ol .944E+00 8.750E+06 =1.575E-Ol .870E+00 8.750E+06 .754E-01 .796E+OO 10125E+O7 -7.212E-Ol .692E+00 1.125E+07

Bl/ 24

1 =9°669E-Ol 2 -9.669E-01 3 -9.669E-01 31/ 25 :22 39588E+00 23 3.588E+OO 24 3.588E+OO

TIME STEP = 11 TIME = 7.625E+O7 DT = 1.375E+07

Nodes and their temperatures

1 19168E+00 2 -1.168E+00 3 -10168E+00 22 3.474E+00 23 3.474E+OO 24 3.474E+OO

TIME STEP = 12 TIME = 9.000E+07 DT = 1.375E+O7

Nedes and their temperatures

1 u1.369E+00 ' 2 -1.369E+00 3 »1.369E+00 22 3.359E+00 23 30359E+00 24 3.359E+00

TIME STEP = 13 TIME = 1.063E+O8 DT = 1.625E+07

Nodes and their temperatures

1 -l.539E+00 2 1.539E+00 3 -l.539E+00 22 30245E+00 23 3.245E+00 24 3.245E+00

TIME STEP = 14 TIME = 1.225E+O8 DT = 1.625E+O7

Nodes and their temperatures

1 -1.709E+OO 2 -1.709E+OO 3 -1.709E+00 22 3.130E+00 23 3.130E+OO 24 3.13OE+OO

TIME STEP = 15 TIME = 1.413E+08 DT = 1.875E+07

Nodes and their temperatures

l -1.857E+00 2 -1.857E+00 3 -l.857E+00 22 3.019E+00 23 3.019E+OO 24 3.019E+OO

TIME STEP = 16 TIME = 1.6OOE+O8 DT = 1.875E+07 Nodes and their temperatures .

1 -2.004E+00 2 -2.004E+OO 3 -2.004E+OO 22 2.908E+00 23 2°908E+00 24 2.908E+OO

Nodes 0134E+00 0803E+OO STEP = 18 Nodes 9264E+OO 0698E+OO STEP = 19 Nodes 0381E+00 9599E+00 STEP = 20 Nodes 0497E+00 0501E+OO STEP = 21 Nodes 0602E+00 0408E+00 STEP = 22 ches 0707E+OO 0316E+00 STEP = 23 o Nodes 6803E+00 0229E+OO their temperatures 2 2.134E+00 23 2.803E+00 TIME = 2.025E+08 their temperatures 2 w2°264E+00 23 2.698E+OO TIME = 2.263E+08 their temperatures 2 020381E+00 23 20599E+00 TIME = 2.500E+08 their temperatures 2 =2°497E+00 23 20501E+00 TIME 2 2.763E+O8 their temperatures 2 m2°602E+00 23 29408E+00 TIME 2 3.025E+08 their temperatures 2 m2°707E+OO 23 20316E+00 TIME = 3.313E+O8 their temperatures 2 -2.803E+OO 23 2.229E+00 0134E+00 0803E+OO 2.125E+07 0264E+OO 0698E+00 2.375E+07 .381E+00 0599E+00 2.375E+O7 0497E+OO 0501E+OO 2.625E+O7 .602E+00 0408E+00 2.625E+07 .707E+00 .316E+00 2.875E+O7 .803E+00 .229E+OO

81/ 26

TIME

1 2

22

2

TIME

1 sz

22 2

TIME

1 3

22 1

TIME

1 ==3

22 1.,

TIME

1 »3

2.2 1

TIME

1 -3

22 1

TIME

1 »3

22 1

STEP = 24 Nodes D899E+00 9142E+OO STEP = 25 Nodes 0987E+00 6060E+00 STEP = 26 Nodes .076E+00 9979E+00 STEP = 27 Nodes 0158E+00 902E+00 STEP = 28 Nodes 0240E+00 .825E+00 STEP = 29 Nodes .316E+00 .753E+OO STEP = 30 Nodes .393E+OO §680E+OO and and TIME = and and and and and TIME = 3.6OOE+08 their temperatures 2 -2.899E+00 23 2.142E+OO TIME = 3.913E+08 their temperatures 2 -2.987E+OO 23 2.060E+OO 4.225E+08 their temperatures 2 3.076E+00 23 1.979E+00 TIME = 4.563E+08 their temperatures 2 -3.158E+OO 23 1.902E+00 TIME = 4.900E+08 their temperatures 2 »3.240E+00 23 10825E+00 TIME = 5.263E+O8 their temperatures 2 -3.316E+OO 23 1.753E+00 TIME = 5.625E+08 their temperatures 2 -3.393E+00 23 1.680E+00 DT DT DT b b . ) h ) 2.875E+O7 .899E+OO .142E+00 3.125E+O7 .987E+00 .060E+00 3.125E+07 .076E+00 °979E+OO 3.375E+07 .158E+OO .902E+00 3.375E+O7 5240E+00 .825E+00 3.625E+07 .316E+OO .753E+00 3.625E+O7 .393E+OO .680E+OO

81/ 27

TIME 1 =3 22 l TIME 1 3 22 l TIME 1 3 22 1 TIME 1 3 22 1. TIME 1 3 22 1 TIME 1 -3 22 1 TIME 1 "3 22 1 STEP = 31 Nodes .464E+00 .612E+00 STEP = 32 Nodes °535E+00 .543E+OO STEP = 33 Nodes .602E+00 .478E+00 STEP = 34 Nodes .669E+OO 413E+OO STEP = 35 Nodes §732E+00 .352E+00 STEP = 36 Nodes 0795E+OO .29OE+OO STEP = ches .855E+00 .231E+OO and and and and and and 37' and TIME = 6.013E+08 their temperatures 2 =3.464E+00 23 1.612E+00 TIME = 604OOE+08 their temperatures 2 »3.535E+00 23 10543E+00 TIME = 6.813E+08 their temperatures 2 09000E+00 23 1°478E+00 TIME = 70225E+08 their temperatures 2 OGOOOE+OO 23 1.413E+00 TIME = 7.663E+08 their temperatures 2 -10898E+00 23 10352E+OO TIME = 8.100E+08 their temperatures -39795E+ " 23 1029OE+UU TIME = 8.563E+08 their temperatures 2 -39855E+00 23 10231E+00 DT DT DT H O H O I H ut ) 3.875E+07 .464E+00 °612E+00 30875E+O7 0535E+OO 0543E+00 4.125E+07 0000E+00 0478E+OO 49125E+O7 .OOOE+OO .413E+00 4.375E+07 0898E+00 0352E+OO 4.375E+O7 ."95E+00 .290E+00 4.625E+07 .855E+00 .231E+OO 81/ 28

TIME STEP = 38 TIME = 9.025E+08 DT =

Nodes and their temperatures

1 -3.915E+00 2 -3.915E+OO 3 -3

22 1.172E+00 23 1.172E+00 24 1

TIME STEP = 39 TIME = 9.513E+08 DT =

Nodes and their temperatures

l 03.971E+00 _ 2 -3.971E+00 3 -3

22 1.116E+OO 23 l.1l6E+OO 24 1

TIME STEP = 40 TIME = 1.000E+09 DT =

Nodes and their temperatures

1 w4°028E+00 2 -4.028E+00 3 -4

22 1°060E+00 23 1.060E+00 24 1

Calculated Temperature Solution

Plotting data stored for TIME = 1.625E+07 TIME = 6.250E+O7 TIME STEP =

81/ 29

4.625E+O7 .915E+OO .172E+00 4.875E+O7 .971E+OO .116E+OO 4.875E+O7 .028E+00 .060E+OO 10

NODE 10 13 16 19 22 NODE 10 13 16 19 22 NODE 10 13 TEMP °669E==01 0442Eø02 .637Eø01 0450E+OO 9058E+00 .681E+00 .207E+00 .588E+00 l l l wwN N l -* U ' I O O KD NODE 2 5 8 11 14 17 20 23 TEMP 9,669E001 .568Em02 9616E-01 .447E+00 .053E+00 .674E+00 .200E+00 .588E+00 S M W N N l -U ' U O D

Plotting data stored før TIME =

Plotting data stored for TIME =

TEMP 2.497E+00 -1.612E+00 9.806E=01 =10017Em01 50515E°01 1.259E+OO 1.924E+00 2.501E+OO TIME = NODE 2.500E+08 TEMP 2.497E+00 510612E+00 9.813E-Ol ©16025E-01 50500Em01 1.257E+00 1.921E+00 2.501E+00

Plotting data stored for TIME =

Plotting data støred før TIME

TEMP 3.393E+00 =2.498E+00 =1.859E+00 9.658Em01 °3.284E°01 3.907E°01 1.068E+00 1.680E+00 TIME = NODE 2 5 8 11 14 17 20 23 5.625E+08 TEMP 3.393E+00 "2.498E+00 *1.859E+OO 9.658E-01 3.289E-01 3.894E-01 1.067E+00 1.680E+00

Plotting data støred for TIME =

Plotting data stored for TIME

TEMP Q4.028E+00 -39102E+00 26441E+00 10517E+00 89564E-01 TIME = NODE 2 5 8 11 14 1.000E+O9 TEMP 40028E+00 "30102E+00 =2°441E+00 =10517E+00 "89578E-01 NODE 3 6 9 12 15 18 21 24 6.250E+07 1.413E+08 2.500E+08 3.913E+08 NODE 3 6 9 12 15 18 21 24 5.625E+08 e f. 73:24:15.? TIME STEP NODE 3 6 9 12 15 9 I D J U O N N H U ' I O O TIME STEP = TIME STEP = TEMP 669Ew01 365E-02 634Ew01 450E+00 058E+00 681E+00 207E+OO 588E+00 20; , TEME .497E+00 °611E+00 .809E=01 °021E=01 .511E-01 .259E+00 .924E+00 .501E+00 30 TEMP .393E+00 .497E+00 .86OE+00 .662E-01 .287E-01 .906E-01 .O68E+00 °680E+00 40 TEMP .028E+OO 9101E+00 9441E+00 0517E+00 .564E-01

B1/ 30

16 -1.271E°Ol 17 1.339E-01 18 -l.266E-01 i 81/ 31 19 5.015E-Ol 20 50118E-Ol 21 5.016E-01 5

22 1,06OE+OO 23 1.060E+OO 24 1.060E+OO

O O O O O O O O O O O O Q O

Bi1aga 2

32/ 1

--- Program FEHEAT --- -- begin --- nu

PROGRAM FEHEAT

This program solves 2 dimensional steady state heat trans fer problems using the finite element method. (Only trianm

gular element may be used).

Dimension the matrices which will not be stored in common storage.

DIMENSION NMJ(300),NEC(300) DIMENSION BT(175)

CHARACTER*6 IDOC1,IDOC2,IDOC3,IDOC4,IDOC5,IDOC6,IDOC7,IDOC8, IDOC9,IDOCll,IDOClZ,IDOClB,IDOCl4,IDOC15

Dimension the remaining matrices into common block storage COMMON/IDY/IDOCI,IDOCZ,IDOC3,IDOC4,IDOC5,IDOC6,IDOC7,IDOC8, IDOC9,IDOC11,IDOClZ,IDOCl3,IDOC14,IDOC15 COMMON/MlI/IMAT(300) ' COMMON/MlR/TKM(300),QI(300) COMMON/M12I/NC1(300),NC2(300) COMMON/M12R/X(l75),Y(175),H(300),FQIE(175) COMMON/M2I/MC1,MC2,NN,NCASE,NBHF(175) COMMON/M2R/BHF(175),TA(300),RHST(175),RHS(300) COMMON/M23I/MC,NBT(175) COMMON/M23R/R1(175) COMMON/M3I/NE,NCN,NDF,NSZF,NBAND COMMON/M3R/TM(30,175) COMMON/M13I/NODE(300,3) C0MM0N/M13R/TME(3,3 DATA IDOC3/'NAMNET'/ OPEN(UNIT=3,FILE=IDOC3,STATUS='OLD') READ(3,3010)IDOC5 READ(3,3010)IDOC6 READ(3,3010)IDOC7 READ(3;3010)IDOC8 READ(3,3010)IDOC9 READ(3,3010)IDOC11 READ(3,3010)IDOC12 READ(3;3010)IDOC13 READ(3,3010)IDOC14 FORMAT (A) OPEN(UNIT=ll,FILE=IDOCll,STATUS='OLD') READ(ll,107,END=1005)NN,NE,MC,MJC,MC1,MC2,NIT CONTINUE

Initialise all matrices and parameters to zero.

NCN=3 NDF=1 NCASE=1 NSZF=NN * NDF TKM(1)=O.833 D0 2 I=1,NN,1 X(I)=OoO Y(I)=O°O NBT(I BT(I

=0.0

)=0

NBHF(I)=

O

O O O O N m 1000 O 0 0 0 4 l 2 BHF(I)=0.0 RHS(I)=O°O FQIE(I)=O°O CONTINUE

82/ 2

DO 28 J=1,NE,1 NC1(J)=O NC2(J)=O H(J)=0.0 TA(J)=0.0 IMAT(J)=O NMJ(J)=O QI(J):O°O CONTINUE

Read in the nodal point data. OPEN(UNIT=5,FILE=IDOC5,STATUS='OLD')

D0 3 J=1,NN,l

READ(5,101,END=1000) I,X(I),Y(I) CONTINUE

CONTINUE

Read in the element data.

OPEN(UNIT=6,EILE=IDOC6,STATUS='OLD')

OPEN(UNIT=l7,FILE='GEOMETRY',STATUS='NEW')

D0 4 I=l,NE,l

READ(6,102,END=1001) J,NODE(J,l),NODE(J,2),NODE(J,3),IMAT(J)

WRITE(17,2222)J,NODE(J,l),X(NODE(J,1)),Y(NODE(J,1)),

NODE(J,2),X(NODE(J,2)),Y(NODE(J,2)),

NODE(J,3),X(NODE(J,3)),Y(NODE(J,3)),IMAT(J) 2222 FORMAT(2I5,2F9.4;IS,2F9.4,I5,2F9.4,15) CONTINUE CONTINUE 1001 0 O 0 0 0 0 O N O O H O O O O S Q O O U B O O O O H O N 0 x]

Read in the fixed temperature bcundary 9conditions. IF (MCQEQOO) GO TO 5

OPEN(UNIT=7,FILE=IDOC7,STATUS='OLD')

D0 6 IM:1,MC,1

READ(7,103,END=1002) M,NBT(M),BT(NBT(M)) CONTINUE CONTINUE

Read in the internal heat generation values for each element. IF(MJC0EQ.O) GO TO l OPEN(UNIT=8,FILE=IDOC8,STATUS='OLD') D0 7 IMJ=1,MJC,1 CONTINUE CONTINUE

Read in the heat flux value for nodes where heat flux is specified.

IF (MC1.EQ.O) GO TO 8

OPEN(UNIT:9,FILE3IDOC9,STATUS='OLD') D0 9 IMC=1,MC1,1

O 0 (D 0 O m O OO O H k D 82 91 83 93 92 84 94 97 85 95 99 86 96 READ(9,103,END=1008) M1,NBHF(M1),BHF(NBHF(M1)) CONTINUE CONTINUE

82/ 3

Fer boundary segments subject to convective heat transfer

read in corresponding nodes ,H, and ambient temperature.

IF (MC2.EQ00) GO T0 10 OPEN(UNIT=14,FILE=IDOC14,STATUS='OLD')

DO 11 M2=1,Mc2,1

READ(14,106,END=1009) NEC(M2),NC1(NEC(M2)),NC2(NEC(M2)),

H(NEC(M2)), TA(NEC(M2))

CONTINUE

CONTINUE

CONTINUE

Print down the control text.

OPEN(UNIT=12,FILE=IDOC12,STATUS='NEW') WRITE(12,120) NN,NE WRITE(12,121)

DO 81 I=1,NN,1 '

WRITE(12,101) I,X(I),Y(I)

CONTINUE

WRITE(12,122) DO 82 J=1,NE,1 WRITE(12,102) J,NODE(J,1),NODE(J,2),NODE(J,3),IMAT(J) CONTINUE IF (MC) 93,93,91 WRITE(12,123) DO 83 M=1,MC,l WRITE(12,103) M,NBT(M),BT(NBT(M)) CONTINUE CONTINUE IF(MJC) 94,94,92 WRITE(12,129) DO 84 I=l,MJC,l WRITE(12,103) I,NMJ(I),QI(NMJ(I)) CONTINUE CONTINUE IF (MCl) 95,95,97 WRITE(12,125) DO 85 Ml=l,MCl,l WRITE(12,103) M1,NBHF(M1),BHF(NBHF(M1)) CONTINUE CONTINUE IF (MCZ) 96,96,99 WRITE(12,126) DO 86 M2=1,MC2,1 WRITE(12,126) NEC(M2),NC1(NEC(M2)),NC2(NEC(M2)), H(NEC(M2)), TA(NEC(M2)) CONTINUE CONTINUE

MX=0

2 4

DO 66 I=1,NE,1 B /

L1=IABS(NODE(I,1)°NODE(I,3))

L2=IABS(NODE(I,3)GNODE(I,2))

L3=IABs(NODE(I,2)-NODE(I,1 )

MX=MAXO(MX,L1,L2,L3)

66

CONTINUE

NBAND = NDF * (MX+1)

WRITE(12,109)NBAND

MUD=NBAND-1

IO=O IOP=1 C CALL FORMK CALL FRI-IS (BT) CALL MCHB(RHS,TM,NSZF,1,MUD,IOP,O.5E-6,IER,IO) C WRITE(12,127) C DO 200 I=l,NN,l WRITE(12,130) I,RHS(I) 200 CONTINUE C IF (NIT.EQ.O) GO TO 228 CALL ISQ§HM(RHS,NODE,X,YyNE,NIT) 228 CONTINUE : CLOSE(UNIT=12) 101 FORMAT(IS,2F10.4) 102 FORMAT(I5;4IÖ) 103 FORMAT(I5,IG,F10.4) 106 FORMAT(316,2F1004) 107 FORMAT(7IG)

109 FORMAT(/,SX,'THE BANDWIDTH IS ',I6)

120 FORMAT(//,5X,'TOTAL NUMBER OF NODES =',I6,3X,

*

'TOTAL NUMBER OF ELEMENTS = ',I6)

121

FORMAT(//,2x,'N0DE',3x,'GLOBAL',4X,'GLOBAL',/,3x,'No.'

,6x,

*

'X',9x,'x',1xy/,i mmmmmm -u--- --'

C

122

FORMAT(//,1X, ELEMENT ,1X,'NODE',2X,'NODE',2X,'NODE ,2

X,

* 'MATLMy ,3xiiNooi,4x, 1 ,5x,I2',5X,'3',4x,'TYPE',/,1x,

* I m a m _ m n u, _ _ _ m n a c., _ _ u m m m D _ m a a _ gaf)

C

123

FORMAT(//,3X,'B.Co',2X,'NODE',2X,'TEMP',/,4X,'NO.',3X,'NO.',2X,

* I m - u u - nm m n u s _ m u _

-C

129 FORMAT(//;2X,'B.C°',2X,'NODE',6X,'QI' ,/,2X,'NO.',3X,'NO.',5X, * '(UNITS)')

125

FORMAT(//,2X,'B°C°',2X,'NODE',6X,'HF'

,/,2X,'NO.',3X,'NO.',5X,

* '(UNITS)')

C 126 FORMAT(//,1X, SEGMENT',1X,'NODE',2X,'NODE',4X,'H',3X,'AMB',/,3X, * 'NO.',4X,'1',5X,'2',3X,'(UNITS)',2X,'TEMP(F)'> C 127 FORMAT(//,13X;'NODE',2X,'TEMPERATURE (F)', * /,12X,' --- ma ---') C 130 FORMAT(10X,I6,F13.2) C STOP END C

C °°°°°°° um FEHEAT program --- -- end --- -a

O O O O O O O O O O 300 O O ' O O O 305 O O O 310 320 330 350 400 500 169 166 167

This subroutine is from FEHEAT program.

SUBROUTINE FORMK

Forms stiffness matrix in upper triangular form. COMMON/M23I/MC,NBT(175) COMMON/M23R/Rl(l75) COMMON/MBI/NE,NCN,NDF,NSZF,NBAND COMMON/M3R/TM(30,175) COMMON/M13I/NODE(300,3) COMMON/M13R/TME(3,3) DIMENSION ST(5250) EQUIVALENCE (ST(1),TM(1,1)) Zerc stiffness matrix.

DOu300 N=1,NSZF DO 300 M?1,NBAND TM(M,N)=0 Scan elements. DO 400 N=1,NE CALL TSM(N)

Return ESTIFM as stiffness matrix.

StOre ESTIFM in SK. First rcws. DO 350 JJ=1,NCN NROWB=(NODE(N,JJ)-l)*NDF DO 350 J=1,NDF IF (NROWB) 350,305,305 NROWB=NROWB+1 I=(JJ-1)*NDE+J Then columns. DO 330 KK=1,NCN NCOLB=(NODE(N,KK)°1) * NDF DO 320 K=1,NDF L=(KK-l) * NDF + K NCOL=NCOLB + K + 1 - NROWB

Skip storing if below band.

IF(NCOL) 320,320,310

TM(NCOL,NROWB) = TM(NCOL,NROWB) + TME(I,L) CONTINUE

CONTINUE CONTINUE CONTINUE

Insert boundary conditions.

DO 500 N=1,MC

I=NBT(N)

TM(1,I) = TM(1,I) * l.E15

R1(I)=TM(1,I)

CONTINUE

D0 1 J=1,NSZF

FORMAT('ROW',I2)

FORMAT(1OE10.2)

FORMAT(//)

BZ/S