2025年8月16日土曜日

NFCアンテナ(6)というかOpenEMS

 

NFCアンテナのシミュレーションをしてみたくなったのでOpenEMSでやってみた。


NFCで主流のISO14443Aってproximity card(近接カード)とかcontactless card(非接触カード)であって、タッチすることが前提の規格だと思うんだけど、なぜか通信距離を求められることが多い。だったらvicinity card(近傍カード=ISO15693)にすればいいのにってのは許されない。いろんなしがらみがあるんだと思う。ほんと大人の世界はいやらしい。は、さておき、こういうのなかなか実験できない場合にシミュレーションで何とかしたくなることもある(実験の方が早いことも多々あるので決断は慎重に)。
では、OpenEMSでシミュレーションしてみよう。OpenEMSの環境構築はまぁあっちこっち検索すればわかるっしょ。Octaveをインストールして、、、云々。
例によって唐突にソースコードはこちら。

gen_small_loop_coil.m
  1. function [CSX port]=gen_small_loop_coil(
  2.   CSX,
  3.   loop_name,
  4.   coil_width,coil_height,coil_nturn,
  5.   metal_width,metal_pitch,metal_thickness,
  6.   sub_thickness,
  7.   series_cap_value,parallel_cap_value,port_n,port_excite,
  8.   x_origin,y_origin,z_origin
  9. )
  10.  
  11. metal_name=[loop_name '_metal'];
  12. CSX=AddMetal(CSX,metal_name);
  13. parallel_cap_name=[loop_name '_pcap_z'];
  14. CSX=AddLumpedElement(CSX,parallel_cap_name,'z','C',parallel_cap_value);
  15. series_cap_name=[loop_name '_scap_z'];
  16. CSX=AddLumpedElement(CSX,series_cap_name,'z','C',series_cap_value);
  17. x_start=coil_width/2-metal_width/2;
  18. x_end=x_start-(coil_nturn)*metal_pitch;
  19. x=x_start:-metal_pitch:x_end;
  20. y_start=coil_height/2-metal_width/2;
  21. y_end=y_start-(coil_nturn)*metal_pitch;
  22. y=y_start:-metal_pitch:y_end;
  23.  
  24. z0=-sub_thickness/2;
  25. z1=sub_thickness/2;
  26.  
  27. start_stop_vector=[metal_width/2 metal_width/2 metal_thickness/2];
  28. origins=[x_origin y_origin z_origin];
  29.  
  30. for i=1:coil_nturn
  31.   if i>1
  32.     point1=[-x(i-1) -y(i) z1]+origins-start_stop_vector;
  33.     point2=[x(i) -y(i) z1]+origins+start_stop_vector;
  34.   else
  35.     point1=[-x(coil_nturn) -y(i) z1]+origins-start_stop_vector;
  36.     point2=[x(i) -y(i) z1]+origins+start_stop_vector;
  37.   endif
  38.   CSX=AddBox(CSX,metal_name,1,point1,point2);
  39.  
  40.   point1=[x(i) -y(i) z1]+origins-start_stop_vector;
  41.   point2=[x(i) y(i) z1]+origins+start_stop_vector;
  42.   CSX=AddBox(CSX,metal_name,1,point1,point2);
  43.  
  44.   point1=[-x(i) y(i) z1]+origins-start_stop_vector;
  45.   point2=[x(i) y(i) z1]+origins+start_stop_vector;
  46.   CSX=AddBox(CSX,metal_name,1,point1,point2);
  47.  
  48.   point1=[-x(i) -y(i+1) z1]+origins-start_stop_vector;
  49.   point2=[-x(i) y(i) z1]+origins+start_stop_vector;
  50.   CSX=AddBox(CSX,metal_name,1,point1,point2);
  51. endfor
  52.  
  53. point1=[-x(coil_nturn) -y(1) z0]+origins-start_stop_vector;
  54. point2=[-x(coil_nturn) -y(coil_nturn+1) z0]+origins+start_stop_vector;
  55. CSX=AddBox(CSX,metal_name,1,point1,point2);
  56.  
  57. point1=[-x(coil_nturn) -y(coil_nturn+1) z0]+origins-start_stop_vector;
  58. point2=[-x(coil_nturn) -y(coil_nturn+1) z1]+origins+start_stop_vector;
  59. CSX=AddBox(CSX,series_cap_name,10,point1,point2);
  60.  
  61. point1=[-x(coil_nturn) -y(1) z0]+origins+start_stop_vector.*[-1 1 1];
  62. point2=[-x(coil_nturn) -y(1) z1]+origins+start_stop_vector.*[1 1 -1];
  63. CSX=AddBox(CSX,parallel_cap_name,10,point1,point2);
  64.  
  65. point1=[-x(coil_nturn) -y(1) z0]+origins+start_stop_vector.*[-1 -1 1];
  66. point2=[-x(coil_nturn) -y(1) z1]+origins+start_stop_vector.*[1 -1 -1];
  67. [CSX port]=AddLumpedPort(CSX,100,port_n,50,point1,point2,[0 0 1],port_excite);
  68.  
  69.  
  70. endfunction

