Рефераты. Решение обратной задачи вихретокового контроля
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
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
begin
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 ); 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].
