Большая База Рефератов - Решение обратной задачи вихретокового контроля - бесплатно рефераты, скачать рефераты, рефераты на тему

Рефераты. Решение обратной задачи вихретокового контроля

Аппроксимация гиперболическим тангенсом (HTG)

В случае задания аппроксимации гиперболическим тангенсом зависимость электропрводности от толщины представляется в виде

SIGMA = si2 + ( si1-si2 )/2*{ 1 + th( ( beta-x )/gamma ) }

Величина и ограничения для параметров si2,beta,gamma задаются начиная со второй строки в "special approximation parameters", для si1 аналогично siI.

П1.3 Результаты расчета

Результаты расчета помещаются в текстовый файл( кодировка ASCII, расширение по умолчанию LST ), при этом результат каждой итерации отбражается строкой вида:

1 <Ф> 0.000353 Rg= 17.003 15.639 9.697

где первая цифра (в данном случае 1) соответствует номеру текущей внутренней итерации, затем после текста "<Ф>" идет значение текущей абсолютной среднеквадратичной невязки по всем гармоникам (в данном случае 0.000353), затем после текста "Rg= ", идут искомые текущие значения переменных минимизации. В случае SPL,PWL,PWC аппроксимаций это непосредственно узловые значения электропроводности для текущего количества узлов, а для EXP,HTG аппроксимаций это параметры { siE, siI, Alfa } или { si1, si2, Beta, Gamma}. B качестве последней строки помещаются nPoints вычисленных значений э/проводности в равномерно расположенных узлах пластины.

П1.4 Основная программа ErIn

(****************************************************************************)

(* ErIn v1.42 *)

(* Eddy current inverse problem solver. *)

(* Moscow Power Engineering Institute , Introscopy dept. *)

{****************************************************************************}

Program ErIn;{23.02.99}

Uses

DOS,CRT, EData, EMath, EDirect, EFile, EMinimum;

Var

m, mLast, i : byte; {loop counters}

procedure about; {Let me to introduce myself}

begin

clrscr;

GetTime( clk1.H, clk1.M, clk1.S, clk1.S100 ); {get start time}

writeln('***********************************************************');

writeln('* ErIn v1.42 Basic * *');

writeln('***********************************************************');

end;

procedure initParameters;

var

apDT : byte; {approximation type for direct task}

begin

apDT := nApprox SHR 4; {XXXXYYYY->0000XXXX}

fHypTg:=(( apDT AND apHypTg ) = apHypTg);

if fHypTg then

begin

si0[ 1 ]:=si[ 1 ]; {si1 - conductivity about bottom of slab}

si0[ 2 ]:=par0[ 2 ]; {si2 - conductivity about top of slab}

si0[ 3 ]:=par0[ 3 ]; {Beta - ratio of approx.}

si0[ 4 ]:=par0[ 4 ]; {Gamma- ratio of approx.}

mCur:=4;

end

else