nfc_trx_analysis_main.m
  1. function [freq s11 s21 s12 s22]=nfc_trx_analysis_main(
  2.   port_direction,
  3.   tx_coil_width,tx_coil_height,tx_coil_nturn,
  4.   rx_coil_width,rx_coil_height,rx_coil_nturn,
  5.   trx_distance,
  6.   disgeomplot
  7. )
  8.  
  9. physical_constants;
  10. unit=1e-3;
  11.  
  12. FDTD=InitFDTD('NrTS',4e6,'EndCriteria',1e-3);
  13. f0=14e6;
  14. fc=12e6;
  15. FDTD=SetGaussExcite(FDTD,f0,fc);
  16. BC={'MUR','MUR','MUR','MUR','MUR','MUR'};
  17. FDTD=SetBoundaryCond(FDTD,BC);
  18.  
  19. CSX=InitCSX();
  20. if port_direction>0
  21.   txport_excite=false;
  22.   rxport_excite=true;
  23. else
  24.   txport_excite=true;
  25.   rxport_excite=false;
  26. endif
  27. [CSX port1]=gen_small_loop_coil(
  28.   CSX,
  29.   'tx',
  30.   tx_coil_width,tx_coil_height,tx_coil_nturn,
  31.   0.2,0.4,0,
  32.   0.2,
  33.   1,1e-18,1,txport_excite,
  34.   0,0,trx_distance/2
  35. );
  36. [CSX port2]=gen_small_loop_coil(
  37.   CSX,
  38.   'rx',
  39.   rx_coil_width,rx_coil_height,rx_coil_nturn,
  40.   0.2,0.4,0,
  41.   0.2,
  42.   1,1e-18,2,rxport_excite,
  43.   0,0,-trx_distance/2
  44. );
  45.  
  46. ## magnetic material
  47. ##CSX = AddMaterial(CSX, 'mag_material');
  48. ##CSX = SetMaterialProperty(CSX, 'mag_material', 'Mue',45);
  49. ##x_mag=rx_coil_width/2-(rx_coil_nturn-1)*0.4-0.2-1;
  50. ##y_mag=rx_coil_height/2-(rx_coil_nturn-1)*0.4-0.2-1;
  51. ##z_mag=0.2/2;
  52. ##mag_origins=[0 0 -trx_distance/2];
  53. ##CSX=AddBox(CSX,'mag_material',1,[-x_mag -y_mag z_mag]+mag_origins,[x_mag y_mag z_mag]+mag_origins);
  54. ##
  55.  
  56. air_box=16;
  57. mesh_box_size=[20 20 16];
  58. mesh_box_resolution=16;
  59. mesh=DetectEdges(CSX);
  60. mesh.x=[-mesh_box_size(1) mesh.x +mesh_box_size(1)];
  61. mesh.y=[-mesh_box_size(2) mesh.y +mesh_box_size(2)];
  62. mesh.z=[-mesh_box_size(3) mesh.z +mesh_box_size(3)];
  63.  
  64. mesh=SmoothMesh(mesh,mesh_box_resolution);
  65.  
  66. mesh.x=[-air_box+mesh.x(1) mesh.x mesh.x(end)+air_box];
  67. mesh.y=[-air_box+mesh.y(1) mesh.y mesh.y(end)+air_box];
  68. mesh.z=[-air_box+mesh.z(1) mesh.z mesh.z(end)+air_box];
  69.  
  70. mesh=SmoothMesh(mesh,c0/(f0+fc)/unit/10,1.5,'algorithm',1);
  71.  
  72. mesh=AddPML(mesh,4);
  73.  
  74. disp(['number of cells: ' num2str(1e-6*numel(mesh.x)*numel(mesh.y)*numel(mesh.z)) 'Mcells'])
  75.  
  76. CSX=DefineRectGrid(CSX,unit,mesh);
  77.  
  78. Sim_Path=['tmp_' mfilename];
  79. Sim_CSX=[mfilename '.xml'];
  80. confirm_recursive_rmdir(0);
  81. [status,message,messageid]=rmdir(Sim_Path,'s');
  82. [status,message,messageid]=mkdir(Sim_Path);
  83. WriteOpenEMS([Sim_Path '/' Sim_CSX],FDTD,CSX);
  84. if disgeomplot==0
  85.   CSXGeomPlot( [Sim_Path '/' Sim_CSX],['--export-polydata-vtk=' Sim_Path ' --RenderDiscMaterial -v']);
  86. endif
  87. RunOpenEMS(Sim_Path,Sim_CSX);
  88.  
  89. port={port1 port2};
  90.  
  91. freq=linspace(f0-fc,f0+fc,10001);
  92. port_calc=calcPort(port,Sim_Path,freq,'RefImpedance',50);
  93.  
  94. s11=port_calc{1}.uf.ref./port_calc{1}.uf.inc;
  95. s21=port_calc{2}.uf.ref./port_calc{1}.uf.inc;
  96. s12=port_calc{1}.uf.ref./port_calc{2}.uf.inc;
  97. s22=port_calc{2}.uf.ref./port_calc{2}.uf.inc;
  98.  
  99. ##[fid,msg]=fopen('out.s2p','w');
  100. ##wdata=[freq' real(s11)' imag(s11)' real(s21)' imag(s21)' real(s12)' imag(s12)' real(s22)' imag(s22)'];
  101. ##dlmwrite(fid,wdata,'delimiter','\t');
  102. ##fclose(fid);
  103.  
  104. endfunction

