VSim螺旋线行波管模拟案例讲义

1、 专业技术资料VSim螺旋线行波管模拟案例讲义第一部分 概述模拟的对象为一个典型的螺旋线行波管。基本结构是在金属管内装有螺旋线和介质夹持杆，然后在两端引入输入和输出端口。当电子束进入螺旋线波导时，就产生微波放大。实际模拟时我们用很多圈的螺旋线，为了看清几何，我们用一组示意图来描述：第二部分 几何生成从上面的图形可以看到，整个几何模型包括下面的主要内容：金属外壳：在外壳上安装了耦合端口，这里使用的耦合端口是一个简单的同轴线针头-壁面耦合，在模拟的时候直接从同轴线送入TEM波。然后是作为波导的螺旋线，同轴线针头直接耦合到螺旋线上。接下来，为了支持螺旋线，需要三根介质夹持杆：剖面图显示出夹持结构：实

2、际计算中，我们需要使用一个圈数很多的螺旋线（N=100 turns以上）。此外，在实践中还需要加入衰减器。我们首先把匝数调整到100，然后可以看到衰减器的结构（去掉外壳以便看清）：中央绿色的是衰减器。这里需要考虑一下VSim的逻辑，金属件都必须被实现为几何结构(gridBoundary)，而夹持杆（介电结构）和衰减器（电导耗散）既可以实现为几何结构，也可以直接用函数引入，它们并不是必须被写成gridBoundary的。确定一个结构是否要被实现为gridBoundary，要看它的具体应用。例如如果设置夹持杆和衰减器都要吸收电子，那么就应该将它们实现为gridBoundary，而如果我们暂时不考虑

3、夹持杆和衰减器的吸收，那么就可以到处理电磁场的时候再写这个部分。除了上述的部分之外，VSim逻辑建议把所有的粒子发射源(emitter)都设置成为gridBoundary。所以，在几何建模部分，需要建立金属部件和发射体的gridBoundary。特别是，如果整个系统中有若干组金属部件，但其中只有一个用作发射体，那么这个必须额外被单独几何建模一次，吸收体同理。gridBoundary可以用常规的对象来做，也可以使用geometry宏来做。下面是用geometry宏实现的金属部件几何，其中复杂的主要是螺旋线结构。基本思路就是从一个实体开始，不断贴上或者挖去部分，最后形成金属部件。$ DM_FRAC

4、 = 0.25$ import geometryfillGeoExpression(tube, geoCylinderXP(x-XBGN_GUN_WALL,y,z,REND_TUBE,LENGTH_TUBE_SHELL) )$ if PERIODIC_GEOMETRY=0fillGeoExpression(inputCoaxShell, geoCylinderXP(y-YBGN_COAX,z,x-XBGN_HELIX,RSHELL_COAX,LEGNTH_COAX) )fillGeoExpression(outputCoaxShell, geoCylinderXP(y-YBGN_COAX,z,

5、x-XEND_HELIX,RSHELL_COAX,LEGNTH_COAX) )$ endifvoidGeoExpression(tubeInterior, geoCylinderXP(x-XEND_GUN_WALL,y,z,RBGN_TUBE,LENGTH_TUBE_INTERIOR) )$ if PERIODIC_GEOMETRY=0voidGeoExpression(inputCoaxInterior, geoCylinderXP(y-YBGN_COAX,z,x-XBGN_HELIX,ROUTER_COAX,LEGNTH_COAX) )voidGeoExpression(outputCoa

6、xInterior, geoCylinderXP(y-YBGN_COAX,z,x-XEND_HELIX,ROUTER_COAX,LEGNTH_COAX) )$ endiffillGeoExpression(helix, geoCylinderXP(x-XBGN_HELIX_WIRE,yHelix(y,z),zHelix(x,y,z),RWIRE_HELIX,LENGTH_HELIX_WIRE) )$ if PERIODIC_GEOMETRY=0voidGeoExpression(trimHelixInput, myHalfCylinderXPZ(x-XBGN_HELIX_WIRE,y,-z,R

7、BGN_TUBE,LENGTH_TRIM_HELIX_INPUT) )voidGeoExpression(trimHelixOutput, myHalfCylinderXPZ(x-XTRIM_HELIX_OUTPUT,y,z,RBGN_TUBE,LENGTH_TRIM_HELIX_OUTPUT) )fillGeoExpression(inputCoaxInner, geoCylinderXP(y-YBGN_COAX,z,x-XBGN_HELIX,RINNER_COAX,LEGNTH_COAX) )fillGeoExpression(outputCoaxInner, geoCylinderXP(

8、y-YBGN_COAX,z,x-XEND_HELIX,RINNER_COAX,LEGNTH_COAX) )$ endifsaveGeoToGridBoundary(metalParts,DM_FRAC)其中有一部分几何用$ if PERIODIC_GEOMETRY=0来设置为条件形式。这些部分主要是两端的切断以及馈电的同轴线接口。涉及到的几个宏的具体描述参考VSim In Depth。emitter几何可以类似地来做，应该指出的是emitter一般来说不必和金属部件有什么关系，但是实践中emitter往往是金属部件的一部分。在这种情况下，emitter必须被独立地实现为一个gridBounda