if(( apDT AND apExp ) = 0 ) then {It's not an EXP approx.}

begin

for i:=1 to nPoints do si0[ i ] :=si [ i ]; {SI data from file}

mCur:=nPoints;

end

else

begin

si0[ 1 ]:=si[ 1 ]; {siI - conductivity about bottom of slab}

si0[ 2 ]:=si[ nPoints ]; {siE - conductivity about top of slab}

si0[ 3 ]:=par0[ 1 ]; {Alfa- ratio of approx.}

mCur:=3;

end;

setApproximationType( apDT ); {approx. type for direct problem}

setApproximationData( si0, mCur ); {approx. data for direct problem}

nApprox := ( nApprox AND $0F ); {XXXXYYYY->0000YYYY}

fHypTg := (( nApprox AND apHypTg ) = apHypTg );

fMulti := (( nApprox AND apExp ) = 0 ) AND NOT fHypTg; {It's not an EXP approx.}

if fMulti then

begin

for i:=1 to nPoints do

begin

Gr[ 1,i ]:=SiMax[ i ];

Gr[ 2,i ]:=SiMin[ i ];

Rg[ i ]:=( Gr[ 1,i ] + Gr[ 2,i ] )/2; {zero estimate of SI}

Rgs[ i ]:=1E33; {biggest integer}

end;

mLast:=nPoints; {loop for every node of approx.}

mCur :=1; {to begin from the only node of approx}

end

else

if fHypTg then

begin

Gr[ 1,1 ]:= siMax[ 1 ]; Gr[ 2,1 ]:= siMin[ 1 ]; Rgs[ 1 ]:=1E33;

Gr[ 1,2 ]:=parMax[ 2 ]; Gr[ 2,2 ]:=parMin[ 2 ]; Rgs[ 2 ]:=1E33;

Gr[ 1,3 ]:=parMax[ 3 ]; Gr[ 2,3 ]:=parMin[ 3 ]; Rgs[ 3 ]:=1E33;

Gr[ 1,4 ]:=parMax[ 4 ]; Gr[ 2,4 ]:=parMin[ 4 ]; Rgs[ 4 ]:=1E33;

for i:=1 to 4 do Rg[ i ]:=( Gr[ 1,i ] + Gr[ 2,i ] )/2;

mLast:=1;

mCur:=4;

end

else

begin

Gr[ 1,1 ]:= siMax[1]; Gr[2,1]:= siMin[1]; Rgs[ 1 ]:=1E33;

Gr[ 1,2 ]:= siMax[nPoints]; Gr[2,2]:= siMin[nPoints]; Rgs[ 2 ]:=1E33;

Gr[ 1,3 ]:= parMax[1]; Gr[2,3]:= parMin[1]; Rgs[ 3 ]:=1E33;

for i:=1 to 3 do Rg[ i ]:=( Gr[ 1,i ] + Gr[ 2,i ] )/2;

mLast:=1;

mCur :=3;

end;

initConst( nLayers, parMaxH, parMaxX , parEps, parEqlB );{set probe params}

end;

procedure directTask; {emulate voltage measurements [with error]}

begin

for i:=1 to nFreqs do

begin

getVoltage( freqs[i], Umr[ i ], Umi[ i ] ); {"measured" Uvn*}

if ( epsU > 0 ) then {add measurement error}

begin

randomize; Umr[ i ]:=Umr[ i ]*( 1 + epsU*( random-0.5 ) );

randomize; Umi[ i ]:=Umi[ i ]*( 1 + epsU*( random-0.5 ) );

end;

writeln('* Voltage measurements have been emulated');

setApproximationType( nApprox ); {approx. type for inverse problem}

setApproximationData( Rg, mCur ); {approx. data for inverse problem}

end;

procedure reduceSILimits; {evaluate SI for m+1 points of approx. using aG}

var

x0, x1, xL, dx, Gr1, Gr2 : real;

j, k : byte;

begin

{----------------------------- get SI min/max for m+1 points of approximation}

dx:=1/( nPoints-1 );

for i:=1 to m+1 do

begin

k:=1;

x1:=0;

x0:=( i-1 )/m;

for j:=1 to nPoints-1 do

begin

xL:=( j-1 )/( nPoints-1 );

if( ( xL < x0 ) AND ( x0 <= xL+dx ) )then

begin

k:=j;

x1:=xL;

end;

Gr[ 1,i ]:=siMax[ k ] + ( siMax[ k+1 ]-siMax[ k ] )*( x0-x1 )/dx;

Gr[ 2,i ]:=siMin[ k ] + ( siMin[ k+1 ]-siMin[ k ] )*( x0-x1 )/dx;

end;

{------------------------------------- get SI for m+1 points of approximation}

for i:=1 to m+1 do

begin

Rg[i]:=getSiFunction( (i-1)/m );

if ( Rg[i] > Gr[1,i] )then Rg[i]:=Gr[1,i];

if ( Rg[i] < Gr[2,i] )then Rg[i]:=Gr[2,i];

if m > 1 then {There're more than 1 point of approx.}

begin

Gr1:= Rg[i]+( Gr[1,i]-Rg[i] )*aG; {reduce upper bound}

Gr2:= Rg[i]-( Rg[i]-Gr[2,i] )*aG; {reduce lower bound}

if ( Gr1 < Gr[1,i] )then Gr[1,i]:=Gr1; {test overflow}

if ( Gr2 > Gr[2,i] )then Gr[2,i]:=Gr2;

end;

setApproximationData( Rg , m+1 );

end;

procedure resultMessage; {to announce new results}

begin

if fMulti then

begin

writeln(' current nodal values of conductivity');

write(' si : ');for i:=1 to m do write(Rg[i] :6:3,' ');writeln;

write(' max: ');for i:=1 to m do write(Gr[1,i]:6:3,' ');writeln;

write(' min: ');for i:=1 to m do write(Gr[2,i]:6:3,' ');writeln;

end

else

begin

for i:=1 to nPoints do si[i]:=getSiFunction( ( i-1 )/( nPoints-1 ) );

if fHypTg then

saveHypTgResults

else

saveExpResults;

end;

procedure clockMessage; {user-friendly message}

begin

writeln('***********************************************************');

write( '* approximation points number :',m:3,' * Time '); clock;

writeln('***********************************************************');

end;

procedure done; {final message}

begin

Sound(222); Delay(111); Sound(444); Delay(111); NoSound; {beep}

write('* Task processing time '); clock; saveTime;

writeln('* Status: Inverse problem has been successfully evaluated.');

end;

Begin

about;

loadData;

initParameters;

directTask;

for m:=1 to mLast do

begin

if fMulti then

begin

mCur:=m;

clockMessage;

end;

doMinimization; {main part of work}

setApproximationData( Rg, mCur ); {set new approx. data}

resultMessage;

if(( fMulti )AND( m < nPoints ))then reduceSILimits;

end;

done;

End.

П1.5 Модуль глобальных данных EData

Unit EData;

Interface

Uses DOS;

Const

maxPAR = 40; {nodes of approximation max number}

maxFUN = 20; {excitation frequencies max number}

maxSPC = 4; {support approximation values number}

iterImax = 50; {max number of internal iterations}

Const

apSpline = 1; {approximation type identifiers}

apHypTg = 3;

apExp = 2;

apPWCon = 4;

apPWLin = 8;

Type

Parameters = array[ 1..maxPAR ] of real; {Si,Mu data}

Functionals = array[ 1..maxFUN ] of real; {Voltage data}

SpecialPar = array[ 1..maxSPC ] of real; {Special data}

Var

hThick : real; {thickness of slab}

nPoints : integer; {nodes of approximation number within [ 0,H ]}

nLayers : integer; {number of piecewise constant SI layers within[ 0,H ]}

nFreqs : integer; {number of excitation frequencies}

nStab : integer; {required number of true digits in SI estimation}

epsU : real; {relative error of measurement Uvn*}

nApprox : byte; {approximation type identifier}

incVal : real; {step for numerical differ.}

parMaxH : integer; {max number of integration steps}

parMaxX : real; {upper bound of integration}

parEps : real; {error of integration}

derivType: byte; {1 for right; 2 for central}

Var

freqs : Functionals; {frequencies of excitment Uvn*}

Umr, Umi : Functionals; { Re(Uvn*),Im(Uvn*) for "measured" data}

Uer, Uei : Functionals; { Re(Uvn*),Im(Uvn*) for estimated data}

mu : Parameters; {relative permeability nodal values}

si, si0 : Parameters; {conductivity approximation nodal values}

siMin, siMax : Parameters; {conductivity nodal values restrictions}

par0 : SpecialPar; {alfa,si2,beta,gamma - for exp&HypTg}

parMin,parMax: SpecialPar; - min/max

zLayer : Parameters; {relative borders of slab layers [0,1]}

Var

aG : real; {scale factor for SImin/max}

Ft : real; {current discrepancy functional value}

fMulti : boolean; {TRUE if it isn't an EXP-approximation}

fHypTg : boolean; {TRUE for Hyperbolic tg approximation}

parEqlB : boolean; {TRUE if b[i]=const}

mCur : integer; {current number of approximation nodes}

inFileName : string; {data file name}

outFileName : string; {results file name}

Var

Rg : Parameters; {current SI estimation}

RgS : Parameters; {previous SI estimation}

Gr : array [ 1..2 ,1..maxPAR ] of real; {SI max/min}

Fh : array [ 1..maxPAR , 1..maxFUN ] of real; {current discrepancies}

Type

TTime=record

H, M, S, S100 : word; {hour,min,sec,sec/100}

end;

Var

clk1, clk2 : TTime; {start&finish time}

procedure loadData; {load all user-defined data from file}

Implementation

procedure loadData;

var

i,eqlB : integer;

FF : text;

begin

assign( FF, outFileName ); {clear output file}

rewrite( FF );

close( FF );

assign( FF, inFileName ); {read input file}

reset( FF );

readln( FF );

readln( FF, hThick, nPoints, nLayers, nFreqs, nStab, epsU, aG, nApprox );

readln( FF );

for i:=1 to nFreqs do read( FF, freqs[i] );

readln( FF );

for i:=1 to nPoints do readln( FF, si[i], siMin[i], siMax[i] );

readln( FF );

readln( FF , incVal, parMaxH, parMaxX, parEps, derivType, eqlB );

readln( FF );

for i:=1 to maxSPC do readln( FF, par0[i] , parMin[i] , parMax[i] );

readln( FF );

if ( eqlB=0 )then

begin

for i:=1 to nLayers+1 do read( FF, zLayer[i] );

parEqlB:=false;

end

else parEqlB:=true;

close( FF );

for i:=1 to maxPAR do mu[i]:=1;

end;

Var

str : string;

Begin

if( ParamCount = 1 )then str:=ParamStr(1)

else

begin

write('Enter I/O file name, please: ');

readln( str );

end;

inFileName :=str+'.txt';

outFileName:=str+'.lst';

End.

П1.6 Модуль работы с файлами EFile

Unit EFile;

Interface

Uses

DOS, EData;

function isStable( ns : integer; var RG1,RG2 ) : boolean;

function saveResults( ns,iter : integer ) : boolean;

procedure saveExpResults;

procedure saveHypTgResults;

procedure clock;

procedure saveTime;

Implementation

Var

FF : text;

i : byte;

function decimalDegree( n:integer ) : real;{10^n}

var

s:real; i:byte;

begin

s:=1;

for i:=1 to n do s:=s*10;

decimalDegree:=s;

end;

function isStable( ns:integer ; var RG1,RG2 ) : boolean;

var

m : real;

R1 : Parameters absolute RG1;

R2 : Parameters absolute RG2;

begin

isStable:=TRUE;

m:=decimalDegree( ns-1 );

for i:=1 to mCur do

begin

if NOT(( ABS( R2[i]-R1[i] )*m )<=ABS( R2[i]) ) then isStable:=FALSE;

RgS[i]:=Rg[i];

end;

function saveResults( ns , iter : integer ) : boolean;

var

sum : real;

begin

sum:=0;

for i:=1 to nFreqs do sum:=sum + Fh[1,i];

sum:=SQRT( sum/nFreqs );

assign( FF , outFileName );

append( FF );

write( iter:2, ' <”>', sum:10:7, ' Rg=' );

write( FF , iter:2, ' <”>', sum:10:7, ' Rg=');

for i:=1 to mCur do

begin

write( Rg[i]:6:3, ' ');

write( FF , Rg[i]:6:3, ' ');

end;

writeln;

writeln( FF );

close( FF );

saveResults:=isStable( ns , Rgs , Rg );

end;

procedure saveExpResults;

begin

assign( FF , outFileName );

append( FF );

writeln( ' siE=',Rg[2]:6:3,' siI=',Rg[1]:6:3,' alfa=',Rg[3]:6:3);

writeln( FF , ' siE=',Rg[2]:6:3,' siI=',Rg[1]:6:3,' alfa=',Rg[3]:6:3);

write( ' SI: ');

write( FF , ' SI: ');

for i:=1 to nPoints do

begin

write( si[i]:6:3,' ');

write( FF , si[i]:6:3,' ');

end;

writeln;

writeln( FF );

close( FF );

end;

procedure saveHypTgResults;

begin

assign( FF , outFileName );

append( FF );

writeln( ' si1=',Rg[2]:6:3,' si2=',Rg[1]:6:3,' beta=',Rg[3]:6:3,' gamma=',Rg[4]:6:3);

writeln( FF , ' si1=',Rg[2]:6:3,' si2=',Rg[1]:6:3,' beta=',Rg[3]:6:3,' gamma=',Rg[4]:6:3);

write( ' SI: ');

write( FF , ' SI: ');

for i:=1 to nPoints do

begin

write( si[i]:6:3,' ');

write( FF , si[i]:6:3,' ');

end;

writeln;

writeln( FF );

close( FF );

end;

procedure clock; {t2 = t2-t1}

var

H1,M1,S1,H2,M2,S2,sec1,sec2 : longint;

begin

GetTime( clk2.H, clk2.M, clk2.S, clk2.S100 ); {current time}

H2:=clk2.H; M2:=clk2.M; S2:=clk2.S; H1:=clk1.H; M1:=clk1.M; S1:=clk1.S;

sec2:= ( H2*60 + M2 )*60 + S2;

sec1:= ( H1*60 + M1 )*60 + S1;

if( sec2 < sec1 )then sec2:=sec2 + 85020; {+23.59.59}

sec2:=sec2 - sec1;

clk2.H := sec2 div 3600; sec2:=sec2 - clk2.H*3600;

clk2.M := sec2 div 60; sec2:=sec2 - clk2.M*60;

Страницы: 1, 2, 3, 4