HThick
|
- толщина пластины,[мм]
|
nPoints
|
- количество узлов аппроксимации электропроводности
для PWL,PWC,SPL аппроксимаций. В случае EXP,HTG аппроксимации вычисление
значений ЭП в них производится по окончании расчетов
|
nLayers
|
- количество интервалов с кусочно-постоянной
электропроводностью, на которые разбивается пластина для непосредственного
расчета вносимой ЭДС по реккурентным формулам для многослойной пластины
|
nFreqs
|
- количество частот возбуждения гармоник вносимой
относительной ЭДС
|
nStab
|
- число стабилизируемых значащих цифр
|
epsU
|
- погрешность измерения
|
aG
|
- коэффициент сжатия ограничений (aG<=1); НЕ
используется при EXP,HTG аппроксимации
|
nApprox
|
- типы аппроксимации прямой и обратной задач
|
si
|
- значения проводимости в узлах аппроксимации
|
siMin, siMax
|
- ограничения на возможные значения проводимости в
узлах аппроксимации в процессе решения ОЗ
|
incVal
|
- величина dx для численного дифференцирования
|
maxSteps
|
- максимальное число отрезков интегрирования
|
maxX
|
- верхний предел интегрирования при расчете Uvn*
|
Eps
|
- погрешность интегрирования при расчете Uvn*
|
dType
|
- тип разностной производной (=1 правая или =2
центральная)
|
eqlB
|
- толщины слоев пластины одинаковы( b=hThick/nLayers)
если eqlB>0, в противном случае используются
координаты слоев из файла
|
П1.2 Используемые аппроксимации
Примечание. Координата ХÎ[0,1]
отсчитывается от дна пластины для всех аппроксимаций.
Сплайн(SPL), кусочно-линейная(PWL), кусочно-постоянная(PWC) аппроксимации.
В процессе расчетов ищутся значения электропроводности
в узлах аппроксимации, причем количество узлов увеличивается от едениицы до nPoints в целях сохранения устойчивости решения.
Начальные значения (узловые значения ²истинной² ЭП для
эмуляции измерений U*вн) задаются в столбце ²si² файла исходных данных, начальные значения ограничений
на узловые значения ЭП в столбцах ²siMin² и ²siMax²(движение
по столбцу сверху вниз соответствует изменению координаты от дна пластины до
обрабатываемой повехности).
Экспоненциальная
аппроксимация(EXP)
В случае задания экспоненциальной аппроксимации
зависимость электропрводности от толщины представляется в виде
SIGMA = ( siE-siI )*EXP( -alfa*(1-x) ) + siI
Варьруемыми параметрами являются эектропроводность на
верхней поверхности siЕ, электропроводность "на бесконечности"
siI и параметр alfa. В файле исходных данных в таблице из nPoints
строк с подзаголовком "si siMin siMax", информация об ограничениях
на параметры siE, siI задается в первой и nPoints-строке.
Величина и ограничения для параметра alfa задаются первой строкой в
"special approximation parameters".
Аппроксимация
гиперболическим тангенсом (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. *)
(*
(C) 1999 by Nikita U.Dolgov *)
(* 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;
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;
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;
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;
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;
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 );
readln( FF,
hThick, nPoints, nLayers, nFreqs, nStab, epsU, aG, nApprox );
readln( FF );
for i:=1 to
nFreqs do read( FF, freqs[i] );
readln( FF );
readln( FF );
readln( FF );
for i:=1 to
nPoints do readln( FF, si[i], siMin[i], siMax[i] );
readln( FF );
readln( FF );
readln( FF ,
incVal, parMaxH, parMaxX, parEps, derivType, eqlB );
readln( FF );
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;
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,' ');
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
{****************************************************************************}
{ ERIN submodule :
EDirect , 15.02.99, (C) 1999 by Nikita U.Dolgov }
{****************************************************************************}
{ 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;
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;
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');
getSiFunction := siSpline;
end;
apExp : begin
writeln('EXP');
getSiFunction := siExp;
end;
apPWCon : begin
writeln('PIECEWISE CONST');
getSiFunction := siPWConst;
end;
apPWLin : begin
writeln('PIECEWISE LINEAR');
getSiFunction := siPWLinear;
end;
apHypTg : begin
writeln('HYPERBOLIC TANGENT');
getSiFunction := siHyperTg;
end;
end;
end;
procedure
setApproximationData( var SIG ; nVal : byte );
var
Sigma :
Parameters absolute SIG; i:byte;
begin
appCount :=
nVal;
for i:=1 to
nVal do appSigma[ i ]:=Sigma[ i ];
if( appType =
apSpline )then Spline( appCount, appSigma, siB, siC, siD);
end;
procedure
getApproximationData( var SIG ; var N : byte );
var
Sigma :
Parameters absolute SIG; i:byte;
begin
N :=
appCount;
for i:=1 to
appCount do Sigma[ i ]:=appSigma[ i ];
end;
procedure
setApproximationItem( SIG:real ; N : byte );
begin
appSigma[ N ]
:= SIG;
if( appType =
apSpline )then Spline( appCount, appSigma, siB, siC, siD);
end;
procedure
functionFi( x:real ; var result:complex );{get boundary conditions function
value}
var
beta : array[
1..maxPAR ]of real;
q : array[
1..maxPAR ]of complex;
fi : array[
0..maxPAR ]of complex;
th , z1 , z2
, z3 , z4 , z5 , z6 , z7 : complex;
i : byte;
begin
mkComp( 0, 0,
fi[0] );
with cInfo do
for i:=1 to
Nlay do
begin
beta[i]:=R*sqrt( w*mu0*sigma[i] ); {calculation of beta}
mkComp(
sqr(x), sqr(beta[i])*m[i], z7 ); {calculation of q, z7=q^2}
SqrtC(
z7, q[i] );
mulComp(
q[i], b[i], z6 ); {calculation of th,z6=q*b}
tanH(
z6, th ); {th=tanH(q*b)}
mkComp(
sqr(m[i])*sqr(x), 0, z6 ); {z6=m2x2}
SubC(
z6, z7, z5); {z5=m2x2-q2}
AddC(
z6, z7, z4); {z4=m2x2+q2}
MulC(
z5, th, z1); {z1=z5*th}
MulC(
z4, th, z2); {z2=z4*th}
mulComp(
q[i], 2*x*m[i], z3 ); {z3=2xmq}
SubC(
z2, z3, z4 );
MulC(
z4, fi[i-1], z5 );
SubC(
z1, z5, z6 ); {z6=high}
AddC(
z2, z3, z4 );
MulC(
z1, fi[i-1], z5 );
SubC(
z4, z5, th ); {th=low}
DivC(
z6, th, fi[i] );
end;
eqlComp(
result, fi[ cInfo.Nlay ] );
end;
procedure
funcSimple( x:real; var result:complex );{intergrand function value}
var
z : complex;
begin
with cInfo do
begin
functionFi( x, result );
mulComp(
result, exp( -x*H ), z );
mulComp(
z, J1( x*Kr ), result );
mulComp(
result, J1( x/Kr ), z );
eqlComp(
result, z );
end;
end;
procedure funcMax(
x:real; var result:complex );{max value;
When Fi[Nlay]=1}
var
z1, z2 :
complex;
begin
with cInfo do
begin
mkComp(
1,0,z1 );
mulComp(
z1, exp(-x*H), z2 );
mulComp(
z2, J1( x*Kr ), z1 );
mulComp(
z1, J1( x/Kr ), result );
end;
end;
procedure
integralBS( func:procFunc ; var result:complex );{integral by Simpson}
var
z , y , tmp :
complex;
hh :
real;
i, iLast :
word;
begin
with cInfo do
begin
hh:=xMax/steps;
iLast:=steps div 2;
nullComp(tmp);
func( 0,
z );
eqlComp(
y, z );
for i:=1
to iLast do
begin
func( ( 2*i-1 )*hh , z );
deltaComp( tmp, z );
end;
mulComp(
tmp, 4, z );
deltaComp( y, z );
nullComp( tmp );
iLast:=iLast-1;
for i:=1
to iLast do
begin
func( 2*i*hh ,z );
deltaComp( tmp, z );
end;
mulComp(
tmp, 2, z );
deltaComp( y, z );
func(
xmax, z );
deltaComp(y,z);
mulComp(
y, hh/3, z );
eqlComp(
result, z );
end;
end;{I = h/3 * [
F0 + 4*sum(F2k-1) + 2*sum(F2k) + Fn ]}
procedure
integral( F:procFunc; var result:complex );{integral with given error}
var
e , e15 :
real;
flag :
boolean;
delta ,
integ1H , integ2H : complex;
begin
with cInfo do
begin
e15
:=eps*15;{ | I2h-I1h |/(2^k -1 )< Eps ; k=4 for Simpson method}
steps:=20;
flag
:=false;
integralBS( F, integ2H );
repeat
eqlComp( integ1H, integ2H );
steps:=steps*2;
integralBS(
F, integ2H );
SubC( integ2H, integ1H, delta );
e:=Leng( delta );
if( e<e15 )then flag:=true;
until( (
flag ) OR (steps>maxsteps) );
if( flag
)then
begin
eqlComp( result, integ2H );
end
else
begin
writeln('Error: Too big number of integration steps.');
halt(1);
end;
end;
end;
procedure
initConst( par1, par2 : integer; par3, par4 : real; par5:boolean );
var
i : byte;
bThick, dl, x
: real;
const
Ri=0.02;
hi=0.005; { radius and lift-off of excitation coil}
Rm=0.02;
hm=0.005; { radius and lift-off of measuring coil}
begin
with cInfo do
begin
Nlay
:=par1;
xMax
:=par3;
maxsteps:=par2;
R
:=sqrt( Ri*Rm );
H
:=( hi+hm )/R;
Kr
:=sqrt( Ri/Rm );
eps
:=par4;
bThick
:=hThick*0.002/R; {2*b/R [m]}
for i:=1
to Nlay do m[i]:= mu[i];
if par5
then
begin
bThick:=bThick/NLay;
for
i:=1 to Nlay do b[i]:=bThick;
dl:=1/NLay;
x:=dl/2; {x grows up from bottom of slab to the top}
for
i:=1 to Nlay do
begin
zCentre[i]:=x;
x:=x + dl;
end;
end
else
for
i:=1 to Nlay do
b[i]:=( zLayer[i+1]-zLayer[i] )*bThick;
zCentre[i]:=( zLayer[i+1]+zLayer[i] )/2;
end;
end;
end;
procedure init(
f:real );{get current approach of conductivity values}
var
i : byte;
begin
with cInfo do
begin
w:=PI23*f;
for i:=1
to Nlay do sigma[i]:=getSiFunction( zCentre[i] )*1E6;
end;
end;
procedure
getVoltage( freq : real ; var ur,ui : real );
var
U, U0, Uvn,
tmp : complex;
begin
init( freq );
integral(
funcSimple, U ); { U =Uvn }
integral(
funcMax , U0 ); { U0=Uvn max }
divComp( U,
Leng(U0), Uvn ); { Uvn=U/|U0| }
mkComp( 0, 1,
tmp ); { tmp=( 0+j1 ) }
MulC( tmp,
Uvn, U ); { U= j*Uvn = Uvn* }
ur := U.re;
ui := U.im;
end;
END.
П1.8 Модуль решения обратной задачи ВТК для НВТП EMinimum
Unit EMinimum;
INTERFACE
Uses EData, Crt,
EFile, EDirect;
procedure
doMinimization;
IMPLEMENTATION
procedure
getFunctional( Reg : byte );
var
ur, ui, dur,
dui, Rgt : real;
ur2, ui2:
Functionals;
i, j, k :
byte;
begin
getApproximationData( si , k );
setApproximationData( Rg, mCur );
case Reg of
0 : for
i:=1 to nFreqs do {get functional F value}
begin
getVoltage( freqs[i], ur, ui );
Uer[ i ]:=ur; {we need it for dU}
Uei[ i ]:=ui;
Fh[1,i] := SQR( ur-Umr[i] ) + SQR( ui-Umi[i] );
end;
{Right:U'(i)= (U(i+1)-U(i))/h}
1 : for
i:=1 to mCur do {get dF/dSI[i] value}
begin
Rgt:=Rg[i]*( 1+incVal ); {si[i]=si[i]+dsi[i]}
setApproximationItem( Rgt, i ); {set new si[i] value}
for j:=1 to nFreqs do
begin {get dUr/dSI,dUi/dSI}
getVoltage( freqs[ j ], ur, ui );
dur:=( ur-Uer[j] )/( Rg[i]*incVal );
dui:=( ui-Uei[j] )/( Rg[i]*incVal );
Fh[i,j]:=2*(dur*(Uer[j]-Umr[j])+dui*(Uei[j]-Umi[j]));
end;
setApproximationItem( Rg[i], i ); {restore si[i] value}
end;
{Central:U'(i)= (U(i+1)-U(i-1))/2h}
2 : for
i:=1 to mCur do {get dF/dSI[i] value}
begin
Rgt:=Rg[i]*( 1+incVal ); {si[i]=si[i]+dsi[i]}
setApproximationItem( Rgt, i ); {set new si[i] value}
for j:=1 to nFreqs do getVoltage( freqs[j],ur2[j],ui2[j] );
Rgt:=Rg[i]*( 1-incVal ); {si[i]=si[i]-dsi[i]}
setApproximationItem( Rgt, i ); {set new si[i] value}
for j:=1 to nFreqs do
begin
{get dUr/dSI,dUi/dSI}
getVoltage( freqs[ j ], ur, ui );
dur:=( ur2[j]-ur )/( 2*Rg[i]*incVal );
dui:=( ui2[j]-ui )/( 2*Rg[i]*incVal );
Fh[i,j]:=2*(dur*(Uer[j]-Umr[j])+dui*(Uei[j]-Umi[j]));
end;
setApproximationItem( Rg[i], i ); {restore si[i] value}
end;
end;
setApproximationData( si , k );
end;
procedure
doMinimization;
const
mp1Max =
maxPAR + 1;
mp2Max =
maxPAR + 2;
m2Max = 2*(
maxPAR + maxFUN );
m21Max =
m2Max + 1;
n2Max =
2*maxFUN;
m1Max =
maxPAR + n2Max;
n1Max =
n2Max + 1;
mm1Max =
maxPAR + n1Max;
minDh : real
= 0.001; {criterion of an exit from golden section method}
var
A : array [ 1
.. m1Max , 1 .. m21Max ] of real;
B : array [ 1
.. m1Max] of real;
Sx: array [ 1
.. m21Max] of real;
Zt: array [ 1
.. maxPAR] of real;
Nb: array [ 1
.. m1Max] of integer;
N0: array [ 1
.. m21Max] of integer;
a1, a2, dh,
r, tt, tp, tl, cv, cv1, cl, cp : real;
n2, n1, mp1,
mp2, mm1, m1, m2, m21 : integer;
ain :
real;
apn :
real;
iq :
integer;
k1 :
integer;
n11 :
integer;
ip :
integer;
iterI :
integer;
i,j,k :
integer;
label
102 ,103 ,104
,105 ,106 ,107 ,108;
begin
n2:=2*nFreqs;
n1:=n2+1; m1:=mCur+n2;
mp1:=mCur+1;
mp2:=mCur+2; mm1:=mCur+n1;
m2:=2*(
mCur+nFreqs ); m21:=m2+1;
for k:=1 to
m1Max do
for i:=1
to m21Max do
A[k,i]:=0;
iterI:=0;
102:
iterI:=iterI+1;
getFunctional( 0 );
for i:=1 to
nFreqs do b[i]:= -Fh[1,i];
getFunctional( derivType );
for k:=1 to
mCur do
begin
Zt[k]:=Rg[k];
for i:=1
to nFreqs do
begin
A[i,k+1]:=Fh[k,i];
A[i+nFreqs,k+1]:=-A[i,k+1];
end;
for i:=1
to nFreqs do B[i]:=B[i]+Rg[k]*A[i,k+1];
end;
for i:=1 to
nFreqs do B[i+nFreqs]:=-B[i];
for i:=n1 to
m1 do B[i]:=Gr[1,i-n2]-Gr[2,i-n2];
for i:=1 to
m1 do
begin
if
i<=n2 then
for
k:=2 to mp1 do B[i]:=B[i]-A[i,k]*Gr[2,k-1];
A[i,1]:=-1;
if
i>n2 then
begin
A[i,1]:=0;
for
k:=2 to mp1 do
if i-n2=k-1 then A[i,k]:=1
else A[i,k]:=0;
end;
for
k:=mp2 to m21 do
if
k-mp1=i then A[i,k]:=1
else
A[i,k]:=0;
end;
k1:=1;
for k:=1 to
n2 do
if
B[k1]>B[k] then k1:=k;
for k:=1 to
mp1 do A[k1,k]:=-A[k1,k];
A[k1,mCur+1+k1]:=0;
B[k1]:=-B[k1];
for i:=1 to
n2 do
if
i<>k1 then
begin
B[i]:=B[i]+B[k1];
for
k:=1 to mm1 do A[i,k]:=A[i,k]+A[k1,k];
end;
for i:=mp2 to m21 do
begin
Sx[i]:=B[i-mp1];
Nb[i-mp1]:=i;
end;
for i:=1 to
mp1 do Sx[i]:=0;
Sx[1]:=B[k1];
Sx[mp1+k1]:=0;
Nb[k1]:=1;
103:
for i:=2 to
m21 do N0[i]:=0;
104:
for i:=m21 downto 2 do
if
N0[i]=0 then n11:=i;
for k:=2 to
m21 do
if
((A[k1,n11]<A[k1,k]) AND (k<>N0[k])) then n11:=k;
if
A[k1,n11]<=0 then goto 105;
iq:=0;
for i:=1 to
m1 do
if
i<>k1 then
begin
if
A[i,n11]>0 then
begin
iq:=iq+1;
if iq=1 then
begin
Sx[n11]:=B[i]/A[i,n11]; ip:=i;
end
else
begin
if Sx[n11]>B[i]/A[i,n11] then
begin
Sx[n11]:=B[i]/A[i,n11]; ip:=i;
end;
end;
end
else
if
iq=0 then
begin
N0[n11]:=n11;
goto 104;
end;
end;
Sx[Nb[ip]]:=0;
Nb[ip]:=n11;
B[ip]:=B[ip]/A[ip,n11];
apn:=A[ip,n11];
for k:=2 to
m21 do A[ip,k]:=A[ip,k]/apn;
for i:=1 to
m1 do
if
i<>ip then
begin
ain:=A[i,n11];
B[i]:=-B[ip]*ain+B[i];
for
j:=1 to m21 do A[i,j]:=-ain*A[ip,j]+A[i,j];
end;
for i:=1 to
m1 do Sx[Nb[i]]:=B[i];
goto 103;
105:
for k:=1 to
mCur do Sx[k+1]:=Sx[k+1]+Gr[2,k];
a1:=0;
a2:=1.;
dh:=a2-a1;
r:=0.618033;
tl:=a1+r*r*dh;
tp:=a1+r*dh;
j:=1;
108:
if j=1 then
tt:=tl else tt:=tp;
106:
for i:=1 to
mCur do Rg[i]:=Zt[i]+tt*(Sx[i+1]-Zt[i]);
getFunctional( 0 );
cv:=abs(Fh[1,1]);
if
nFreqs>1 then
for k:=2
to nFreqs do
begin
cv1:=abs(Fh[1,k]);
if
cv<cv1 then cv:=cv1;
end;
if (j=1) or
(j=3) then cl:=cv
else cp:=cv;
if j=1 then
begin
j:=2;
goto
108;
end;
if
dh<MinDh then goto 107;
if cl>cp
then
begin
a1:=tl;
dh:=a2-a1; tl:=tp; tp:=a1+r*dh ; tt:=tp; cl:=cp; j:=4;
end
else
begin
a2:=tp;
dh:=tp-a1; tp:=tl; tl:=a1+r*r*dh; tt:=tl; cp:=cl; j:=3;
end;
goto 106;
107:
if (iterI
< iterImax)AND(NOT saveResults( nStab,iterI )) then goto 102;
end;
End.
Приложение 2 - Удельная электрическая проводимость материалов
Приведем сводку справочных данных согласно[7-9].
Материал
|
smin ,[МСм/м]
|
smax
,[МСм/м]
|
Немагнитные стали
|
0.4
|
1.8
|
Бронзы (БрБ, Бр2, Бр9)
|
6.8
|
17
|
Латуни (ЛС59, ЛС62)
|
13.5
|
17.8
|
Магниевые сплавы (МЛ5-МЛ15)
|
5.8
|
18.5
|
Титановые сплавы (ОТ4, ВТ3-ВТ16)
|
0.48
|
2.15
|
Алюминиевые сплавы (В95, Д16, Д19)
|
15.1
|
26.9
|
Приложение 4 - Abstract
The inverse eddy current problem can be described as the task of
reconstructing an unknown distribution of electrical conductivity from
eddy-current probe voltage measurements recorded as function of excitation
frequency. Conductivity variation may be a result of surface processing with
substances like hydrogen and carbon or surface heating.
Mathematical
reasons and supporting software for inverse conductivity profiling were
developed by us. Inverse problem was solved for layered plane and cylindrical
conductors.
Because
the inverse problem is nonlinear, we propose using an iterative algorithm which
can be formalized as the minimization of an error functional related to the
difference between the probe voltages theoretically predicted by the direct
problem solving and the measured probe voltages.
Numerical
results were obtained for some models of conductivity distribution. It was
shown that inverse problem can be solved exactly in case of correct
measurements. Good estimation of the true conductivity distribution takes place
also for measurement noise about 2 percents but in case of 5 percent error
results are worse.