附程序如下：
unit Nquasidynamic;
interface
uses
Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms,
Dialogs, StdCtrls, ExtCtrls,math;
type
TForm1 = class(TForm)
Label1: TLabel;
Label4: TLabel;
Label5: TLabel;
Label6: TLabel;
Label8: TLabel;
Image1: TImage;
COUNT: TLabel;
COUNT3: TLabel;
Label9: TLabel;
Button1: TButton;
Edit1: TEdit;
Edit4: TEdit;
Edit5: TEdit;
Edit6: TEdit;
Edit8: TEdit;
Button2: TButton;
OpenDialog1: TOpenDialog;
derivative: TLabel;
Edit2: TEdit;
Label2: TLabel;
Edit3: TEdit;
Label3: TLabel;
procedure Button1Click(Sender: TObject);
procedure Button2Click(Sender: TObject);
procedure FormCreate(Sender: TObject);
private
{ Private declarations }
public
{ Public declarations }
end;
var
Form1: TForm1;
ex:boolean;
implementation
{$R *.dfm}
procedure TForm1.Button1Click(Sender: TObject);
var
t,tau,theta,vel,disp:double;
tauh,thetah,velh,disph:double;
tau2,theta2,vel2,disp2,tauk:double; //tauk: the spring force
taui,thetai,veli,dispi:double;
vc,v0,dc:double; //vc: dimensionless steady state velocity
tauold,dispold:double;
R_tau,R_theta,R_vel:double;
R_tau2,R_theta2,R_vel2:double;
twoquasi,twoquasih:double; // 状态演化方程的第二项表达式
delta,deltai,deltah,delta2:double;
Tsquare:double;
quasi:boolean;
i,m,n:integer;
B_A,kcr,k,alfa:double;
dispmax:double;
h,hh,eps:double;
fname:string; // the tfile variants
data:textfile;
px,py,pz:longint; //plot variants
color:tcolor;
{.........................N chi equation.........................................}
n_index:double; //n次方
const
Tsystem:double=5;
maxx:double=300;
minx:double=0;
maxy:double=50;
miny:double=-50;
begin
if opendialog1.Execute then
begin
fname:=opendialog1.FileName ;
assignfile(data,fname);
rewrite(data);
end
else
exit;
writeln(data,'time',' ','disp',' ','delta0',' ','tau',' ','theta',' ','vel',' ','state');
m:=0;
i:=0;
n:=0; //i,n the count value
dispi:=0; //initial mass displacement
t:=0; //time
taui:=0; //initial shear stress
thetai:=2; //initial state valve
alfa:=strtofloat(edit6.Text); //critical stiffness times临界刚度的倍数
dc:=strtofloat(edit8.Text); //characteristic displacement
vc:=strtofloat(edit4.text); //initial steady state velocity dimentionless =v*/dc
veli:=vc; //initial mass velocity
v0:=strtofloat(edit5.text)*vc; // load poit velocity
dispmax:=strtofloat(edit1.Text); //displacement upper limit
B_A:=strtofloat(edit2.text); //b/a
kcr:=B_A-1;
k:=alfa*kcr; //the stiffness
h:=5.0E-5*(1/veli); //step size
eps:=1E-8; // the precision
tauold:=0;
dispold:=0;
deltai:=0;
Tsquare:=sqr(2*pi/Tsystem);
quasi:=true;
n_index:=strtofloat(edit3.Text);
image1.Canvas.Pen.color:=clred; //plot a line
form1.Image1.canvas.MoveTo(0,trunc(Image1.height/2));
form1.Image1.Canvas.lineTo(Image1.width,trunc(Image1.height/2));
repeat
{...............................quasistatic motion portion..................................}
twoquasi:=vc*power(veli*thetai/vc/n_index,n_index); //状态演化方程的第二项表达式
if ((veli*Tsystem/2/pi)<5e-4) and quasi then
begin
R_tau:=k*(v0-veli);
R_theta:=vc-twoquasi;
R_vel:=veli*(k*(v0-veli)-B_A/thetai*(vc-twoquasi));
repeat
h:=h/2;
hh:=h/2;
tau:=taui+h*R_tau;
tauh:=taui+hh*R_tau;
theta:=thetai+h*R_theta;
thetah:=thetai+hh*R_theta;
vel:=veli+h*R_vel;
velh:=veli+hh*R_vel;
disp:=dispi+h*veli;
disph:=dispi+hh*veli;
delta:=deltai+h*v0;
deltah:=deltai+hh*v0;
twoquasih:=vc*power(velh*thetah/vc/n_index,n_index);
R_tau2:=k*(v0-velh);
R_theta2:=vc-twoquasih;
R_vel2:=velh*(k*(v0-velh)-B_A/thetah*(vc-twoquasih));
tau2:=tauh+hh*R_tau2;
theta2:=thetah+hh*R_theta2;
vel2:=velh+hh*R_vel2;
disp2:=disph+hh*velh;
delta2:=deltah+hh*v0;
until (abs(tau2-tau)<=abs(eps*tau2)) and (abs(theta2-theta)<=abs(eps*theta2)) and
(abs(vel2-vel)<=abs(eps*vel2)) and (abs(disp2-disp)<=abs(eps*disp2))and (abs(delta2-delta)<=abs(eps*delta2));
end
{.........................dynamic motion portion..............................}
else
begin
quasi:=false;
R_theta:=vc-twoquasi;
R_tau:=1/veli*Tsquare*(deltai-dispi-taui/k)+B_A/thetai*(vc-twoquasi);
R_vel:=Tsquare*(deltai-dispi-taui/k);
repeat
h:=h/2;
hh:=h/2;
tau:=taui+h*R_tau;
tauh:=taui+hh*R_tau;
vel:=veli+h*R_vel;
velh:=veli+hh*R_vel;
theta:=thetai+h*R_theta;
thetah:=thetai+hh*R_theta;
disp:=dispi+h*veli;
disph:=dispi+hh*veli;
delta:=deltai+h*v0;
deltah:=deltai+hh*v0;
twoquasih:=vc*power(velh*thetah/vc/n_index,n_index);
R_theta2:=vc-twoquasih;
R_tau2:=1/velh*Tsquare*(deltah-disph-tauh/k)+B_A/thetah*(vc-twoquasih);
R_vel2:=Tsquare*(deltah-disph-tauh/k);
tau2:=tauh+hh*R_tau2;
vel2:=velh+hh*R_vel2;
theta2:=thetah+hh*R_theta2;
disp2:=disph+hh*velh;
delta2:=deltah+hh*v0;
until(abs(tau2-tau)<=abs(eps*tau2)) and (abs(theta2-theta)<=abs(eps*theta2)) and
(abs(vel2-vel)<=abs(eps*vel2)) and (abs(disp2-disp)<=abs(eps*disp2))and (abs(delta2-delta)<=abs(eps*delta2));
if abs(1-k*(delta2-disp2)/tau2)<=1e-8 then quasi:=true; //一个判定条件，在动态时overshoot弹簧力=剪切力时转到静态。
end;
tauk:=k*(delta2-disp2);
t:=t+h; //t plus the final step size
if (abs(disp2-dispold)>=0.05) or (abs(tau2-tauold)>=0.05) then //the disp and tau record interval
begin
px:=trunc(form1.Image1.Width/(maxx-minx)*(disp2-minx)); // zoom in or zoom out
py:=trunc(form1.Image1.height-form1.Image1.height/(maxy-miny)*(tau2-miny));
pz:=trunc(form1.Image1.height-form1.Image1.height/(maxy-miny)*(tauk-miny));
if (veli*Tsystem/2/pi<5e-4) and quasi then
begin
writeln(data,t,' ',disp2,' ',delta2,' ',tau2,' ',theta2,' ',vel2,' ','quasistatic',' ',R_vel,' ',R_theta,' ',R_tau); //write file
image1.canvas.Pixels[px,py]:=clblack;
end
else
begin
writeln(data,t,' ',disp2,' ',delta2,' ',tau2,' ',theta2,' ',vel2,' ','dynamic');
image1.canvas.Pixels[px,py]:=clred;
image1.Canvas.Pixels[px,pz]:=clblack;
end;
dispold:=disp2;
tauold:=tau2;
Application.ProcessMessages;
end;
h:=5.0E-5/vel2;
taui:=tau2;
thetai:=theta2;
dispi:=disp2;
deltai:=delta2;
veli:=vel2;
until (dispi>dispmax) or ex;
closefile(data);
end;
procedure TForm1.Button2Click(Sender: TObject);
begin
ex:=true;
end;
procedure TForm1.FormCreate(Sender: TObject);
begin
ex:=false;
end;
end.
7 个解决方案
#1
运行时问题提示是出在dynamic portion 之后的twoquasih:=vc*power(velh*thetah/vc/n_index,n_index);
#2
设一下断点看看那几个变量的值是不是有0出现了
#3
无效的浮点操作？
n_index:=strtofloat(edit3.Text);
这句设个断点，看看赋值内容如何
n_index:=strtofloat(edit3.Text);
这句设个断点，看看赋值内容如何
#4
我也觉得这里要下断点
或者楼主试试 加默认值
strtofloatdef(edit3.Text,1);
或者楼主试试 加默认值
strtofloatdef(edit3.Text,1);
#5
问题出在你的power()函数，请仔细看看它的参数的数据类型。
#6
n_index:double;
这句 类型申明 改成Real类型 看看
这句 类型申明 改成Real类型 看看
#7
断点跟踪，法宝
#1
运行时问题提示是出在dynamic portion 之后的twoquasih:=vc*power(velh*thetah/vc/n_index,n_index);
#2
设一下断点看看那几个变量的值是不是有0出现了
#3
无效的浮点操作？
n_index:=strtofloat(edit3.Text);
这句设个断点，看看赋值内容如何
n_index:=strtofloat(edit3.Text);
这句设个断点，看看赋值内容如何
#4
我也觉得这里要下断点
或者楼主试试 加默认值
strtofloatdef(edit3.Text,1);
或者楼主试试 加默认值
strtofloatdef(edit3.Text,1);
#5
问题出在你的power()函数，请仔细看看它的参数的数据类型。
#6
n_index:double;
这句 类型申明 改成Real类型 看看
这句 类型申明 改成Real类型 看看
#7
断点跟踪，法宝