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

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;

clk2.S := sec2;

writeln( clk2.H:2, ':', clk2.M:2, ':', clk2.S:2 );

end;

procedure saveTime;

begin

assign( FF , outFileName );

append( FF );

write( FF ,'* Processing time ',clk2.H:2, ':', clk2.M:2, ':', clk2.S:2 );

close( FF );

end;

End.

П1.7 Модуль решения прямой задачи ВТК для НВТП EDirect

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

{ Estimates Uvn* for Eddy current testing of inhomogeneous multilayer slab }

{ with surface( flat ) probe. }

{ It can do it using one of five types of conductivity approximation : }

{Spline, Exponential, Piecewise constant, Piecewise linear,Hyperbolic tangent}

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

{$F+}

Unit EDirect;

Interface

Uses EData, EMath;

Type

siFunc = function( x:real ) : real;

Var

getSiFunction : siFunc; {for external getting SI estimate}

procedure initConst( par1,par2:integer; par3,par4:real; par5:boolean );

procedure getVoltage( freq : real ; var ur,ui : real ); { Uvn* = ur + j*ui }

procedure setApproximationType( approx : byte ); { type of approx. }

procedure setApproximationItem( SIG:real ; N : byte ); { set SIGMA[ N ]}

procedure setApproximationData( var SIG; nVal : byte ); { SIGMA[1..nVal] }

procedure getApproximationData( var SIG ; var N : byte ); { get SIGMA[ N ]}

Implementation

Const

PI23 = 2000*pi; {2*pi*KHz}

mu0 = 4*pi*1E-7; {magnetic const}

Var

appSigma : Parameters; {conductivity approximation data buffer}

appCount : byte; {size of conductivity approximation data buffer}

appType : byte; {conductivity approximation type identifier}

Type

commonInfo=record

w : real; {cyclical excitation frecuency}

R : real; {equivalent radius of probe}

H : real; {generalized lift-off of probe}

Kr : real; {parameter of probe}

eps : real; {error of integration}

xMax : real; {upper bound of integration}

steps : integer; {current number of integration steps}

maxsteps: integer; {max number of integration steps}

Nlay : integer; {number of layers in slab}

sigma : Parameters; {conductivity of layers}

m : Parameters; {relative permeability of layers}

b : Parameters; {thickness of layers}

zCentre : Parameters; {centre of layer}

end;

procFunc = procedure( x:real; var result:complex);

Var

siB, siC, siD : Parameters; {support for Spline approx.}

cInfo : commonInfo; {one-way access low level info}

function siSpline( x:real ) : real;{Spline approximation}

begin

if( appCount = 1 )then

siSpline := appSigma[ 1 ]

else

siSpline:=Seval( appCount, x, appSigma, siB, siC, siD);

end;

function siExp( x:real ) : real;{Exponential approximation}

begin

siExp:=(appSigma[2]-appSigma[1])*EXP( -appSigma[3]*(1-x) ) + appSigma[1];

end;

function siPWConst( x:real ) : real;{Piecewise constant approximation}

var

dx, dh : real; i : byte;

begin

if( appCount = 1 )then siPWConst := appSigma[ 1 ]

else

begin

dh:=1/( appCount-1 );

dx:=dh/2;

i:=1;

while( x > dx ) do

begin

i:=i + 1;

dx:=dx + dh;

end;

siPWConst:=appSigma[ i ];

end;

function siPWLinear( x:real ) : real;{Piecewise linear approximation}

var

dx, dh : real;

i : byte;

begin

if( appCount = 1 )then siPWLinear := appSigma[ 1 ]

else

begin

dh:=1/( appCount-1 );

dx:=0;

i:=1;

repeat

i:=i + 1;

dx:=dx + dh;

until( x <= dx );

siPWLinear:=appSigma[i-1]+( appSigma[i]-appSigma[i-1] )*( x/dh+2-i);

end;

function siHyperTg( x:real ) : real;{Hyperbolic tangent approximation}

begin

siHyperTg:=appSigma[2]+(appSigma[1]-appSigma[2])*(1+th((appSigma[3]-x)/appSigma[4]))/2;

end;

procedure setApproximationType( approx : byte );

begin

appType := approx;

write('* conductivity approximation type : ');

case approx of

apSpline : begin

writeln('SPLINE');

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

МЕНЮ

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