ON THE ROBUSTNESS OF A MULTIGRID METHOD FOR ANISOTROPIC REACTION-DIFFUSION PROBLEMS
Robust method to retrieve the constitutive effective parameters of metamaterials
Robust method to retrieve the constitutive effective parameters of metamaterials Xudong Chen,Tomasz M.Grzegorczyk,Bae-Ian Wu,Joe Pacheco,Jr.,and Jin Au Kong Research Laboratory of Electronics,Massachusetts Institute of Technology,Cambridge,Massachusetts02139,USA(Received25February2004;published26July2004)We propose an improved method to retrieve the effective constitutive parameters(permittivity and perme-ability)of a slab of metamaterial from the measurement of S parameters.Improvements over existing methods include the determination of thefirst boundary and the thickness of the effective slab,the selection of the correct sign of effective impedance,and a mathematical method to choose the correct branch of the real part of the refractive index.The sensitivity of the effective constitutive parameters to the accuracy of the S parameters is also discussed.The method has been applied to various metamaterials and the successful retrieval results prove its effectiveness and robustness.DOI:10.1103/PhysRevE.70.016608PACS number(s):41.20.Jb,42.25.Bs,78.20.CiI.INTRODUCTIONLeft-handed(LH)structures have been realized so far asmetamaterials[1–3]and very quickly,researchers have beenworking on retrieving their effective permittivity and perme-ability to better characterize them[4–6].Known methods todate[7,8]use S parameters calculated from a wave incidentnormally on a slab of metamaterial,from which the effectiverefractive index n and impedance z arefirst obtained.Thepermittivity⑀and permeabilityare then directly calculated from=nz and⑀=n/z.It is also known that this retrieval process may fail insome instances,such as when the thickness of the effectiveslab(exhibiting bulk properties)is not accurately estimated [4]or when reflection͑S11͒and transmission͑S21͒data are very small in magnitude[6].Although these issues have beenaddressed to some extent in previous works,we have foundthat the retrieved results in some cases are still unsatisfac-tory.In particular,the determination of thefirst boundary ofthe effective homogeneous slab,the selection of the sign ofthe roots when solving for the impedance z,the determina-tion of the branch of the real part of refractive index n,andthe origin of the spikes appearing in the retrieved impedancez,are many problems that deserve further investigation.Theaforementioned issues are addressed in the next sections ofthis paper and some typical retrieval results are presented toshow the robustness and effectiveness of the method.Notethat the values of⑀,,and z are relative to those in free space,thus dimensionless.II.RETRIEV AL METHODA.Theoretical retrieval equationsIn order to retrieve the effective permittivity and perme-ability of a slab of metamaterial,we need to characterize it as an effective homogeneous slab.In this case,we can retrieve the permittivity and permeability from the reflection͑S11͒and transmission͑S21͒data.For a plane wave incident nor-mally on a homogeneous slab of thickness d with the origin coinciding with thefirst face of the slab,S11is equal to the reflection coefficient,and S21is related to the transmission coefficient T by S21=Te ik0d,where k0denotes the wave num-ber of the incident wave in free space.The S parameters are related to refractive index n and impedance z by[7,9,10]:S11=R01͑1−e i2nk0d͒1−R012e i2nk0d,͑1a͒S21=͑1−R012͒e ink0d1−R012e i2nk0d,͑1b͒where R01=z−1րz+1.As it has been pointed out in[4,5],the refractive index n and the impedance z are obtained by inverting Eqs.(1a)and (1b),yieldingz=±ͱ͑1+S11͒2−S212͑1−S11͒2−S212,͑2a͒e ink0d=X±iͱ1−X2,͑2b͒where X=1ր2S21͑1−S112+S212͒.Since the metamaterial un-der consideration is a passive medium,the signs in Eqs.(2a) and(2b)are determined by the requirementzЈജ0,͑3a͒nЉജ0,͑3b͒where͑·͒Јand͑·͒Љdenote the real part and imaginary part operators,respectively.The value of refractive index n can be determined from Eq.(2b)asn=1k0d͕͓͓ln͑e ink0d͔͒Љ+2m͔−i͓ln͑e ink0d͔͒Ј͖,͑4͒where m is an integer related to the branch index of nЈ.As mentioned above,the imaginary part of n is uniquely deter-mined,but the real part is complicated by the branches of the logarithm function.Equations(2a)and(2b)can be directly applied in the case of a homogeneous slab for which the boundaries of the slab are well defined and the S parameters are accurately known. However,since a metamaterial itself is not homogeneous,the two apparently straightforward issues mentioned above need to be carefully addressed.First,the location of the twoPHYSICAL REVIEW E70,016608(2004)boundaries of the effective slab need to be determined,whichwe do here by ensuring a constant impedance for various slab thicknesses.Second,the S parameters obtained from numerical computation or measurements are noisy which can cause the retrieval method to fail,especially at those frequen-cies where z and n are sensitive to small variations of S 11and S 21.These two problems are examined in detail in the fol-lowing sections.B.Determination of the first boundary and the thicknessof the effective slabA homogeneous slab of material can be characterized bythe fact that its impedance does not depend on its thickness.Our understanding of the physical meaning of the first effec-tive boundary is a plane beyond which the reflected wave behaves like a plane wave for a plane wave incidence.When a plane wave is incident on a metamaterial,currents will be induced on the metals creating a scattered field.The field produced by the induced currents is not uniform:It is stron-gest around the metal and decay at a certain distance.The first effective boundary is located where the reflected wave behaves like a plane wave,and it has to be determined.We use z 1and z 2to represent the impedances of two slabs ofmetamaterial of different thicknesses.The reflection S 11de-pends on the position of the first boundary and the transmis-sion S 21depends on the thickness of the slab.In addition,since the impedance z is also a function of S 11and S 21,z depends on the first boundary and the thickness of the slab as well.Taking into account the above-mentioned properties,we propose a method whereby the first boundary and the thickness of the sample are determined by optimization so that z 1matches z 2at all frequencies.Figure 1illustrates the configuration of the problem for a metamaterial made of two cells in the propagation direction (x direction ).The geometry of the metamaterial has been taken from [11,12],in which the dimensions have been slightly modified for ease of dis-cretization in finite-difference time-domain (FDTD )simula-tions.With the split-ring resonator (SRR )and rod in the cen-ter of the unit cell,the periodicity along the propagation direction is d 0.The first boundary of the effective homoge-neous medium is located at x 1below ͑x 1ജ0͒or above ͑x 1Ͻ0͒the first unit-cell boundary,and the thickness of the effective medium is 2d 0+x 2−x 1.The optimization model is set up to minimize the mismatch of impedances of different numbers of cells of metamaterial:min f ͑x ¯͒=1N f ͚i =1Nf͉z 1͑f i ,x ¯͒−z 2͑f i ,x ¯͉͒max ͕͉z 1͑f i ,x ¯͉͒,͉z 2͑f i ,x ¯͉͖͒,s.t.:−0.5d 0ഛx 1,x 2ഛ0.5d 0,x ¯=͑x 1,x 2͒,͑5͒where N f is the total number of sample frequencies and z j ͑f i ͒represents the impedance of slab j ͑j =1,2͒at frequency f i .In the ideal case,z 1matches z 2for all frequencies with the objective function value equal to zero.We use the differential evolution algorithm [13]to optimize the objective function,and the optimized solution is x ¯opt =͑3.8565ϫ10−4d 0,1.0479ϫ10−4d 0͒.The corresponding values of impedance are shown in Fig.2.It can be seen that the retrieved impedances for one,two,and three cells of SRR-rod metamaterial match well for most frequencies,while matching was not as satis-factory when the method in Ref.[4]was used (which corre-sponds to x 1=0.5d 0in our formulation ).We also calculated the impedance z for the case of x ¯=͑0,0͒and found that the results are almost the same as the optimized ones.We have corroborated this result with many other metamaterial thick-nesses and geometries to eventually conclude empirically that the first effective boundary of a symmetric one-dimensional (1D )metamaterial [1,4,14]coincides with the first unit-cell boundary and the second effective boundary coincides with the last unit-cell boundary.Fortwo-FIG.1.Illustration of the effective boundaries of a two-cellmetamaterial.The SRRs and rods are periodic along yˆand z ˆdirec-tions with periodicity a y =4mm,a z =3mm.The unit-cell thickness ͑d 0͒in the direction of wave incidence is 4mm.We choose the first and the last unit-cell boundary as the reference plane for the param-eters x 1and x 2,respectively.The value of x 1and x 2are positive (negative )if the dashed lines are below (above )their reference planes.The thickness of the effective homogeneous medium is 2d 0+x 2−x 1͑mm ͒.FIG.2.Optimized impedance z for one,two,and three cells of SRR-rod metamaterial of Fig.1in the direction of propagation.CHEN et al.PHYSICAL REVIEW E 70,016608(2004)dimensional (2D )[11,14]and asymmetric 1D metamaterials,no such rule could be found and the effective boundaries of the slab need to be determined from optimization.C.Determination of n and z from S 11and S 21It is a common method to determine z and n from Eqs.(2a )and (2b )with the requirement of Eqs.(3a )and (3b ),where z and n are determined independently.However this method may fail in practice when z Јand n Љare close to zero:A little perturbation of S 11and S 21,easily produced in experi-mental measurements or numerical simulations,may change the sign of z Јand n Љ,making it unreliable to apply the re-quirement of Eqs.(3a )and (3b ),as discussed in Ref.[6].In fact,z and n are related and we should use their relationship to determine the signs in Eqs.(2a )and (2b ).In order to determine the correct sign of z ,we distinguish two cases.The first is for ͉z Ј͉ജ␦z ,where ␦z is a positive number,for which we apply Eq.(3a ).In the second case,for ͉z Ј͉Ͻ␦z ,the sign of z is determined so that the corresponding refractive index n has a non-negative imaginary part,or equivalently ͉e ink 0d ͉ഛ1,where n is derived from Eqs.(1a )and (1b ):e ink 0d =S 211−S 11z −1z +1.͑6͒Note that once we obtain the value of z ,the value of e ink 0d is obtained from Eq.(6),so that we avoid the sign ambiguity in Eq.(2b )[it can be proven that only one sign of the imagi-nary part in Eq.(2b )makes it equivalent to Eq.(6)].Figure 3shows the retrieved impedance using the method presented in this paper and using only the condition of Eq.(3a ).It is noted that the discontinuities obtained when only applying the criterion z Јജ0are removed.D.Determination of the branch of n ЈWe have presented in the previous sections a method of solving for z and n Љ,but n Јremains ambiguous because ofthe branches of logarithm function as seen in Eq.(4).In order to address this problem,it has been suggested to use a slab of small thickness and applying the requirement that ⑀͑f ͒and ͑f ͒are continuous functions of frequency [4,5].However,no further details on the continuity argument were provided.In our method,we determine the proper branch by using the mathematical continuity of the parameters,with special attention to possible discontinuities due to reso-nances.The method is an iterative one:Assuming we have obtained the value of the refractive index n ͑f 0͒at frequency f 0,we obtain n ͑f 1͒at the next frequency sample f 1by ex-panding the function e in ͑f 1͒k 0͑f 1͒d in a Taylor series:e i n ͑f 1͒k 0͑f 1͒d Ϸe i n ͑f 0͒k 0͑f 0͒d ͩ1+⌬+12⌬2ͪ,͑7͒where ⌬=i n ͑f 1͒k 0͑f 1͒d −i n ͑f 0͒k 0͑f 0͒d and k 0͑f 0͒denotes the wave number in free space at frequency f 0.In Eq.(7),the branch index [m in Eq.(4)]of the real part of n ͑f 1͒is the only unknown.Since the left-hand side of Eq.(7)is obtained from Eq.(6),Eq.(7)is a binomial function of the unknown n ͑f 1͒.Out of the two roots,one of them is an approximation of the true solution.Since we have obtained n Љ͑f 1͒,we choose the correct root among the two by com-paring their imaginary parts with n Љ͑f 1͒.The root whose imaginary part is closest to n Љ͑f 1͒is the correct one,and we denote it as n 0.Since n 0is a good approximation of n ͑f 1͒,we choose the integer m in Eq.(4)so that n Ј͑f 1͒is as close to n 0Јas possible.The branch of n Јat the initial frequency is determined as follows:From =nz and ⑀=n /z ,we haveЉ=n Јz Љ+n Љz Ј,͑8a ͒⑀Љ=1͉z ͉2͑−n Јz Љ+n Љz Ј͒.͑8b ͒The requirements Љജ0and ⑀Љജ0leadtoparison of the retrieved imped-ance z (real and imaginary parts )for one cell of metamaterial shown in Fig.1by the method pre-sented in this paper and a traditional method us-ing only the requirement z Јജ0.ROBUST METHOD TO RETRIEVE THE CONSTITUTIVE …PHYSICAL REVIEW E 70,016608(2004)͉n Јz Љ͉ഛn Љz Ј.͑9͒In particular,when n Љz Јis close to zero but z Љis not,n Јshould be close to zero.At the initial frequency,we solve for the branch integer m satisfying Eq.(9).If there is only one solution,it is the correct branch.In case of multiple solu-tions,for each of the candidate branch index m ,we deter-mine the value of n Јat all subsequent frequencies using the above-mentioned iterative method.Because the requirement of Eq.(9)applies to n Јat all frequencies,we use it to check the validity of n Јat all frequencies produced by the candi-date initial branch.Note in the special case when n Љz Јis close to zero but z Љis not,the checking process can easily be carried out.Therefore,the initial branch that satisfies Eq.(9)at both the initial frequency and the subsequent frequencies is the correct one.For the SRR-rod structure,we found that there is a fre-quency region at which there is no branch index m satisfying Eq.(9).We call this frequency region the resonance band.The properties of the resonance band are still disputed by researchers.Some papers [15–17]mention the existence of multiple modes in this region since the real part of n is large,yielding a wavelength comparable to or smaller than the unit size of the metamaterial thereby rendering the retrieval of the effective parameters of the metamaterials impossible.Other papers [5,18]state that retrieval is possible and the retrieved permittivity ⑀has a negative imaginary part in the resonance band.In this paper,we do not address this issue and for this reason the retrieved results presented here are interrupted in frequency by the resonance region (see,for example,Fig.4).In this case,since n ͑f ͒is not continuous through all frequen-cies,we have to determine the initial branches for two fre-quency regions:Below and above the resonance band.Note that below the resonant band,the retrieved branch index is zero,which confirms the validity of the traditional method used for low-frequency retrieval.The retrieved refractive in-dexes n for one,two,and three cells in the propagation di-rection are shown in Fig.4,where the resonance band is seen to extend between 11GHz and 12GHz.We observe that thevalues of n for different cell numbers are identical above the resonant region.Below the resonant band,however,the re-trieved n for one and two cells match well,but the result for three cells differs significantly from the previous two.This discrepancy is due to the small magnitude of S 21in this fre-quency band,as we shall discuss in the next section.E.Sensitivity analysisA close examination of the retrieved z and n for one,two,and three cells of metamaterial presented so far shows that the three results do not always match.There are two cases of discrepancies.The first is that the retrieved refractive index n for three cells of metamaterial does not match the value for one and two cells at low frequencies (5GHz–11GHz in Fig.4).The second is that the impedance z appears to spike at some frequencies (around 12GHz,17GHz,and 19.5GHz in Fig.2).We shall show here that these discrepancies are due to the sensitivity of z and n to the accuracy of S 11and S 21.The first case appears when ͉S 21͉is close to zero.In the region below the resonance band,the transmission is usually small,especially for thicker metamaterials.From Eq.(2b ),we see that S 21appears in the denominator,so that the values of n are very sensitive to small perturbations of S 21.Yet,a small transmission has little influence on the retrieval of z ,which can be seen by computing:ץz 2ץS 21=8S 21S 11͓͑1−S 11͒2−S 212͔2,͑10͒from which it is clear that a small ͉S 21͉makes ץz 2րץS 21small (approximately zero ).In addition,we can see from Eq.(1b )that if n Љis not small,S 21will decrease exponentially with d ,i.e.,with an increasing number of cells in the propa-gation direction.Therefore,the smaller S 21,the higher the computation and measurement relative errors,which leads to less accurate retrieval results.The second case appears when S 212is close to unity while S 11is close to zero.Similar to the first case,the denominator in the expression of z [see Eq.(2a )]approaches zero,thus making it difficult to retrieve z .However,in this case,the value of n is stable.In this situation,instead of solving fornFIG.4.Retrieved refractive index n (real and imaginary parts )for one,two,and three cells of the metamaterial structure shown in Fig.1.FIG.5.Range of z (real and imaginary parts )for tolerance ␦r =0.015and ␦t =0.0in Eqs.(11a )and (11b ).The impedance is for a three-cell metamaterial shown in Fig.1.CHEN et al.PHYSICAL REVIEW E 70,016608(2004)and z which exactly satisfy Eqs.(1a )and (1b ),we solve for the following inequalities:ͯS 11−R 01͑1−e i2nk 0d ͒1−R 012e i2nk 0dͯഛ␦r ,͑11a ͒ͯS 21−͑1−R 012͒e ink 0d 1−R 012ei 2nk 0dͯഛ␦t ,͑11b ͒where ␦r and ␦t are small positive numbers.Figure 5shows the range of z satisfying Eqs.(11a )and (11b )for ␦r =0.015and ␦t =0.0.At each frequency,all z having a real and imagi-nary parts between the bounds shown in Fig.5satisfy Eqs.(11a )and (11b ).It can be seen that the magnitude of the spikes is within the tolerance error,which implies that they are due to the noise in the S 11and S 21data.Finally,note that although the retrieved n and z for mul-tiple cells may be different from that for one cell at some specific frequencies,the calculated S 11and S 21for multiple cells using the retrieved ⑀and from one cell data match well with the S 11and S 21data computed for multiple cells directly from numerical simulation,as is illustrated in Fig.6.F.ResultsThe retrieved permittivity ⑀and permeability of a one cell of SRR-rod structure of Fig.1are shown in Fig.7.Notethat although the results satisfy the condition ⑀Љജ0and Љജ0,the positive energy requirement ץ͑⑀͒/ץϾ0[19,20]is violated in the frequency band 12GHz–12.2GHz.As a re-sult,the resonance band is extended to 11GHz–12.2GHz,as shown by the vertical dashed lines in Fig.7(a ).The value of ⑀and are both negative in the frequency range 12.2GHz–12.8GHz,showing that in this band,the metamaterial exhibits a LH behavior.We also retrieved the effective parameters of four and five cells of metamaterial shown in Fig.1,and the retrieval results are close to those for one,two,and three cells.In addition,we also applied our method to retrieve the effective parameters of the structure taken from [14,21],as shown in the inset of Fig.8(a ).For a 1D structure,by match-ing the impedance z for one and two cells of the metamate-rial using the previously described method,we obtain x ¯opt=͑2.2053ϫ10−3d 0,1.1356ϫ10−3d 0͒,where d 0is the length of unit cell.Again,we find that the two boundaries of the effective homogeneous medium are close to the outer unit-cell boundaries of the 1D metamaterial.Figure 8shows the retrieved z ,n ,⑀,and for one cell of this metamaterial.It can be seen that the frequency range of 13.8GHz–14.5GHz is a LH band,which agrees with the conclusion in Ref.[14].It should be noted,however,that for a 2D version of this metamaterial,the effective boundaries should be obtained from the optimization process,as they do notnecessarilyFIG.6.S 11and S 21(real and imaginary parts )for three cells:Comparison between FDTD re-sults (dot line with *)and calculated S parameters based on the retrieved ⑀and (solid line with ᮀ)for a one-cell metamaterial shown in Fig.1.FIG.7.Retrieved ⑀and (real and imaginary parts )for a one-cell metamaterial shown in Fig.1.The vertical dashed lines denote the limits of the resonance band.ROBUST METHOD TO RETRIEVE THE CONSTITUTIVE …PHYSICAL REVIEW E 70,016608(2004)match the unit-cell boundaries of the metamaterial.Indeed,in this specific case,we obtain x ¯opt =͑0.33110d 0,0.30342d 0͒.III.CONCLUSIONWe have proposed an improved method to retrieve the effective parameters (index of refraction,impedance,permit-tivity,and permeability )of metamaterials from transmission and reflection data.The successful retrieval results for vari-ous metamaterial structures show the effectiveness of the method.Our main conclusions are as follows:(1)The first boundary and the thickness of the effective media can be determined by matching z through all sample frequencies for different lengths of the slabs in the propaga-tion direction.For symmetric 1D metamaterials,we have drawn the empirical conclusion that the first boundary coin-cides with the first boundary of the unit cell facing the inci-dent wave,and the thickness of the effective medium is ap-proximately equal to the number of unit cells multiplied by the length of a unit cell.For 2D and asymmetric 1D metama-terials,the effective boundaries have to be determined by optimization.(2)The requirement z Јജ0cannot be used directly for practical retrievals when z Јis close to zero because the nu-merical or measurement errors may flip the sign of z Ј,mak-ing the result unreliable.In this case,we have to determinethe sign of z by the value of its corresponding n so that n Љജ0.(3)There is a resonance band characterized by the fact that the requirement Љജ0and ⑀Љജ0cannot be satisfied at those frequencies.On each side of the resonance,the branch of n Јcan be obtained by a Taylor expansion approach con-sidering the fact that the refractive index n is a continuous function of frequency.Since the refractive index n at the initial frequency determines the values of n Јat the subse-quent frequencies,we determine the branch of the real part of n at the initial frequency by requiring that Љand ⑀Љare non-negative across all the frequency band.(4)Due to the noise contained in the S parameters,the retrieved n and z at some specific frequencies are not reli-able,especially for thicker metamaterials at lower frequen-cies.In spite of this,the fact that S 11and S 21for multiple cells of metamaterial calculated from the retrieved ⑀and for a unit-cell metamaterial match the S 11and S 21computed directly from numerical simulation confirms that the metamaterials can be treated as an effective homogeneous material.ACKNOWLEDGMENTSThis work was supported by DARPA (Contract No.N00014-03-1-0716)and ONR (Contact No.N00014-01-1-0713).[1]D.R.Smith,W.J.Padilla,D.C.Vier,S.C.Nemat-Nasser,and S.Schultz,Phys.Rev.Lett.84,4184(2000).[2]R.A.Shelby,D.R.Smith,and S.Schultz,Science 292,77(2001).[3]L.Ran,X.Zhang,K.Chen,T.M.Grzegorczyk,and J.A.Kong,Chin.Sci.Bull.48,1325(2003).[4]D.R.Smith,S.Schultz,P.Markoš,and C.M.Soukoulis,Phys.Rev.B 65,195104(2002).[5]P.Markošand C.M.Soukoulis,Opt.Express 11,649(2003).[6]R.W.Ziolkowshi,IEEE Trans.Antennas Propag.51,1516FIG.8.Retrieved z ,n ,⑀,and (real and imaginary parts )for a one-cell metamaterial structure taken from Refs.[14,21]and shown in the inset (a ).The vertical dashed lines denote the limits of the resonance band.CHEN et al.PHYSICAL REVIEW E 70,016608(2004)(2003).[7]A.M.Nicolson and G.F.Ross,IEEE Trans.Instrum.Meas.19,377(1970).[8]J.Baker-Jarvis,E.J.Vanzura,and W.Kissick,IEEE Trans.Microwave Theory Tech.38,1096(1990).[9]J.A.Kong,Electromagnetic Wave Theory(EMW,Cambridge,MA,2000).[10]J.A.Kong,Prog.Electromagn.Res.35,1(2002).[11]D. A.Shelby, D.R.Smith,S. C.Nemat-Nasser,and S.Schultz,Appl.Phys.Lett.78,489(2001).[12]C.D.Moss,T.M.Grzegorczyk,Y.Zhang,and J.A.Kong,Prog.Electromagn.Res.35,315(2002).[13]R.Storn and K.Price,J.Global Optim.11,341(1997).[14]T.M.Grzegorczyk,C.D.Moss,J.Lu,and J.A.Kong,NewRing Resonator for the Design of Left-Handed Materials at Microwave Frequencies,Proceedings of Progress in Electro-magnetics Research Symposium,Honolulu,Hawaii,2003 (EMW,Cambridge,MA,2003),p.286.[15]C.R.Simovski,P.A.Belov,and S.He,IEEE Trans.AntennasPropag.51,2582(2003).[16]J.E.Sipe and J.V.Lranendonk,Phys.Rev.A9,1806(1974).[17]P.A.Belov,S.A.Tretyakov,and A.J.Viitanen,Phys.Rev.E66,016608(2002).[18]T.Koschny,P.Markoš,D.R.Smith,and C.M.Soukoulis,Phys.Rev.E68,065602(2003).[19]ndau,E.M.Lifshitz,and L.P.Pitaevskii,Electrody-namics of Continuous Media,2nd ed.(Pergamon,Oxford, 1984).[20]V.Veselago,p.10,509(1968).[21]C.D.Moss,Ph.D.thesis,Massachusetts Institute of Technol-ogy,2004.ROBUST METHOD TO RETRIEVE THE CONSTITUTIVE…PHYSICAL REVIEW E70,016608(2004)。
ANSYS workbench核工业多相流动和传热
― 用户扩展到几乎所有的蒸汽发电公司
― 所有主要的OEM都是ANSYS的用户
― 几乎所有的核电厂配套设施主要设备的开发过程中都使用过
ANSYS软件
― 美国及其它的监管机构都是用ANSYS软件
5
© 2017 ANSYS, Inc.
August 3, 2017
ANSYS UGM 2017
1
ANSYS 与核工业
Facility located at IRSN Saclay
Volume: 7 m3, height: 4.8 m Wall temperature controlled Test 101: depressurization test
21
© 2017 ANSYS, Inc.
August 3, 2017
22
© 2017 ANSYS, Inc.
August 3, 2017
ANSYS UGM 2017
硼稀释
• 欧盟的EUBORA和FLOWMIX-R项目使用ANSYS Fluent产品 • ANSYS-CFX被应用在VVER-1000机组的相关分析中 • 该分析需要模拟完整的一回路系统，要求能够合理模拟硼在冷
同时能够直接打开和编辑其它商业软件的几何文件网格网格划分难度大需要灵活的网格划分策略turbogrid短时间内对复杂的叶片和叶栅通道自动生成高质量的结构化网格fluentmeshing使用网格包裹技术针对大规模复杂网格生成高质量的网格文件求解器计算规模大流场复杂多数情况下涉及到多物理场耦合ansysworkbench协同仿真平台能够快速完成不同软件之间的数据传递实现多物理场间的单向和双向耦合11august32017ansysugm2017?2017ansysinc
ANSYS Fluent Overset Mesh应用介绍
12
重叠最小化
基于donor优先的网格体积单元 − 网格过渡位置的网格分辨率是一致的 − 简化网格设计 − 如果非结构网格尺寸分布不均匀或者
网格尺寸增长变化不是单调变化的， 那么创建的重叠网格交界面形状是不 规则的
14
重叠最小化
基于donor 优先的边界距离 − 网格过渡尽可能的远离边界；最好在两
边界中间
− 因为网格的过渡不考虑网格的分辨率，
所以需要更加仔细的进行网格设计
15
重叠最小化
网格优先顺序 − 按照设定的优先级顺序生成重叠网格
对背景网格和组件网格是有利的
• 包含局部供体单元的网格优先级高
− 网格过渡不考虑网格的分辨率
0
0
0 1
1
1
0
0
16
重叠拓扑
组件网格可任意重叠 − 重叠网格边界是允许的
solve
donor dead
receptor
Background grid
Component grid
10
域的连通性
为了避免出现orphan cells，就需要有足 够多的网格重叠 − 边界临近问题 − 为了避免孤立网格单元，重叠区域最少要有4
个网格单元
wall
To display Receptor Cells:
（donor），查找的donor cells是用来为插值提供数据信息的。Donor cells是solve cells的子集。Donor cells查找失败的receptor cells被标记为 orphan cells。
9
域的连通性
重叠单元类型 − Solve cells [1] − Receptor cells [0] − Donor cells [2] − Orphan cells [-1] − Dead cell [-2]
网格与CFD求解精度的关系
• Hex, Tri, Quad: 大多数网格的EquiAngle Skew 应该低于0.85。 • Tetrahedral: 大多数网格的EquiAngle Skew 应该低于0.9， 对于某些简单物理问题可以更大。 对于关键区域网格的Size Change应该低于2， 实际尺寸变化要视所研究的物理现象而定。
R
c
Circumscribed circle
b
e2
ac 2 bd 2
21
网格质量 – Diagonal Ratio
Diagonal Ratio 指标 (QDR) 仅用于Quad和Hex网格，定义如下:
QDR max d1, d 2 , , d N min d1, d 2 , , d N
c1 c0
node值是cell值的加权平均
7
FLUENT中的有限体积法(6)
面值计算
Green-Gauss
FLUENT默认 相邻网格中心点的代数平均 计算量较小，准确度可以接受
Cell-Based method
Green-Gauss
Node-Based method
所求解的面上各节点值的代数平均 节点值取作周围网格单元值的代数加权平均 当非结构化网格时，比Cell-Based 更准确
Coarse mesh (228 cells)
medium mesh (912 cells)
ansys软件问答合集(二)
47 在Ansys中，碰到提示“Volume 1 cannot be meshed. 208 location(s) found where non-adjacent boundary triangles touch. Geometry configuration may not be valid or smaller element size definition may be required.”。这是什么问题？ 回答：提示就是告诉你需要更小的单元，可能单元太大的时候出现的网格有有问题，比如狭长 的网格，计算的时候集中应力太大。
48 在Ansys中，碰到错误Volume11 could not be swept because a source and a target area could not be determined automatically。please try again...，这是什么原因？ 回答：体不符合SWEEP的条件，把体修改成比较规则的形状，可以分割试试。 49 在Ansys中，碰到警告和错误：“*** WARNING *** SUPPRESSED MESSAGE CP = 1312.641 TIME= 16:51:48 An error has occurred writing to the file = 12 which may imply a full disk. The system I/O error = 28. Please refer to your system documentation on I/O errors. ”，这是什 么错误和警告？ 回答：1.I/O设备口错误，I/O=26，错误，告诉你磁盘已满，让你清理磁盘。但是实际问题的解 决不是这样，是你的磁盘格式不对，将你的磁盘格式从FAT26改称NTFS的就可以了。因为 FAT26格式的要求你的单一文件不能大于4G。但是我们一旦做瞬态或者是谐相应的时候都很 容易超过这个数，所以系统抱错。Байду номын сангаас2.I/O设备口错误，I/O=9,错误，和上一个一样告诉你磁盘已满，让你清理磁盘。但是实际问题 是由于你的磁盘太碎了造成的，你只要进行磁盘碎片整理就可以了，这个问题就迎刃而解。
5.4多相流模型5.5分散相模型2
4 α p ρppVol α p ρ ppVol 2t n 1 n 1 n 1 n 1 Anb nb p SU Sp p
n 1 n
3 α p ρppVol
n 1
(5.383)
将上式重写为
App Anbnb S
(2) 体积分数方程
a. 体积分数方程 通过求解一相或多相体积分数的连续性方程，可以追踪各相之间的界面。第 q 相体 积分数的连续性方程为
n 1 pq m qp (5.385) α ρ α ρ v S m q q q q q α q ρq t p 1 qp 为从相q向相p的传质； m pq 为从相 其中，ρq为第q相的物理密度；vq 为第q相的速度； m
在 VOF 模型中，各相流体共享一个方程组，每一相的体积分数在整个计算域内被追踪。 适用 VOF 模型的多相流应用包括分层流、有自由表面流动、液体灌注、容器内液体振 荡、液体中大气泡运动、堰流、喷注破碎的预测和气 - 液界面的稳态与瞬态追踪等。
114
沈阳航空工业学院
(2) 混合模型
混合模型的相可以是流体或颗粒，并被看作互相穿插的连续统一体。混合模型求解 混合物动量方程，以设定的相对速度描述弥散相。适用混合模型的应用包括低载粉率的 带粉气流、含气泡流、沉降过程和旋风分离器等。混合模型还可以用于模拟无相对速度 的匀质弥散多相流。
117
沈阳航空工业学院
则有下列 3 中可能：
αq = 0：单元中没有第q相流体。 αq = 1：单元中充满第q相流体。 0 < αq < 1：单元中有第q相流体与其它一相或多相流体间的界面。
根据局部αq值，计算域内每一控制容积被赋予适当的物性和变量值。
重力异常二维正演中的无网格方法
重力异常二维正演中的无网格方法李俊杰;严家斌【摘要】无网格法是一类新型数值算法,具有精度高、高阶形函数构造与物性加载便利等特点,在计算力学领域应用广泛.将无网格方法(PIM、RPIM及EFGM)用于重力异常场二维正演计算:首先从重力异常二维变分问题出发,利用Galerkin法结合高斯积分公式推导了对应的无网格离散系统矩阵表达式;其次通过数值试验得出了RPIM-MQ、RPIM-exp及EFGM-exp形状参数的建议值,最后比较分析了最优形状参数下不同无网格法的计算效果.结果表明:无网格法适用于介质物性分布变化较大的重力异常二维正演,exp函数形状参数αc最优取值区间为[1.5,1.7],β建议值为0.6,MQ函数g取值区间为-4.1～1.9;EFGM较PIM及RPIM具有更高的计算精度.【期刊名称】《煤田地质与勘探》【年(卷),期】2018(046)006【总页数】6页(P181-186)【关键词】无网格法;点插值法;径向基点插值法;无单元Galerkin法;重力异常【作者】李俊杰;严家斌【作者单位】浙江省水利水电勘测设计院,浙江杭州310002;中南大学地球科学与信息物理学院,湖南长沙410083【正文语种】中文【中图分类】P631重力异常二维正演一般采用解析积分法，然而对于非均匀、几何形状复杂的地质体需要利用数值积分进行近似的正演计算，常用的方法是基于三角剖分的有限元法[1-4]，但高质量的三角剖分程序实现过程较复杂。
无网格法[5]作为一类新型节点数值算法，虽然其发展历史仅20余年，但其具有形函数构造无需网格，精度高，计算一般复杂模型及连续介质模型较网格方法便利等优势，引起了国内外学者浓厚的兴趣，现已成为数值计算领域一大研究热点。
无网格法在弹性波场[6-8]、雷达波场[9]、大地电磁场[10-17]、直流电场[18]及磁场[19]的正演模拟中取得了一定的成效，研究结果表明当地下介质为一维时，数值方法能取得很高的计算精度[10-11]，但对应的高维问题不存在解析解，计算复杂模型时只能通过断面图异常的形态来近似反映算法的优劣。
超级计算机
¾ 超级计算机...............................................................................................................................1 超级计算环境 2007 年 3 季度运行情况简报 .........................................................................1
深腾6800
● 共197名用户，3季度增加用户6名。 ● 有134名用户利用LSF提交作业，共完成.51000多个作业，用户作业平均规模为5.9个CPU，累计 使用机时112万CPU小时（按Walltime计算）。 ● 2007年3季度，深腾6800的磁盘阵列系统与QsNet网络系统先后发生故障，导致深腾6800的平均 整体使用率有所下降，为83.5%（按Walltime计算），平均CPU利用率69.1%（按CPUtime计算）。CPUtim e与Walltime之比平均为82.7%。 ● 2007年3季度，作业平均等待时间为23.3小时。 ● 已完成作业按规模分布情况：串行作业数量占62.6%，4处理器节点内并行作业数量占21.1%。 而根据作业使用的CPU小时计算，占用机时最多的并行作业规模分别为16处理器、33-63处理器、32处理 器、64处理器，其比例分别为21.8%，18.1%，17.9%和12.4%，串行作业仅使用总机时的2.1%，表明深腾6 800的计算机时还是主要用于较大规模的并行作业计算。
1. Introduction ...................................................................................................................2
经验总结
附件里的文件中是6月份在网上下载的几个作业题，我自己做了一下，曾经应一个网友求助给过他。这些例子是介绍精馏塔从简洁计算到严格计算，再关联vinyl chloride反应及吸收的过程，这些题目的设计其实就是引导学生大概领略一下几个计算领域的情况。今天拿出来原封给大家，然后我就在我的电脑上删除了。文件的.txt改成.rar，解压。
第二是电解液的选择使用了重复的设定和求解，导致SOLVER01/02互为条件反复求解。我将在FLOWSHEETING OPTION中删除了这些指标。
另外，不要将没有使用的物料的选择也选择上。
有时候一些人设置了一些求解的约束条件，但是有时候这些条件互为条件或者矛盾。在模拟设计中很忌讳一个条件控制两个以上的有因果关联的参数。
SRKK方法也同样适用于水－烃类系统。
2.PROII
原题：PROII出现如下错误，请教高人是什么原因？
2 error(s), 0 warning(s), and 0 message(s) detected:
*** error *** Invalid specification in unit 1, 't1' - "primary" section refers to a non-existent stream ''.
谁知道这个计算公式．．
我说一下自己的看法，不知道是不是准确，供大家参考。
我们首先要知道为什么找到一个精馏塔的最佳回流比？就是象找到在实际的塔操作中用准确的回流量来控制塔顶的产品质量。如果不是这样的话大家谁会要这样做呢，你说是不是。但是我劝大家不要这样做，这样的命题只有在专业论文里面才有用。
为什么这样说呢。因为在实际的塔顶产品质量控制算法中，为了简化运算都要加入人工外来数据干预。就是说采用在线分析仪表实时的分析塔顶组分，然后数据实时被运算程序采集，作为运算程序的纠偏。如果没有这个，仅仅凭借软件自己来无限级别的回归运算，有跑偏不说，只有核爆炸模拟才那么做吧（我猜测）。我在美国学习装置模拟优化的时候，美国人不这么做。如果凭借完全的软件来沙盘推演的话，大概过程应该这样（以PRO_II为例）：在建立到精确模型计算后，在模块中建立CASE STUDY,建立在线分析（你可以建立无数个），然后数据反馈给ASPEN模型，至运行OK。然后启用动态模型（就像本坛紫色月亮的Invensys Simsci Dynsim），在动态中找到最佳回流比。用ASPEN也是一样过程。你还愿意做吗？
ANSYS R19.0 FLUENT 多相流改进
正位移泵中的空化
• 业界长期使用ANSYS/FLUENT来模拟正位移泵 • 客户使用特定的齿轮，齿轮转子和活塞泵模型进行气蚀测试。 • 随着R19空化模型的改进，齿轮泵和活塞泵的流量/角速度曲线和
填充速度预测得到了显著提高。
Discharge
Fill speed
Angular velocity
泄露
相间面密度框架的改进
- 相间面积密度可用于mixture/drift-flux模型 - 之前只有欧拉多相流模型可用
- 适用于所有的相对组合 - 之前只适用于连续-分散系统
- 引入了基于相间面密度的梯度方法
欧拉多相流
14
© 2018 ANSYS, Inc.
April 25, 2018
Mixture 多相流/drift-flux 模型
毛细管压 力影响
水相速度矢量
无毛细管压力
26
© 2018 ANSYS, Inc.
April 25, 2018
ANSYS 新品发布会和技术巡展 2018
毛细管压力将水向反重力方向拉动
多相流-非迭代时间推进法(NITA)
• 多相流模型 NITA • 算法改进（R18.x) • 升力和壁面润滑力的鲁棒性改进（R19，隐藏） • 在动量子迭代间相对速度的显式处理 • 默认启用 • 可用（ rpsetvar ‘mp/nita/mom-force-expl? #f ）禁用此功能 • 体积力加权的压力插值作为TUI的选项（R19） • 用 solve>set>multiphase-numerics/nita-controls/face-pressure-options调用
3
© 2018 ANSYS, Inc.