У меня нет "Математики", но я написал численный интегратор и получил неожиданные эффекты. Самый яркий из них - ненулевая постоянная составляющая тока рамки.
Интегратор:
// gcc da_float.c -o da_float.cgi -Wall -Werror -O3 -lm
// ./da_float.cgi
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <unistd.h>
#include <math.h>
struct s_tab{ // значения для одного момента времени
double i; // ток
double a; // угол
double v; // скорость
double m; // момент
};
struct s_scale{ // структура автоматического масштабирования
double k1; // ток
double k2; // угол
double k3; // скорость
double k4; // момент
};
int wr_bmp_24_RGB(FILE *fp, const char *pic, int w, int h);
//----------- sign_cos -----------
int sign_cos(double a)
{
double cos_a = cos(a*3.14/180);
if(cos_a > 0) return(1);
return(-1);
}
// масштаб времени (целые положительные числа)
#define T1 8 // T1>=1. Растяжка времени. На T1 делятся приращения.
#define T2 10 // сжатие времени на рисунке, на точность не влияет
// Начальные значения
#define I0 0 // Ток
#define A0 0 // Угол
#define V0 0 // Скорость
#define FI 0 // угол установки шеток. (0 или +90)
// Параметры
#define E0 2 // напряжение батарейки
#define E1 20 // ЭДС-индукции
#define R0 10 // сопротивление потерь
#define M1 0 // момент вязкого трения M=M1*v;
#define M2 0 // момент сухого трения
#define J 9000 // момент инерции, J>0, он в знаменателе
// Размеры рисунка. Можно менять, пошире сделать - больше видно
#define W 600 // ширина
#define H 200 // высота
// Масштабные коэффициенты для графиков
#define AVT 1 // 0 ручная; 1 автоматическая установка кроме T1 и T2
#define K1 0.2 // ток
#define K2 0.9 // угол
#define K3 2 // скорость
#define K4 3 // момент
//----------- avit ---------------
// Один шаг интегрирования
int avit(struct s_tab *s, int t)
{
double di, dv, cos_a, m, m1, m2, v1, z;
double i = s[t].i; // ток
double a = s[t].a; // угол
double v = s[t].v; // скорость вращения
cos_a = cos(a*3.14/180);
z = sign_cos(a+FI); // FI -угол установки щеток. Знак коммутатора.
// +1 батарейка согласно направлению обхода рамки;
// -1 встречно.
// приращение тока. E=L(di/dt); i=интеграл((E/L)dt); di=E/L;
// dt здесь всегда равно 1.
di = E0 * z - E1 * v * cos_a * (3.14/180) - R0*i/1000;
// T1 -растяжка времени, уменьшение дифференциалов
i = i + di/T1; // ток
m = i * cos_a; // момент электромагнитный
m1 = M1 * v; // момент вязкого трения
if(v<0) m2 = -M2; else m2=M2; // m2 -момент сухого трения
dv = (m - m1 - m2)/J; // приращение скорости
v1 = v; // запомнил скорость предыдущего шага
v = v + dv/T1; // скорость
if(M2>0 && abs(m)<=M2){ // сухое трение
if((v1<=0 && v>0) || (v1>=0 && v<0)) v=0;
}
a = a + v/T1; // угол
if(a > 360) a=a-360; // Это для ненакопления ошибки Pi, и чтобы
if(a < -360) a=a+360; // график угла не уходил за край рисунка.
t = t + 1;
s[t].i = i;
s[t].a = a;
s[t].v = v;
s[t].m = m;
return(0);
}
//------------ tab -------------------
// Интегратор, заполняет таблицу
int tab(struct s_tab *s, int w)
{
int t=0, k=w*T2;
s[t].i = I0;
s[t].a = A0;
s[t].v = V0;
s[t].m = 0;
for(t=0; t<k; t++){
avit(s, t);
}
return(0);
}
//---------- avt_scale -------------
// автоматическое масштабирование
int avt_scale(struct s_tab *s, struct s_scale *s1, int h1, int h2, int w)
{
double imin, imax, amin, amax, vmin, vmax, mmin, mmax, k1, k2, k3, k4;
int t, n;
imin=imax=amin=amax=vmin=vmax=mmin=mmax=0;
n = w*T2; h1 += 5; h2 -= 5;
for(t=0; t<n; t++){
if(s[t].i < imin) imin = s[t].i;
if(s[t].i > imax) imax = s[t].i;
if(s[t].a < amin) amin = s[t].a;
if(s[t].a > amax) amax = s[t].a;
if(s[t].v < vmin) vmin = s[t].v;
if(s[t].v > vmax) vmax = s[t].v;
if(s[t].m < mmin) mmin = s[t].m;
if(s[t].m > mmax) mmax = s[t].m;
}
if((-imin) > imax && (-imin) > 0.1) k1 = -h2/imin;
else if(imax > 0.1) k1 = -h1/imax; else k1 = 1;
if((-amin) > amax && (-amin) > 0.1) k2 = -h2/amin;
else if(amax > 0.1) k2 = -h1/amax; else k2 = 1;
if((-vmin) > vmax && (-vmin) > 0.1) k3 = -h2/vmin;
else if(vmax > 0.1) k3 = -h1/vmax; else k3 = 1;
if((-mmin) > mmax && (-mmin) > 0.1) k4 = -h2/mmin;
else if(mmax > 0.1) k4 = -h1/mmax; else k4 = 1;
k4 = k4*0.7; // дизайн рисунка
if(AVT==1) s1->k1 = k1; else s1->k1 = K1;
if(AVT==1) s1->k2 = k2; else s1->k2 = K2;
if(AVT==1) s1->k3 = k3; else s1->k3 = K3;
if(AVT==1) s1->k4 = k4; else s1->k4 = K4;
printf("t=%d\n", w*T2/T1);
printf("Imax=%0.3f; Imin=%0.3f\n", imax, imin);
printf("Vmax=%0.3f; Vmin=%0.3f\n", vmax, vmin);
printf("Mmax=%0.3f; Mmin=%0.3f\n", mmax, mmin);
return(0);
}
//--------- set_pixel -------------------
void set_pixel(char *p, int r, int g, int b)
{
*p++ = r; *p++ = g; *p++ = b;
return;
}
//------------ set_img ----------------
// Преобразование таблицы в рисунок
int set_img(struct s_tab *s, char *pic, int w, int h, FILE *fp)
{
int w3, h1, h2, t, t1, val1, val2, val3, val4;
double k1, k2, k3, k4;
char *p1, *p2;
struct s_scale s1[1];
w3 = w*3;
h2 = (h-2)/2;
h1 = -h2;
p1 = pic + (-h1*w3);
avt_scale(s, s1, h1, h2, w);
k1 = s1->k1;
k2 = s1->k2;
k3 = s1->k3;
k4 = s1->k4;
printf("k1=%0.3f; k2=%0.3f; k3=%0.3f; k4=%0.3f\n", k1, k2, k3, k4);
memset(pic, 255, w3*h);
for(p2=p1, t=0; t<w; t++){
*p2++ = 0; *p2++ = 0; *p2++ = 0;
}
for(t=0; t<w; t++){
for(t1=0; t1<T2; t1++){ // сжатие по t; все в один столбик пикселей.
val1 = (int)(s[t*T2+t1].i * k1);
val2 = (int)(s[t*T2+t1].a * k2);
val3 = (int)(s[t*T2+t1].v * k3);
val4 = (int)(s[t*T2+t1].m * k4);
if(1 && val1 > h1 && val1 < h2){
p2 = p1 - val1*w3;
set_pixel(p2, 255, 0, 0);
}
if(1 && val2 > h1 && val2 < h2){
p2 = p1 - val2*w3;
set_pixel(p2, 0, 255, 0);
}
if(1 && val3 > h1 && val3 < h2){
p2 = p1 - val3*w3;
set_pixel(p2, 0, 0, 255);
}
if(1 && val4 > h1 && val4 < h2){
p2 = p1 - val4*w3;
set_pixel(p2, 180, 100, 0);
}
}
p1 += 3;
}
k1 = wr_bmp_24_RGB(fp, pic, w, h);
return(k1);
}
//------------ set -------------------
int set(FILE *fp)
{
int w=W, h=H;
struct s_tab *s=NULL;
char *pic=NULL;
s = (struct s_tab*)malloc(w * T2 * sizeof(struct s_tab));
if(s==NULL){ printf("s = (struct s_tab*)malloc()=NULL\n"); exit(0);}
pic = (char*)malloc(w*h*3);
if(pic==NULL){ printf("pic = (char*)malloc()=NULL\n"); exit(0);}
tab(s, w);
set_img(s, pic, w, h, fp);
exit(0);
}
//--------- main --------------------
int main()
{
const char *f="out.bmp";
FILE *fp=NULL;
fp=fopen(f, "w");
if(fp==NULL){ printf("fp=fopen(%s)=NULL\n", f); exit(0);}
set(fp);
exit(0);
}
//------------------- wr_bmp_24_RGB ------------------
// записать RGB-буфер в BMP-файл
int wr_bmp_24_RGB(FILE *fp, const char *pic, int w, int h)
{
char hdr[54];
int i, j, dl, pad, w1;
const char *p1;
w1=3*w;
if(w1%4) pad = 4 - w1%4;
else pad=0;
dl = (w1 + pad)*h + 54;
memset(hdr, 0, sizeof(hdr));
hdr[0]='B'; hdr[1]='M';
hdr[2]=dl; hdr[3]=(dl>>8); hdr[4]=(dl>>16); hdr[5]=(dl>>24);
hdr[10]=54; hdr[14]=40;
hdr[18]=w; hdr[19]=(w>>8); hdr[20]=(w>>16); hdr[21]=(w>>24);
hdr[22]=h; hdr[23]=(h>>8); hdr[24]=(h>>16); hdr[25]=(h>>24);
hdr[26]=1; hdr[28]=24;
if(fwrite(hdr, 54, 1, fp) != 1) return(-1);
p1=pic + w1*h;
for(i=0; i<h; i++){
p1 -= w1;
for(j=0; j<w; j++){
if(fputc(*(p1+2), fp) < 0) return(-1);
if(fputc(*(p1+1), fp) < 0) return(-1);
if(fputc(*(p1), fp) < 0) return(-1);
p1 += 3;
}
p1 -= w1;
if(pad>0) if(fwrite("\0\0\0", pad, 1, fp) != 1) return(-1);
}
return(0);
}
Привожу некоторые графики. Вот настройки, они относятся к первому рисунку. Т. к. настроек много, то эти настройки я назвал исходными, а к остальным рисункам я буду приводить лишь те настройки, которые отличаются от исходных или облегчают смотрение графиков, для удобства. Снизу через пунктирную черту привожу выдачу программы, это для сравнения рисунков.
Рис. 1.
Рамка с сопротивлением потерь, с не повернутыми щетками, с ЭДС индукции, без начальной скорости, без нагрузки.
Код:
// масштаб времени (целые положительные числа)
#define T1 8 // T1>=1. Растяжка времени. На T1 делятся приращения.
#define T2 10 // сжатие времени на рисунке, на точность не влияет
// Начальные значения
#define I0 0 // Ток
#define A0 0 // Угол
#define V0 0 // Скорость
#define FI 0 // угол установки шеток. (0 или +90)
// Параметры
#define E0 2 // напряжение батарейки
#define E1 20 // ЭДС-индукции
#define R0 10 // сопротивление потерь
#define M1 0 // момент вязкого трения M=M1*v;
#define M2 0 // момент сухого трения
#define J 9000 // момент инерции, J>0, он в знаменателе
// Размеры рисунка. Можно менять, пошире сделать - больше видно
#define W 600 // ширина
#define H 200 // высота
// Масштабные коэффициенты для графиков
#define AVT 1 // 0 ручная; 1 автоматическая установка кроме T1 и T2
#define K1 0.2 // ток
#define K2 0.9 // угол
#define K3 2 // скорость
#define K4 3 // момент
--------------
t=750
Imax=147.619; Imin=-103.816
Vmax=1.896; Vmin=0.000
Mmax=105.349; Mmin=-35.976
k1=0.637; k2=0.261; k3=49.586; k4=0.625
В выдаче время всего графика, значения тока, скорости и момента, максимальные на рисунке, а также масштабные коэффициенты в порядке: ток, угол, скорость, момент.
На рисунке:
Красный - ток рамки.
Зеленый - угол рамки, когда 360 градусов - скачок вниз.
Синий - скорость рамки.
Оранжевый - электромагнитный вращательный момент.
Рис. 2.
Тоже, что на рис. 1, но сильно сжатое по времени.
Код:
#define T2 10000 // сжатие времени на рисунке, на точность не влияет
-------------
t=750000
Imax=147.619; Imin=-103.816
Vmax=7.302; Vmin=0.000
Mmax=105.349; Mmin=-35.976
k1=0.637; k2=0.261; k3=12.873; k4=0.625
Видно установившуюся скорость, она максимальная на рисунке. Из выдачи:
Vmax=7.302, запомню её, это скорость холостого хода.
Рис. 3.
Тоже, что на рис. 1, но с другим масштабом времени и с начальной скоростью, превышающей скорость холостого хода
Код:
#define T2 1 // сжатие времени на рисунке, на точность не влияет
#define V0 8 // Скорость
-------
t=75
Imax=3.224; Imin=-2.849
Vmax=8.000; Vmin=0.000
Mmax=2.204; Mmin=-2.547
k1=29.153; k2=0.261; k3=11.750; k4=25.830
Рис. 4.
То же, что на предыдущем рисунке, но сильно сжатое по времени
Код:
#define T2 10000 // сжатие времени на рисунке, на точность не влияет
#define V0 8 // Скорость
--------
t=750000
Imax=5.193; Imin=-5.503
Vmax=8.000; Vmin=0.000
Mmax=2.679; Mmin=-2.810
k1=17.082; k2=0.261; k3=11.750; k4=23.419
Видно, что скорость снижается примерно до скорости холостого хода и такой и остается. Когда она снижается, наверно батарейка подзаряжается макроскопически, т. е. в среднем за много оборотов.
Рис. 5.
То же, что на рис. 1, но с другим масштабом времени и с нагрузкой
Код:
#define T2 100 // сжатие времени на рисунке, на точность не влияет
#define M1 2 // момент вязкого трения M=M1*v;
-----------
t=7500
Imax=147.914; Imin=-106.218
Vmax=2.407; Vmin=0.000
Mmax=105.485; Mmin=-35.573
k1=0.636; k2=0.261; k3=39.058; k4=0.624
Видно, что скорость установилась
Vmax=2.407, и она сильно меньше скорости холостого хода.
Рис. 6.
Сильно отличается от всех прежних, теперь рамка без потерь, и щетки сдвинуты. Здесь самое важное, что ток всё время нарастает макроскопически, т. е. появляется постоянная составляющая.
Код:
// Начальные значения
#define I0 0 // Ток
#define A0 -45 // Угол
#define V0 0 // Скорость
#define FI 90 // угол установки щеток. (0 или +90)
#define T2 16 // сжатие времени на рисунке, на точность не влияет
#define R0 0 // сопротивление потерь
#define M1 0 // момент вязкого трения M=M1*v;
---------
t=1200
Imax=326.616; Imin=0.000
Vmax=4.752; Vmin=0.000
Mmax=326.615; Mmin=-241.398
k1=0.288; k2=0.261; k3=19.782; k4=0.201
Вот ещё два рисунка, на них то же самое, что на предыдущем, меняется только масштаб времени.
Рис. 7.
Код:
#define T2 200 // сжатие времени на рисунке, на точность не влияет
#define A0 -45 // Угол
#define FI 90 // угол установки щеток. (0 или +90)
#define R0 0 // сопротивление потерь
-----------
t=15000
Imax=1117.377; Imin=0.000
Vmax=15.451; Vmin=0.000
Mmax=1106.786; Mmin=-1084.272
k1=0.084; k2=0.261; k3=6.084; k4=0.059
Рис. 8.
Код:
#define T2 20000 // сжатие времени на рисунке, на точность не влияет
#define A0 -45 // Угол
#define FI 90 // угол установки щеток. (0 или +90)
#define R0 0 // сопротивление потерь
---------
t=1500000
Imax=30673.043; Imin=-48121.550
Vmax=610.793; Vmin=-559.168
Mmax=48093.708; Mmin=-48108.711
k1=0.002; k2=0.261; k3=0.154; k4=0.001
Последние два рисунка выглядят нехорошо. Предполагаю, что или рамка без потерь на самом деле макроскопически неустойчива, или это из-за накопления погрешностей. Рисунки сомнительны, и их можно отбраковать. Но рис. 6, там вроде бы всё в порядке. Мне удавалось получить подобный эффект - постоянную составляющую - и в рамке с потерями, но там всплеск был значительно меньше и плавно сходил на нет. Я поначалу не вёл записи, тыкал наугад, и сейчас мне не воспроизвести.
Предположим, что моё предположение о постоянной составляющей верно. Тогда вот что получается. В рамке течет большой постоянный ток, значительно превышающий размах пилы. Рамка крутится. Когда она повернётся своей положительной стороной к северному полюсу магнита, она сильно тормозит, а когда повернётся к южному полюсу - сильно ускоряется. Возможно, это причина сильных микроскопических колебаний скорости.