nfc_trx_analysis.m
  1. close all
  2. clear all
  3. clc
  4.  
  5. distance=[96 20 16 12 8 4 2 1];
  6.  
  7. for i=1:length(distance)
  8.   [freq_1 s11_1 s21_1 s12_1 s22_1]=nfc_trx_analysis_main(0,36,32,4,26,12,4,distance(i),1);
  9.   [freq_2 s11_2 s21_2 s12_2 s22_2]=nfc_trx_analysis_main(1,36,32,4,26,12,4,distance(i),1);
  10.  
  11.   fn=['out_' num2str(distance(i)) '.s2p']
  12.   #fn=['out_m45_' num2str(distance(i)) '.s2p']
  13.   [fid,msg]=fopen(fn,'w');
  14.   fprintf(fid,"# Hz S RI R 50\n");
  15.   wdata=[freq_1' real(s11_1)' imag(s11_1)' real(s21_1)' imag(s21_1)' real(s12_2)' imag(s12_2)' real(s22_2)' imag(s22_2)'];
  16.   dlmwrite(fid,wdata,'delimiter','\t');
  17.   fclose(fid);
  18. endfor

コードの構成がまぁまぁいびつなんだけど、しゃーない。OpenEMSで集中定数も含めてシミュレーションできるので、そのようなコードにしているけど、今回は無効になるような値を入れて無効化している。OpenEMSで2portのシミュレーションをするには、port1をexcite=trueにしたケースとport2をexcite=trueにしたケースの2回のシミュレーションが必要になる。んで、コイル間距離を8パタンやってしまおうという上記のコードを我が家のCorei5 Ubuntu24.02マシンで実行すると6.5時間くらいで完了。まぁ思ったより短い。てか短くするためにいろいろごまかしている。メタルはperfect electric conductor(PEC)で、厚さを0としている。誘電体は入れてないので真空にメタルのコイルが浮かんでるモデル。
で、結果を見てみる。
コイル間距離96mm(まぁ十分遠いのでほぼ結合していないと見ていい状態)のSパラファイルをQucs-sで見てみて、で、同調も取ってみる。
ってなるので、TXコイルは1548nH、RXコイルは686nHってことになる。あってるかどうか、、、例のサイトによると、TXコイルは1530nH、RXコイルは670nHってなるっぽい。まぁまぁあってるすごい。ちなみに、同調とってインダクタンスを推定するのはシミュレーションデータだからこそたぶんやっても大丈夫なやつで、実測データには寄生容量が含まれているので、このやり方ではなく、過去、日曜技術者がやったように低周波で \[ L=\frac{X}{\omega} \] をインダクタンス値とすべきである(はず)。

で、実際のものはどんな感じなのかって思う、、、完全に同じものは無理だけど、こんな感じで大体で作っちゃう。
んで、
TXアンテナは
ってことで1540nH。
RXアンテナは
ってことで690nH。
すごいね、よくあってる。むしろ工作が正確なのにびっくり。
では、2つのコイルを結合させてみる。とりあえず、20mm。
こんなかんじ。スペーサーボルトを駆使してちょうどコイル間が20mmになっている。
これを我が家のNanoVNA H4で測定するにはport1とport2を入れ替えて2回測定する必要がある。、、、やっぱまともなネットアナが欲しいなって思っちゃうけど、がまんも大事。
んで、見てみるわけだが、まずはOpenEMSのシミュレーション結果がこちら(注:Qucs-s回路図内の同調Cは無効化している)。
んで、実測結果はこちら。
いや、これがまた、よくあってるんよ。OpenEMSがまぁまぁ使えるってのと日曜技術者の器用さがまだまだ損なわれていないことに超びっくり。

んで、実測データに対してQucs-sで狙いの共振周波数を13.56MHzとして同調をとってみる。
送信側は76pF、受信側は191pFでそれぞれ直列共振周波数が13.56MHz、並列共振周波数が13.56MHzあたりになる。
で、この状態でACシミュレーションをしてみる。
これって入力している電圧源が1Vrmsに対して、出力が7.5Vrmsでるってことだと思うんだけど、ほんとかな?と思いながら今回をクローズする。これの実測はやらんといかんのはわかるけど、やる気が出るかどうか、、、

いやーさぼりにさぼった、、、このブログを購読している人はレアだろうから、まいっか。