API tools faq
paste
Advertisement
Guest User

Gauss Femm calculator

a guest
Feb 7th, 2012
473
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
Lua 20.32 KB | None | 0 0
raw download clone embed print report
  1. -- Это lua скрипт для FEMM 4.0
  2.  
  3. -- Функция вычисления (переделана из gauss_ver_2_4_04)
  4. function calculate(C,U,Dpr,Tiz,Lkat,Dkat,Lmag,Nom_mat_magnitoprovoda,Nom_mat_puli,Lpuli,Dpuli,Lsdv,Dstvola,Vel0,delta_t)
  5.     params=C..","..U..","..Dpr..","..Tiz
  6.     params=params..","..Lkat..","..Dkat..","..Lmag..","..Nom_mat_magnitoprovoda
  7.     params=params..","..Nom_mat_puli..","..Lpuli..","..Dpuli..","..Lsdv..","..Dstvola..","..Vel0..","..delta_t
  8.     kRC = 140      -- Постоянная константа RC для распространённых электрколитических нденсаторов, Ом*мкФ
  9.     Rcc = (kRC/C)  -- Внутреннее сопротивление конденсатора
  10.     Rv = 0.35+Rcc     -- Сопротивление подводящих проводов + сопротивление тиристора + внутреннее сопротивление конденсатора, Ом
  11.     Vol = 3                 -- Кратность свободного пространства вокруг модели (рекомендуется значение от 3 до 5)
  12.     Coil_meshsize = 0.5    -- Размер сетки катушки, мм
  13.     Proj_meshsize = 0.35    -- Размер сетки пули, мм
  14.     max_segm      = 5     -- Максимальный размер сегмента пространства, град
  15.     t = 0           -- Начальный момент времени, секунды
  16.     sigma = 0.0000000175    -- Удельное сопротивление меди, Ом * Метр
  17.     ro = 7800       -- Плотность железа, Кг/Метр^3
  18.     pi = 3.1415926535  
  19.     --------------------------------------------------------------------------------
  20.     -- Начинаем
  21.     --------------------------------------------------------------------------------
  22.     Start_date= date()
  23.     create(0)                           -- создаем документ для магнитных задач
  24.     mi_probdef(0,"millimeters","axi",1E-8,30)           -- создаем задачу
  25.     mi_saveas("temp.fem")                       -- сохраняем файл под другим именем
  26.     mi_addmaterial("Air",1,1)                   -- добавляем материал воздух
  27.     mi_addmaterial("Cu",1,1,"","","",58,"","","",4,"","",1,Dpr) -- добавляем материал медный провод диаметром Dpr проводимостью 58
  28.     mi_addcircprop("katushka",0,0,1)                -- добавляем катушку
  29.     dofile ("func.lua") -- компилируем функции для дальнейшего доступа к ним                
  30.     vvod_materiala (Nom_mat_magnitoprovoda,"M") -- вводим материал магнитопровода назовем его как М и номер материала
  31.     Material_magnitoprovoda="M" .. Nom_mat_magnitoprovoda
  32.     vvod_materiala (Nom_mat_puli,"P")       ---- вводим материал пули назовем его как Р и номер материала
  33.     Material_puli="P" .. Nom_mat_puli  
  34.     --------------------------------------------------------------------------------
  35.     -- Располагаем объекты
  36.     --------------------------------------------------------------------------------
  37.     --Создаем пространство в Vol раз большее чем катушка
  38.     mi_addnode(0,(Lkat+Lmag)*-Vol)              -- ставим точку
  39.     mi_addnode(0,(Lkat+Lmag)*Vol)               -- ставим точку
  40.     mi_addsegment(0,(Lkat+Lmag)*-Vol,0,(Lkat+Lmag)*Vol)     -- рисуем линию
  41.     mi_addarc(0,(Lkat+Lmag)*-5,0,(Lkat+Lmag)*Vol,180,max_segm)  -- рисуем дугу
  42.     mi_addblocklabel((Lkat+Lmag)*0.7*Vol,0)             -- добавляем блок 
  43.     mi_clearselected()                      -- отменяем все
  44.     mi_selectlabel((Lkat+Lmag)*0.7*Vol,0)               -- выделяем метку блока
  45.     mi_setblockprop("Air", 1, "", "", "",0)             -- устанавливаем свойства блока с имнем Air и номером блока 0
  46.     mi_zoomnatural()    -- устанавливаем зум так что бы было видно на весь экран
  47.     -------------------------------------------------------------------------- Создаем пулю
  48.     if Dstvola < Dpuli then Dstvola = Dpuli+0.1 end -- защита от неумех
  49.     -- если длина пули равна диаметру значит шар
  50.     if Lpuli==Dpuli then
  51.         mi_addnode(0,Lkat/2-Lsdv)
  52.         mi_addnode(0,Lkat/2+Lpuli-Lsdv)
  53.         mi_clearselected()
  54.         mi_selectnode (0,Lkat/2-Lsdv)
  55.         mi_selectnode (0,Lkat/2+Lpuli-Lsdv)
  56.         mi_setnodeprop("",1)
  57.         mi_addarc(0,Lkat/2-Lsdv,0,Lkat/2+Lpuli-Lsdv,180,5)
  58.     else    -- иначе просто цилиндр
  59.         mi_addnode(0,Lkat/2-Lsdv)
  60.         mi_addnode(Dpuli/2,Lkat/2-Lsdv)
  61.         mi_addnode(Dpuli/2,Lkat/2+Lpuli-Lsdv)
  62.         mi_addnode(0,Lkat/2+Lpuli-Lsdv)
  63.         mi_clearselected()
  64.         mi_selectnode(0,Lkat/2-Lsdv)
  65.         mi_selectnode(Dpuli/2,Lkat/2-Lsdv)
  66.         mi_selectnode(Dpuli/2,Lkat/2+Lpuli-Lsdv)
  67.         mi_selectnode(0,Lkat/2+Lpuli-Lsdv)
  68.         mi_setnodeprop("",1)
  69.         mi_addsegment(Dpuli/2,Lkat/2-Lsdv,Dpuli/2,Lkat/2+Lpuli-Lsdv)
  70.         mi_addsegment(Dpuli/2,Lkat/2+Lpuli-Lsdv,0,Lkat/2+Lpuli-Lsdv)
  71.         mi_addsegment(0,Lkat/2+Lpuli-Lsdv,0,Lkat/2-Lsdv)
  72.         mi_addsegment(0,Lkat/2-Lsdv,Dpuli/2,Lkat/2-Lsdv)
  73.     end
  74.     mi_addblocklabel(Dpuli/4,Lkat/2+Lpuli/2-Lsdv)
  75.     mi_clearselected()
  76.     mi_selectlabel(Dpuli/4,Lkat/2+Lpuli/2-Lsdv)
  77.     mi_setblockprop(Material_puli, 0, Proj_meshsize, "", "",1)          -- номер блока 1
  78.     ------------------------------------------------------------------------- Создаем катушку
  79.     mi_addnode(Dstvola/2,Lkat/2)            -- основание
  80.     mi_addnode(Dstvola/2,-Lkat/2)           -- основание
  81.     mi_addnode(Dkat/2,Lkat/2)               -- внешняя начальная часть
  82.     mi_addnode(Dkat/2,-Lkat/2)              -- внешняя конечная часть
  83.     mi_addsegment(Dstvola/2,-Lkat/2,Dstvola/2,Lkat/2)
  84.     mi_addsegment(Dstvola/2,Lkat/2,Dkat/2,Lkat/2)
  85.     mi_addsegment(Dkat/2,Lkat/2,Dkat/2,-Lkat/2)
  86.     mi_addsegment(Dkat/2,-Lkat/2,Dstvola/2,-Lkat/2)
  87.     mi_clearselected()
  88.     mi_selectnode(Dstvola/2,Lkat/2)         -- основание
  89.     mi_selectnode(Dstvola/2,-Lkat/2)        -- основание
  90.     mi_selectnode(Dkat/2,Lkat/2)               
  91.     mi_selectnode(Dkat/2,-Lkat/2)              
  92.     mi_setnodeprop("",2)
  93.     mi_addblocklabel(Dstvola/2+(Dkat/2-Dstvola/2)/2,0)
  94.     mi_clearselected()
  95.     mi_selectlabel(Dstvola/2+(Dkat/2-Dstvola/2)/2,0)
  96.     mi_setblockprop("Cu", 0, Coil_meshsize, "katushka", "",2) -- номер блока 2
  97.     -------------------------------------------------------------------------- Создаем внешний магнитопровод
  98.     if (Lmag > 0) and (Nom_mat_magnitoprovoda > 0) then
  99.         mi_addnode(Dstvola/2,Lkat/2+Lmag)
  100.         mi_addnode(Dkat/2+Lmag,Lkat/2+Lmag)
  101.         mi_addnode(Dkat/2+Lmag,-Lkat/2-Lmag)   
  102.         mi_addnode(Dstvola/2,-Lkat/2-Lmag)
  103.         mi_addsegment(Dstvola/2,Lkat/2,Dstvola/2,Lkat/2+Lmag)
  104.         mi_addsegment(Dstvola/2,Lkat/2+Lmag,Dkat/2+Lmag,Lkat/2+Lmag)
  105.         mi_addsegment(Dkat/2+Lmag,Lkat/2+Lmag,Dkat/2+Lmag,-Lkat/2-Lmag)
  106.         mi_addsegment(Dkat/2+Lmag,-Lkat/2-Lmag,Dstvola/2,-Lkat/2-Lmag)
  107.         mi_addsegment(Dstvola/2,-Lkat/2-Lmag,Dstvola/2,-Lkat/2)
  108.         mi_addblocklabel(Dkat/2+Lmag/2,0)
  109.         mi_clearselected()
  110.         mi_selectlabel(Dkat/2+Lmag/2,0)
  111.         mi_setblockprop(Material_magnitoprovoda, 1, "", "", "",3)       -- номер блока 3
  112.     end
  113.     mi_clearselected()
  114.     --------------------------------------------------------------------------------
  115.     -- система СИ - метр, Фарад
  116.     --------------------------------------------------------------------------------
  117.     C = C/1000000
  118.     Dpriz = Dpr+Tiz -- Диаметр провода в изоляции
  119.     Dpr = Dpr/1000
  120.     Dpriz = Dpriz/1000     
  121.     Lpuli  = Lpuli/1000
  122.     Dpuli = Dpuli/1000
  123.     Dstvola = Dstvola/1000             
  124.     Lkat = Lkat/1000
  125.     Dkat = Dkat/1000
  126.     Lsdv = Lsdv/1000
  127.     Lmag = Lmag/1000
  128.     --------------------------------------------------------------------------------
  129.     -- Анализируем и запускаем постпроцессор
  130.     mi_analyze(1)               -- анализируем (скрывая окно анализа "1") 0 - будет видно окно и будет работать медленее
  131.     mi_loadsolution()           -- запускаем саму программу пост процессора
  132.     mo_groupselectblock(2)
  133.     Skat = mo_blockintegral(5)      -- Площадь сечения катушки, Метр^2
  134.     Vkat = mo_blockintegral(10)     -- Объем катушки, Метр^3
  135.     mo_clearblock()
  136.     mo_groupselectblock(1)
  137.     Vpuli = mo_blockintegral(10)        -- Объем пули, Метр^3 
  138.     mo_clearblock()            
  139.     Mpuli=ro*Vpuli              -- Масса пули, кг
  140.     N=Skat*0.94/(Dpriz*Dpriz)       -- Количество витков в катушке уточнённое
  141.     DLprovoda=N * 2 * pi * (Dkat + Dstvola)/4   -- Длина обмоточного провода уточнённая, м
  142.     Rkat=sigma*DLprovoda/(pi*(Dpr/2)^2) -- Сопротивление всего обмоточного провода катушки, Ом
  143.     R=Rv+Rkat               -- Полное сопротивление системы
  144.     --Устанавливаем число витков, а силу тока 100 А для оценки индуктивности
  145.     mi_clearselected()
  146.     mi_selectlabel(Dstvola*1000/2+(Dkat/2-Dstvola/2)*1000/2,0)
  147.     mi_setblockprop("Cu", 0, Coil_meshsize, "katushka", "",2,N) -- последнее значение - число витков
  148.     mi_clearselected()
  149.     mi_modifycircprop("katushka",1,100)
  150.     -- Анализируем и запускаем постпроцессор
  151.     mi_analyze(1)               -- анализируем (скрывая окно анализа "1")  
  152.     mo_reload()             -- перезапускаем программу пост процессора         
  153.     current_re,current_im,volts_re,volts_im,flux_re,flux_im=mo_getcircuitproperties("katushka") -- получаем данные с катушки
  154.     L=flux_re/current_re            -- Оценочная индуктивность, Генри
  155.     --------------------------------------------------------------------------------
  156.     -- НАчало симуляции
  157.     --------------------------------------------------------------------------------
  158.     dt = delta_t/1000000 -- перевод приращения времени в секунды
  159.     x=0     -- начальная позиция пули
  160.     I0=0.00000001   -- достаточно малое значение тока
  161.     t=0     -- общее время
  162.     Vel=Vel0
  163.     Vmax=Vel
  164.     Uc = U
  165.     I=I0        -- начальное значение тока
  166.     Force = 0
  167.     Fii = 0
  168.     Fix = 0
  169.     KC=1        -- счетчик циклов, для вывода в файл
  170.     T_I={}      -- создаем массив (таблицу как её называют в Lua)
  171.     T_F={}     
  172.     T_Vel={}   
  173.     T_x={}     
  174.     T_t={}
  175.    
  176.     maxI=0
  177.     Isum=0
  178.     Ic=0
  179.    
  180.     repeat      ------------------------------------------------------------ начинаем цикл
  181.         t = t+dt
  182.         --- Рассчитываем dFi/dI при I и силу
  183.                 mi_modifycircprop("katushka",1,I)   -- Устанавливает ток
  184.                 mi_analyze(1)           -- анализируем (скрывая окно анализа "1") 0 - будет видно окно и будет работать медленее  
  185.                 mo_reload()             -- перезапускаем программу пост процессора
  186.                 mo_groupselectblock(1)
  187.         Force = mo_blockintegral(19)        -- Сила действующая на пулю, Ньютон 
  188.         Force=Force*-1              -- ставим "-" из за координат (направление силы в сторону уменьшения координаты)
  189.         current_re,current_im,volts_re,volts_im,flux_re,flux_im=mo_getcircuitproperties("katushka") -- получаем данные с катушки
  190.         Fi0=flux_re                     -- магнитный поток
  191.             mi_modifycircprop("katushka",1,I*1.001) -- Устанавливает ток, увельченный на 1.001
  192.             mi_analyze(1)               -- анализируем (скрывая окно анализа "1") 0 - будет видно окно и будет работать медленее  
  193.             mo_reload()             -- перезапускаем программу пост процессора         
  194.         current_re,current_im,volts_re,volts_im,flux_re,flux_im=mo_getcircuitproperties("katushka") -- получаем данные с катушки
  195.         Fi1=flux_re                     -- магнитный поток при I=I+0.001*I, dI=0.001*I
  196.         Fii=(Fi1-Fi0)/(0.001*I)                              -- Рассчитываем dFi/dI
  197.         apuli = Force / Mpuli           -- Ускорение пули, Метр/секунда^2
  198.         dx = Vel*dt+apuli*dt*dt/2       -- Приращение координаты, метр
  199.         x = x+dx                -- Новая позиция пули
  200.         Vel = Vel+apuli*dt          -- Скорость после приращения, метр/секунда
  201.         if tonumber(Vmax)<tonumber(Vel) then Vmax=tonumber(Vel) end
  202.         mi_selectgroup(1)           -- Выделяем пулю
  203.         mi_movetranslate(0,-dx*1000)        -- Перемещаем её на dx
  204.         --- Рассчитываем dFi/dx при x
  205.         Fix= Force/I
  206.         ------- Расчитываем ток и напряжение на конденсаторе
  207.         I=I+dt*(Uc-I*R-Fix*Vel)/Fii            
  208.         Uc = Uc-dt*I/C
  209.         if Uc< 0 then  Uc=0 end   --если стоит паралельный диод
  210.         if x > (Lpuli+Lkat-Lsdv) and Force >-1  then I=0 end -- Пуля вылетела за пределы катушки и сила очень маленькая
  211.         if x < 0 then I=0 end -- Пуля вылетела назад  
  212.         T_I[KC]=I       -- записываем данные в массив
  213.         T_F[KC]=Force
  214.         T_Vel[KC]=Vel      
  215.         T_x[KC]=x*1000     
  216.         T_t[KC]=t*1000000  
  217.         KC=KC+1
  218.        
  219.         if I>maxI then
  220.             maxI=I
  221.         end
  222.         Isum=Isum+I
  223.         Ic=Ic+1
  224.     until (I <= 0) or (Vel < 0 ) -- повторяем расчет, пока не будет тока
  225.     Isr=Isum/Ic
  226.  
  227.     Epuli = (Mpuli*Vel^2)/2 - (Mpuli*Vel0^2)/2
  228.     EC= (C*U^2)/2
  229.     KPD = Epuli*100/EC
  230.  
  231.     print ("Симуляция завершена. Данные записаны в файл "..File_name.." result.csv")
  232.     ----------------------------------------------------------------------------------------------------
  233.     -- Записываем всё в csv файл
  234.     ----------------------------------------------------------------------------------------------------
  235.     handle = openfile(File_name .. " result.csv", "a")-- создаем файл а - будем дописывать в конец файла w - записать стерев всё что было перед тем
  236.     res=params ..",".. Start_date..","..R..",".. Rv..","..Rkat..","
  237.     res=   res..","..N..","..tostring(Dpr*1000)..","..DLprovoda..","..tostring(Lkat*1000)
  238.     res=   res..","..tostring(Dkat*1000)..","..tostring(L*1000000)..","..tostring(Mpuli*1000)..",".. tostring(t*1000000)
  239.     res=   res..","..Epuli..",".. EC..",".. KPD..",".. Vel..",".. Vmax..","..maxI..","..Isr..",".. date().."\n"
  240.     write(handle,res)
  241.     closefile(handle)
  242.  
  243.     -- Удаляем промежуточные файлы
  244.     remove ("temp.fem")
  245.     remove ("temp.ans")
  246. end
  247.  
  248. -- Чтение блока допустимых значений параметра
  249. function readBlock(handle)
  250.     result={}
  251.     local stop=0
  252.     i=1
  253.     while stop==0 do
  254.         line = read(handle, "*l")
  255.         if line=="[/]" then stop=1
  256.                        else result[i]=tonumber(line)
  257.         end
  258.         i=i+1
  259.     end
  260.     return result
  261. end
  262.  
  263. function calculateParams(params)
  264.     calculate(params[1 ],params[2 ],params[3 ],
  265.               params[4 ],params[5 ],params[6 ],
  266.               params[7 ],params[8 ],params[9 ],
  267.               params[10],params[11],params[12],
  268.               params[13],params[14],params[15])
  269. end
  270.  
  271. -- Подбор параметров (рекуррентный)
  272. function change(params,num,other)
  273.     if num==16 then
  274.         calculateParams(other) -- если все значения параметров подобраны, то считаем
  275.     else
  276.         local param=params[num] -- Допустимые значения меняемого параметра
  277.         for key,value in param do
  278.             local my = other
  279.             my[num]=value
  280.             change(params,num+1,my)
  281.         end
  282.     end
  283. end
  284.  
  285.  
  286. showconsole()                           -- показываем окно вывода промежуточных данных
  287. clearconsole()
  288.  
  289. File_name=prompt ("Введите имя файла с данными гауса, без расширения .txt")
  290.  
  291. print("Чтение данных...")
  292. handle = openfile(File_name .. ".txt","r")
  293.  
  294. params={}
  295. --  1. Ёмкость конденсаторов
  296. --  2. Напряжение конденсаторов
  297. --  3. Диаметр обмоточного провода катушки, милиметр
  298. --  4. Удвоенная толщина изоляции провода (разница диаметра в изоляции и диаметра голого), мм
  299. --  5. Длина катушки (не задавать меньше диаметра обмот. провода катушки), милиметр             
  300. --  6. Внешний диаметр катушки, милиметр
  301. --  7. Толщина внешнего магнитопровода, по форме повторяет катушку, если ноль то его нет, милиметр
  302. --  8. материал из которого сделан магнитопровод катушки см. таблицу
  303. --  9. материал из которого сделана пуля см. таблицу
  304. -- 10. Длина пули, милиметр       
  305. -- 11. Диаметр пули, милиметр
  306. -- 12. Расстояние, на которое в начальный момент вдвинута пуля в катушку или находится до катушки с минусом, милиметр
  307. -- 13. Внешний диаметр ствола (не задавать меньше диаметра пули), милиметр
  308. -- 14. Начальная скорость пули, м/с (Вместо 0 лучше какое-то небольшое значение, иначе долго на месте стоит)
  309. -- 15. Приращение времени, мкС
  310.  
  311. comb_count=1
  312. for i=1,13 do
  313.     params[i]=readBlock(handle)
  314. end
  315. params[14]={read(handle, "*l")}
  316. params[15]={read(handle, "*l")}
  317.  
  318. closefile(handle)
  319.  
  320. handle=openfile(File_name .. " result.csv", "w")
  321. write (handle,"Ёмкость батареи Ф,")
  322. write (handle,"Напряжение батареи В,")
  323. write (handle,"Диаметр провода мм,")
  324. write (handle,"Удвоенная толщина изоляции мм,")
  325. write (handle,"Длина катушки мм,")
  326. write (handle,"Внешний диаметр катушки мм,")
  327. write (handle,"Толщина внешнего магнитопровода мм,")
  328. write (handle,"материал магнитопровода,")
  329. write (handle,"материал пули,")
  330. write (handle,"Длина пули мм,")
  331. write (handle,"Диаметр пули,")
  332. write (handle,"Расстояние, на которое в начальный момент вдвинута пуля в катушку,")
  333. write (handle,"Внешний диаметр ствола,")
  334. write (handle,"Начальная скорость пули,")
  335. write (handle,"Приращение времени,")
  336. write (handle,"Начало симуляции ,")
  337. write (handle,"Сопротивление общее Ом,")
  338. write (handle,"Внешнее сопротивление Ом,")
  339. write (handle,"Сопротивление катушки Oм,")
  340. write (handle,"Количество витков в катушке,")
  341. write (handle,"Диаметр обмоточного провода катушки милиметр,")
  342. write (handle,"Длина провода в катушке метр,")
  343. write (handle,"Длина катушки милиметр,")
  344. write (handle,"Внешний диаметр катушки милиметр,")
  345. write (handle,"Индуктивность катушки с пулей в начальном положении микроГенри,")
  346. write (handle,"Масса пули грамм,")
  347. write (handle,"Время процесса (микросек),")
  348. write (handle,"Энергия пули Дж,")
  349. write (handle,"Энергия конденсатора Дж,")
  350. write (handle,"КПД гауса(%),")
  351. write (handle,"Скорость пули на выходе из катушки м/с,")
  352. write (handle,"Максимальная скорость которая была достигнута м/с,")
  353. write (handle,"Максимальный ток А,")
  354. write (handle,"Средний ток А,")
  355. write (handle,"Окончания симуляции\n")
  356. closefile(handle)
  357.  
  358. change(params,1,{}) -- Начинаем подбор
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!