9、ry：resetGeoToVoid()fillGeoExpression(emitterDisk, geoCylinderXP(x-XBGN_GUN_WALL,y,z,RBGN_TUBE,THICKNESS_ENDCAPS) )saveGeoToGridBoundary(emitter,DM_FRAC)第三部分：场模型和频率分析基本框架产生金属部件模型后，就可以启用电磁场模型了。电磁场模型看上去很复杂，实际上基本都是从一个模板产生的，最简单的模板是以前用的dMYee宏。对于微波器件，一般我们可以从VSim6的例子中选择一个类似的模型，然后逐步修改成我们需要的模型。这里我们使用的是一个比较通用的

10、模板，它基本来自于dMYee宏。# 场定义.#推进公式# - - - - - - - - - - - - - - - - - - - - - - - - -# update operations# 实际的推进步骤通常几乎不会有人真的从头写出整个电磁模块，而都是在一个类型比较相似的模块上进行修改。当然，我们可以使用以前遗留的一个宏deymittra.mac,此外现在的TWT输入文件也是一个不错的模板。这个模板是利用MultiField写成的。为了在需要的时候修改这个模板，我们需要解释其基本内容。一般来说，MultiField框架要包含三个主要的部分：场定义，更新器和算法包装。由于MultiFie

11、ld是一个具有相当自由的场计算框架，原则上我们应该首先列出场方程。实践中因为微波器件的场模型是非常类似的，一般只要简单修改就可以，无需如此复杂。这里为了解释清楚，我们列出基本概念。微波器件计算的电磁场，一般包含静态外场（主要是外磁场）和微波时变电磁场。时变电磁场的方程是和在计算电磁学中的方法一样,电场强度（和电流强度）被定义在网格的棱边上，而磁场强度B被定义在面心，如同下面这个图：和标准计算电磁学稍有区别的是，由于存在粒子电流J，为了计算的精度，上述方程是被拆成三步来进行FDTD推进。按照VSim的逻辑，我们引入一个参数toDtFrac，它的取值从0.0到1.0。我们假设每个推进步推进中，迭代

12、器的绝对时间从推进到，每个电磁场推进步的时间步长是。实际中MultiField总是要由一组“更新步”组成，而每个更新步代表从“当前时间”推进到某个值。现在我们可以列出PIC算法里面最标准的电磁场推进步，推进开始的时刻是，所有场的“内部时间”都是：第一步：进行磁场B的推进。设置toDtFac参数为0.5，也就是推进目标时间然后按照法拉第电磁感应定律将磁场推进为其中，其中是磁场在这个更新步之前的内部时间。这里显然。于是，于是磁场被推进了半个时间步，推进结束后，磁场内部时间被推进到。第二步：进行电场E的推进，这是利用广义的安培定律：由于电场只有一步推进，因此这里的toDtFrac=1.0，所以，并且

13、在执行这个推进后，电场的内部时间被推进到。第三步：再次推进磁场，仍然使用法拉第定律，并且设置toDtFrac=1.0，于是并且。执行这个推进后，磁场内部时间被推进到。现在我们可以分析行波管的MultiField模型了。在开始的部分，必须定义这个MultiField框架所蕴含的场：numComponents = 3offset = edgeoverlap = 1 2numComponents = 3offset = faceoverlap = 1 2kind = depFieldnumComponents = 3offset = edgeoverlap = 1 2$ if INCLUDE_PAR

14、TICLES=1numComponents = 3offset = noneoverlap = 1 2dumpPeriod = 0numComponents = 3offset = noneoverlap = 1 2dumpPeriod = 0$ endif几个场的语法是一样的，都是要定义场的分量，位置(offset参数用来约定场量定义是在顶点，面心还是棱边上），存盘周期以及鬼网格结构。在这里首先定义了电场，磁场和电流场。接下来两个节点电磁场nodalE和nodalB是用来将计算出的电磁场传递给带电粒子计算。所以这个部分是使用了条件编译。接下来是更新算法，我们看到一系列的FieldUpdate

15、r和FieldMultiUpdater：gridBoundary = metalPartsinteriorness = deymittrakind = yeeFaradayUpdatercomponents = 0 1 2expandToTopInComponentDir = 1lowerBounds = 0 0 0upperBounds = NX NY NZreadFields = edgeEwriteFields = faceBgridBoundary = metalPartsinteriorness = deymittrakind = deyMittraUpdaterlowerBound

16、s = 0 0 0upperBounds = NX NY NZreadFields = edgeEwriteFields = faceBgridBoundary = metalPartsinteriorness = deymittrakind = yeeAmpereUpdatercomponents = 0 1 2contractFromBottomInNonComponentDir = 1lowerBounds = 0 0 0upperBounds = NX NY NZreadFieldCompShifts = 0 0readFields = faceB edgeJwriteFields =

17、 edgeE$ if INCLUDE_PARTICLES=1kind = edgeToNodeVeclowerBounds = 0 0 0upperBounds = $NX+1$ $NY+1$ $NZ+1$readFields = edgeEwriteFields = nodalEkind = faceToNodeVeclowerBounds = 0 0 0upperBounds = $NX+1$ $NY+1$ $NZ+1$readFields = faceBwriteFields = nodalB$endif这里面有五个Updater。一般来说前三个都是标准的：yeeFaradayUpdat

18、er和deyMittraUpdater用来推进电磁场的法拉第感应定律：由于导体表面处没有切向电场，需要对这个方程做特殊考虑，所以法拉第定律包含了两个updater，yeeFaradayUpdater用于计算不和导体接触的网格，而deyMittraUpdater用于计算导体边界处的网格，两个updater一般都是在一起应用的。yeeAmpereUpdatery用来推进安培定律：上面两种updater足以完成FDTD推进步骤，为了把电磁场输出给粒子，需要把edgeE和faceB的值传送给nodalE和nodalB，这是后面两个updater完成的工作。一般情况下我们基本不需要修改这些updater

19、，只要补充我们需要的updater并且封装就可以了。最后是包装成完整的推进步：toDtFrac = 0.5messageFields = faceBupdaters = yeeFaraday deyMittraFaradaytoDtFrac = 1.0updaters = yeeAmperemessageFields = edgeEtoDtFrac = 1.0messageFields = faceBupdaters = yeeFaraday deyMittraFaraday$ if INCLUDE_PARTICLES=1toDtFrac = 1.0updaters = nodalEupdat

20、emessageFields = nodalEtoDtFrac = 1.0updaters = nodalBupdate messageFields = nodalB$ endif按照我们说的那样，首先是finishFarady（半步法拉第定律），然后是Ampere(一步安培定律），再startFarady（另外半步法拉第定律）。实际上完成这三步之后，一个FDTD推进步就完成，剩下的两个更新步用于将E和B输出到nodalE和nodalB。加入介质和损耗上面的计算只包含了导体部分，而实际的行波管还包含夹持杆和衰减器。从模拟的角度来说夹持杆是一块电介质，而衰减器是一块部分导体。对于电介质和衰减介质

21、的建模，VSim有两种方案，一种是使用现成的两个求解器（更新器），另一种是干脆自己写一个更新器。两者的区别是前者可以处理复杂外形（包括CAD导入），但是只能是均匀的介电函数或者电导函数；后者只能处理用函数描述的结构，但可以处理介电/电导函数的演化。这里因为外形很简单，我们使用第二种方法。首先看一下处理介电和电导问题的方法，在考虑介电和电导的情况下，安培定律是使用半点差分，这个方程可以写成化简后就是或者写成其中第二个方程完全等价于原来的安培更新步，所以我们只需要在原来的安培更新步之前之后各插入一个转换步，就可以实现带电导和介电的器件的求解。为了实现这两个步骤，首先要写出介电和电导的表达式，这两个

22、表达式是： 1.0 + (EPS_ROD-1.0)*geoCylinderXP(x-XBGN_HELIX_WIRE,y-YROD1,z-ZROD1,RADIUS_ROD,LENGTH_HELIX_WIRE) + (EPS_ROD-1.0)*geoCylinderXP(x-XBGN_HELIX_WIRE,y-YROD2,z-ZROD2,RADIUS_ROD,LENGTH_HELIX_WIRE) + (EPS_ROD-1.0)*geoCylinderXP(x-XBGN_HELIX_WIRE,y-YROD3,z-ZROD3,RADIUS_ROD,LENGTH_HELIX_WIRE)# definit

23、ions for attenuator$ VPHASEX_HELIX = LIGHTSPEED*SINPHI$ XBGNTURN_ATTENUATOR = XBGN_HELIX + BGNTURN_ATTENUATOR*PITCH_HELIX$ LENGTH_ATTENUATOR = (ENDTURN_ATTENUATOR - BGNTURN_ATTENUATOR)*PITCH_HELIX$ SIGMA = (VPHASEX_HELIX*EPS0/LENGTH_ATTENUATOR)*ln(10.0*(DB_ATTENUATOR/10.0)$ ATTENUATE_FACTOR = 0.5*DT

24、*SIGMA/EPS0 ATTENUATE_FACTOR*geoPipeXP(x-XBGNTURN_ATTENUATOR,y,z,RMID_HELIX,REND_TUBE,LENGTH_ATTENUATOR)然后需要在MultiField框架内加入updater和updaterStep: # = epsRel - 0.5*sigma*dt/eps0numComponents = 3overlap = 1 2offset = edgekind = variablelowerBounds = 0 0 0upperBounds = NX NY NZcomponents = 0 1 2kind = e

25、xpressionexpression = epsRelTerm(x,y,z) - attenuateTerm(x,y,z)kind = expressionexpression = epsRelTerm(x,y,z) - attenuateTerm(x,y,z)kind = expressionexpression = epsRelTerm(x,y,z) - attenuateTerm(x,y,z)# = epsRel + 0.5*sigma*dt/eps0numComponents = 3overlap = 1 2offset = edgekind = variablelowerBound

26、s = 0 0 0upperBounds = NX NY NZcomponents = 0 1 2kind = expressionexpression = epsRelTerm(x,y,z) + attenuateTerm(x,y,z)kind = expressionexpression = epsRelTerm(x,y,z) + attenuateTerm(x,y,z)kind = expressionexpression = epsRelTerm(x,y,z) + attenuateTerm(x,y,z)kind = fieldBinOpUpdaterlowerBounds = 0 0

27、 0upperBounds = NX NY NZbinOp = multiplyaCoeff = 0.0bCoeff = 0.0readFields = edgeE epsRelMinus writeFields = edgeEkind = fieldBinOpUpdaterlowerBounds = 0 0 0upperBounds = NX NY NZbinOp = divideaCoeff = 0.0bCoeff = 0.0readFields = edgeE epsRelPlus writeFields = edgeE上面是电磁场的基本部分。按照我们知道的，行波管除了等离子体自生电磁场

28、之外，还有一个导向磁场，为此我们还要再加入一个更新器：kind = STFuncUpdateroperation = addlowerBounds = 0 0 0upperBounds = $NX+1$ $NY+1$ $NZ+1$writeFields = nodalBcomponent = 0kind = expressionexpression = BSTATIC然后我们需要把这些updater安装在正确的位置，首先我们装入EtoD和DtoE：toDtFrac = 1.0updaters = EtoD yeeAmpere DtoEmessageFields = edgeE然后装入静磁场：t

29、oDtFrac = 1.0updaters = nodalBupdate addInBstaticmessageFields = nodalB最后还有一个问题，即行波管具有输入和输出端口。在VSim6中，最简单的实现输入和输出端口的方法是使用Port和incident宏。这两个宏必须使用在ampereUpdater附近：applyPortAtBoundary_beforeAmpere(1,1,LIGHTSPEED,DY,DT,edgeE,NX,NY,NZ).applyPortAtBoundary_afterAmpere(1,1,LIGHTSPEED,DY,DT,edgeE,NX,NY,NZ)a

30、ddIncidentWaveAtBoundaryPort(1,1,LIGHTSPEED,DY,DT,edgeE,NX,NY,NZ,0,0,inputWaveX(x,y,z,t)addIncidentWaveAtBoundaryPort(1,1,LIGHTSPEED,DY,DT,edgeE,NX,NY,NZ,2,0,inputWaveZ(x,y,z,t)此外，这两个宏都包含在port文件里面，所以要在文件的开头加入$import port原理上，这样行波管的cold部分就设计完成了，接下来，我们可以加入一些诊断来检查电磁波会不会在行波管内正常传播。例如，加入出射信号诊断：# input sign

31、als$ IXVINPUTBGN = int(XBGN_HELIX-RSHELL_COAX-XBGN)/DX)$ IXVINPUTMID = int(XBGN_HELIX -XBGN)/DX)$ IXVINPUTEND = int(XBGN_HELIX+RSHELL_COAX-XBGN)/DX)kind = pseudoPotentialfield = yeeEM.edgeEreferencePoint = IXVINPUTMID $NY-1$ $int(NZ/2)$measurePoint = $IXVINPUTEND+1$ $NY-1$ $int(NZ/2)$kind = fieldPoy

32、nlowerBounds = IXVINPUTBGN $NY-1$ 0upperBounds = $IXVINPUTEND+1$ $NY-1$ NZfields =yeeEM.edgeE, yeeEM.faceB# output signals$ IXVOUTPUTBGN = int(XEND_HELIX-RSHELL_COAX-XBGN)/DX)$ IXVOUTPUTMID = int(XEND_HELIX -XBGN)/DX)$ IXVOUTPUTEND = int(XEND_HELIX+RSHELL_COAX-XBGN)/DX)$ IYVINPUT = NY-1kind = pseudo

33、Potentialfield = yeeEM.edgeEreferencePoint = IXVOUTPUTMID $NY-1$ $int(NZ/2)$measurePoint = $IXVOUTPUTEND+1$ $NY-1$ $int(NZ/2)$kind = fieldPoynlowerBounds = IXVOUTPUTBGN $NY-1$ 0upperBounds = $IXVOUTPUTEND+1$ $NY-1$ NZfields = yeeEM.edgeE, yeeEM.faceB诊断的内容包括输入输出功率/电压。其中电压需要特殊注意。严格的说电磁波器件并没有电压，实践中VSim

34、用“伪势”量度两点之间的压差。在现在的情况下，输入输出电压就是在输入和输出端口同轴电缆芯部和外皮上各任选一点然后测量的伪势。输入电压输出电压第四部分：粒子加入和放大性能物理建模的最后一步是加入粒子。由于没有背景等离子体，所以只要加入粒子发射和吸收就可以了。引入粒子的基本模式就是先写一个标准的Species： $ NPARTICLES_PER_CELL_PER_TIME_STEP = 0.1 $ VEMIT = sqrt(2.0*BEAM_VOLTAGE/511.0e3)*LIGHTSPEED $ BEAM_AREA = PI*BEAM_RADIUS*2 $ JDENS = BEAM_CURRE

35、NT/BEAM_AREA $ NSTEPS_ACROSS_CELL = DX/(VEMIT*DT) $ NOMINAL_DENSITY = JDENS/(VEMIT*1.602176487e-19) $ NOMINAL_PARTICLES_PER_CELL = NPARTICLES_PER_CELL_PER_TIME_STEP*NSTEPS_ACROSS_CELL kind = relBorisVW charge = -1.602176487e-19 mass = 9.10938215e-31 nominalDensity = NOMINAL_DENSITY nomPtclsPerCell =

36、 NOMINAL_PARTICLES_PER_CELL fields = yeeEM.nodalE yeeEM.nodalB currDeps = Jdepostion.这里有一个需要特殊注意的地方，是我们必须将粒子和电磁场耦合起来。耦合包含两部分，一个是让粒子接受电磁场力的作用，这是通过fields = yeeEM.nodalE yeeEM.nodalB 这样语法实现的；另一点则是将粒子的电流密度传递给电磁场，这是通过currDeps语句。为了让currDeps生效，我们还要对其中的Jdeposition进行定义： kind = areaWeighting depField = yeeEM.

37、edgeJ 其中areaWeighting 表示使用面积权重计算，具体到这个语法的详细内容可以研究手册，我们这里只指出一点：如果使用的是MultiField类型的场，那么我们的电流是自己定义的，因此必须使用一个独立的VectorDepositior；如果按照以前的方法使用emMultiField，那么电流场自动定义成SumRhoJ，也就无需自己再定义VectorDepositior，当然对应地也不需要定义edgeJ。定义了粒子类型后，需要引入一个粒子源。在现在的情况下，发射体已经定义，所以直接使用一个cutCellPosGen的xvLoaderEmitter就好。一般我们总是xvLoaderE

38、mitter: kind = xvLoaderEmitter positionFunction = bitReversed useCornerMove = true emitBasedOnLocalForce = false load = false applyTimes = 0.0 1.0 kind = cutCellPosGen numMacroPtclsPerStep = $int(NPARTICLES_PER_CELL_PER_TIME_STEP*BEAM_AREA/(DY*DZ)$ emitterBoundary = emitter emissionOffset = 0.01 kind=expression expression = H(x-(XEND_GUN_WALL-DX)*H(BEAM_RADIUS*2-y*2-z*2)-0.5 kind = funcVel