Home

ZINC Tutorial Manual

1. 29 2 vscan eps GSview Eile Edit Options View Orientation Media Help 1 Scan distance m 1 Scan distance m 25 Dzmap eps GSview Ele Edit Options View Orientation Media Help ee i Garin CL QQ EE Dz C m 2 3e 012 3 5e 012 4e 012 4 5e 012 5e 012 5 5e 012 6e 012 6 5e 012 7e 012 30 Fig 30 Encapsulated postscript graphs generated by Zinc using key plot 1 Ifyou don t like gnuplot You don t have to use it Since the data being plotted is created and stored as out files you can plot it using any spreadsheet program or graph plotter Or you could use scripting languages like Matlab to read the out files and prepare graphs It s also possible to get Matlab to write Zinc s input files call Zinc and read in and plot the resulting graphs Also in the directory Dielslap zin Exactly the same simulation but with key sim 3 That is diagonally scaled GMRES solver Example 2 piezoelectric example Directory examples piezo Introduction We now consider a piezoelectric problem Here there are 4 variables to solve for the x y z components of elastic displacement ux uy uz and the electrostatic potential V This is our first multiphysics problem However note that the problem is linear Z ZMesh 31 Fig 1 Piezoelectric simulation consisting of two piezoelectrics cube in a cube Region 3 nodes cyan are clamped and fixed at 1V Region 4
2. SB21 23 1 3 B21 21 2 1 SB22 2 2 2 2 SB22 2 3 2 3 SB22 col 3 row 1 Q block 11 3 1 SOL 1 2 3 2 BD 1 3 3 3 Q1 col 3 row 2 Q block 21 3 1 502 2 2 3 2 502 2 3 3 3 Q2 E section at bottom right 3 1 3 I SE 3 20 De SE 3 3 3 3 SE region 1 elements 2 values f 1 Sf1 2 Sf2 L EAE AE AE AE FE AE AE AE AE E FE AE AE AE AE E FE AE AE AE AE AE FE AE AE AE AE AE FE FE AE AE AE AE FE FE AE E AE AE E FE AE AE AE AE E FE AE AE AE AE AE FE AE AE AE AE AE FE AE AE AE AE AE FE HE FE HE region 2 nodes 3 values 1 xic 2 x2c 3 pc region 3 nodes 1 value 3 pref 50 1 xic 2 x2c 3 pctpref 2 The main concept here is token replacement whereby elements in the C and f matrices are defined not by expression but by tokens beginning Any string can be used when Zinc comes across the string it calls a user supplied function with the token name passed in as an argument Thus for the C functions Zinc calls function cfun for the f tokens it calls ffun Here is the cfun function written in Fortran cathodetoken f90 function Cfun label x y z ur dur nvar istep use nlinc implicit none character label integer nvar istep double precision Cfun x y z ur nvar dur nvar 3 double precision D11 D12 D13 D21 D22 D23 D31 D32 D33 double precision x1 x2 p x3 wl w2 w3 m denom call NLinit x1 ur x2 ur p ur 1 2 3 x3 1 x1 x2 m ml x1 m2 x2 m3 x3
3. dielectric05 out lt lt linescan100 0 0 00 10 00 00 1 0 4 1 eps1 2 eps2 vz dielectric6 out lt lt planescan 200 1 0 1 0 1 0 10 40 40 v dielectrict out lt lt planescan 2 0 0 1 0 1 0 1 0 1 0 40 40 41 eps1 2 eps2rvz Status READY Open Run Cancel Fig 24 The gnuplot script file used to generate the graphs The gnuplot file consists of a series of commands most of which revolve around the plot command There are also lots of set commands In this case we are setting the term inal to be emf so that gnuplot will output enhanced metafiles Let s consider the first plot command in this file set output d mgc zincprog examples dielectric dielectric0l emf set title 0 10100E 00 0 30000E 00 0 10000E 01 to 0 10100E 00 0 30000E 00 0 10000E 01 set xlabel scan distance set ylabel v plot d mgc zincprog examples dielectric dielectricOl out u 1 5 wl This instructs gnuplot to plot the file dielectricO1 out to make graph file dielectricOl emf u 1 5 means use columns 1 and 5 resp for the x and y axes w 1 means with lines ie datapoints are connected using lines For more information about the gnuplot commands see the gnuplot manual and tutorials However a lot can be gleaned from the default gnuplot script file generated by zppwin So how do we improve the graph quality All gnuplot does it plot the out files Since these file
4. 1 epsl 2 eps2 Vz ll Q The 1 eps1 2 eps2 expressions means if we are in region 1 use value eps1 if we are in region 2 use value eps2 Here eps and eps2 are parameters appearing in the con file eps0 8 854e 12 epsl eps0 eps2 eps0 2 As we scan in the plane indicated we will encounter region 1 and region 2 inner sphere areas The appropriate value of permittivity needs to be taken in each case The expression is then eps Vz Vz means dV dz This is the standard notation used by Zinc since we have given potential the name V the derivatives can be accessed using Vx Vy Vz Thus the whole expression is eps Ez Dz the z component of displacement field The parameters after the planescan command indicate what sort of scan is needed 2 means perpendicular to the y direction 1 0 1 0 means scan x between 1 0 1 0 The next 1 0 1 0 means z 1 0 1 0 40 40 is the number of points requested along each direction 23 The meaning of the other scans should now be apparent Note that quite arbitrary functions involving the variables solved for and the constants in con can be used The second line scan plots V 2 which is the potential squared All the usual operators are allowed and the functions sin cos tan exp log etc can be used also see the Zinc User Manual Chapter 5 for more information Traceability Before continuing this tutorial a quick word about traceability In running the
5. 3e 013 0 10000E 01 0 10000E 02 0 40000E 02 to 0 10000E 01 0 10000E 02 0 40000E 02 0 0 005 u variation x scan UZ 1e 012 1 1e 012 1 2e 012 1 3e 012 1 4e 012 1 5e 012 1 6e 012 1 7e 012 0 01 0 015 0 02 scan distance 0 10000E 01 0 10000E 02 0 40000E 02 to 0 10000E 01 0 10000E 02 0 40000E 02 0 0 005 0 01 0 015 0 02 scan distance 38 Voltage variation x scan 0 10000E 01 0 10000E 02 0 40000E 02 to 0 10000E 01 0 10000E 02 0 40000E 02 0 69765 0 6976 0 69755 0 6975 0 69745 gt 0 6974 0 69735 0 6973 0 69725 0 6972 0 69715 E variation x scan 0 005 0 01 scan distance 0 015 0 10000E 01 0 10000E 02 0 40000E 02 to 0 10000E 01 0 10000E 02 0 40000E 02 49 15 49 2 49 25 49 3 VZ 49 35 49 4 49 45 49 5 U Variation z scan 0 005 0 01 scan distance 39 0 015 0 02 5e 012 4 5e 012 4e 012 3 5e 012 3e 012 2 5e 012 UX 2e 012 1 5e 012 1e 012 5e 013 5e 013 0 50000E 02 0 00000E 00 0 10000E 01 to 0 50000E 02 0 00000E 00 0 10000E 01 0 0 005 uz Variation Z scan 2 5e 012 2e 012 1 5e 012 1e 012 5e 013 UZ o 5e 013 1e 012 1 5e 012 2e 012 2 5e 012 0 01 scan distance 0 015 0 02 0 50000E 02 0 00000E 00 0 10000E 01 to 0 50000E 02 0 00000E 0
6. 4 Initial conditions parameters Oi 0 5 Mass fraction of H in the anode Oo 0 5 Mass fraction of H20 in the anode O2 1 3 Mass fraction of O gt in the cathode Wino 1 3 Mass fraction of H20 in the cathode WO 1 3 Mass fraction of N gt in the cathode xe 0 9 Molar fraction of H in the anode 0 1 Molar fraction of H20 in the anode xd 0 2550607 Molar fraction of O2 in the cathode oe 0 4534413 Molar fraction of H20 in the cathode X2 0 2914980 Molar fraction of N in the cathode p 1 013x10 Pressure in the anode Pa p 1 013x10 Pressure in the cathode Pa 47 Zinc implementation We now need to see how to implement equation 4 in Zinc The Zinc PDE specification is V CVu au f 5 n CVu gu g We will set a q g 0 and put non zero values only in the C and f matrices Let us split the C matrix into 3 parts 1 C 2 A COA 3 Then comparison of 5 with 4 shows that B 0 0 B O0 0 Q0 0 C 0 B 0 0 B 0 0 Q ol 6 0 0 B 0 0 B 0 0 Q Ce B 0 0 B 0 0 Q 0f 7 0 B 0 0 B 0 0 0 000000 E 0 0 C 0 00000 0 E ol 8 000000 0 0 E with Tn X r Br OD jik 1 n 1 9 n 1 O 5 Y ALC k 1 n 1 10 jl pores 11 RT n with A and C defined as 48 m x m m Ay RR LR m 12 m x m C D KER pt 13 KR m Sa This completes the specification of the C matrix The f matrix is given by L f f 0 with c PE I
7. C55 C54 C53 15 25 35 C41 46 Cas C46 C42 C44 C45 C44 C43 14 24 34 C31 C36 C35 C36 C32 C34 C35 C34 C33 13 23 33 11 16 15 16 12 14 15 14 13 Kit A12 A13 21 26 25 26 22 24 25 24 23 7K21 k22 K23 31 36 35 36 32 C34 35 34 33 7K31 kK3 k33 where c etc are the elastic stiffnesses N m2 e1 etc are the piezoelectric constants C m2 K etc are the permittivities F m The upper left block of this matrix gives the elastic behaviour the most complicated part of the simulation the lower right block gives the dielectric behaviour and is essentially the same as the matrix we derived in Example 1 except we allow for anisotropy The off diagonal blocks specify the elastic electric coupling See User Manual chapter 4 for the full derivation of this matrix We will use C for the Zinc matrix and c for elastic stiffnesses Although the two are related it s important not to get confused In our problem all the other matrices can be set to zero as we have no body surface forces or charges If these were required the f matrix would be used for body forces and charges while the g matrix would be used for surface forces and charges In the case of steady state elastic vibration the a matrix can be used to encode the material densities and vibration frequency in fact the a vector would be given as 32 a diag p
8. TD 7073 1 20 0 101 0 3 1 0 V linescan 100 0 101 0 3 1 0 Qu LOG 0 3 FL0 V 2 linescan 100 0 101 0 3 1 0 0101 2013 1 0 Vz linescan 100 0 0 O20 L 0 0 0 0 0 1 0 Vz linescan 100 0 0 0 20 1150 0 0 0 0 1 0 amp 1 epsl 2 eps2 Vz yplane y x1 x2 z1 z2 Nx Nz planescan 2 10 50 1 0 1 0 1 0 T0 40 40 lt V planescan 2 020 1 0 1 0 1 0 1 0 40 40 amp 1 epsl 2 eps2 Vz Here we are requesting 7 scans 5 linescans and 2 planescans Linescans are cuts through the data along a line The first one is a line scan with 100 points between coordinates 0 101 0 3 1 0 and 0 101 0 3 1 0 i e it is a scan along the z direction The function to be plotted is shown after the equals sign V in the case the electrostatic potential Zppwin knows about V because it also reads the zin input file in which we recall we put the command Labels V The labels command is optional if you did not use it you d have to refer to ul instead of V In multi variable simulations the default names are ul u2 etc Clearly it is more readable to specify a mnemonic name using labels The file illustrates the use of comments with and the continuation character amp for splitting lines which are uncomfortably long Whether split or not each scan command should not be more than 1000 characters long Having prepared the zpp file launch ZPP and select the file using the Open button as shown in Fig 19 19
9. fit that category See the file cathode zin for the full set of expressions Running the simulation After a lengthy preamble we are now ready to run the simulation Fig 2 shows the simulation geometry which may be compared to Fig 1 The geometry is essentially 2 D so we have only two planes of elements in the third direction Note that the distances in Fig 2 are in mm whereas the equations we have derived are in SI units m therefore a scale factor of 0 001 is needed and this is in the command scale le 3 in the cathode zin files Finally let s look at the node fixing Dirichlet boundary part of the ZIN file and the initial condition ion 2 nodes 3 values x1c x2c pc e oi i Q r 1 2 3 region 3 nodes 1 value 3 pref init 1 xic 2 x2c 3 pctpref 2 This says that region 2 nodes will be completely fixed see Fig 2 but region 3 nodes will be fixed only in terms of the pressure The two mass fractions are allowed to vary The initialisation block init fixes the mass fractions and pressures at x1c x2c and pc pref 2 These constants are all in cathode con 54 Z ZMesh File View Help 6 56606E 61 8 z 0 50000E701 Fig 2 Fuel cell geometry Mass fractions and pressures are fixed on region 2 nodes green Pressure only is fixed on region 3 nodes cyan The Zinc window is shown in Fig 3 Note that Zinc has noticed the presence of the cathodetoken f90 file and allows you to view it
10. gt Open ZOU option as shown in Fig 17 We can see from the figure that the boundary conditions V 1 and V 0 Volts have been implemented In systems with more than one variable you can switch between the different variables using buttons at the top left of the screen The usual features to hide regions and take slices are still available Fig 18 shows a slice through the data Note that the inner region region 2 has been set to show the region while the outer region is showing the data Thus the colour key elements on the right of the screen now have three states on DATA off 17 Z ZMesh File View Help Showing variable 1 of 1 Variable name U 1 6068 Fig 17 Viewing the simulation output ZOU file Z iMesh Eile View Help Showing variable 1 of 1 Variable name U 1 6668 Fig 18 Slice view through the data showing simulation result in region 1 Post processing We now wish to examine the output of the simulation using ZPP the Zinc post processor This program allows us to take various scans through the data and prepare graphs either of the simulated quantities electric potential in this case or quantities derived from these In a particular we will show how to calculate the E field and D field 18 Open a new text file called zpp Zinc post processor and type in the following key plot 2 output emf files ftol le 3 tolerance for fitting N x1 yl z1 x2 y2 z2 linescan 100 K
11. in region I and c12inc inclusion when we are in region 2 Note that uxx is simply du dx which is the necessary strain etc The last two linescans make use of constants in piezo con eps0 8 854e 12 c11 12 6e10 c33 11 7e10 c44 2 3e10 c12 7 95e10 c13 8 41e10 c66 2 347e10 e31 6 5 e33 23 3 el5 17 k11 1700 eps0 k33 1470 eps0 c22 c11 c23 c13 c55 c44 e32 e31 e24 e15 k22 k11 cllinc 12 6e11 c33inc 11 7e11 c44inc 2 3e11 36 cl2inc 7 95e11 cl3inc 8 4lell c66inc 2 347el1 e3linc 6 5 e33inc 23 3 el5inc 17 k1llinc 1700 eps0 k33inc 1470 eps0 c22inc cllinc c23inc cl3inc c55inc c44inc e32inc e3linc e24inc el5inc k22inc kllinc Note that parameters can be defined in terms of previously specified parameters Here we use this to enforce the symmetry equivalence of certain coefficients assuming the piezoelectric has 6mm symmetry and is poled along the z direction Results The following are the enhanced metafiles EMF supplied with the Zinc distribution for this example piezo01 emf piezo02 emf etc You should check you get the same results U Variation x scan 0 10000E 01 0 10000E 02 0 40000E 02 to 0 10000E 01 0 10000E 02 0 40000E 02 1 5e 011 T 1e 011 5e 012 5e 012 1e 011 1 5e 011 T l 0 0 005 0 01 0 015 0 02 scan distance 37 uy Variation x scan UY 1 2e 012 1 1e 012 1e 012 9e 013 8e 013 7e 013 6e 013 5e 013 4e 013
12. problem we will consider a cathode simulation as shown in Fig 1 We will take n 3 species namely O2 H2O and No Since we are dealing with mass fractions which sum to unity we need only consider the first two species Thus there are three variables to be solved Xo2 XH20 p K p yo yc p U 0237 H20 Key 2 mm Cathode 1 mm Insulation Inf low Outflow Pret gt gt 0 25 mm Figure 1 Geometry and boundary conditions for the cathode model The D coefficients in 4 are not constants but themselves depend on the mass fractions as 2 2 0 0 Q 4 QA D XD xD X3Dp 11 X X X EE E a D lt 22 X X X D0 2 3 DoD DD D D 2 2 2 o U UU D xD xD xD 33 X X X L 4 2 T 2 alo o olto o D xD xD xD Et rn DoD D D D D olo o O 0 0 _ 0 D x D xD XD ETE EE D D D D 13 23 0 Sz o ma 0 o o Due xD xD x D EE EE The D coefficient matrix is symmetrical Tables 1 4 summarise the full set of material properties used in the simulation These are present in the constants file cathode con cathodetoken con Table 1 Constants for the electrode flow models Name Value Description Units R 8 314 Gas constant J K mol T 353 Temperature K KH2 3 9x10 Henry s constant for H2 dissolved in H2O Pa m3 mol KO2 3 2x10 Henr
13. they are not an error the output GO GGG GGG 544522E 03 544522E 03 544522E 03 544522E 03 544522E 03 544522E 03 544522E 03 544522E 03 544522E 03 544522E 03 544522E 03 544522E 03 544522E 03 544522E 03 544522E 03 544522E 03 544522E 03 544522E 03 544522E 03 9 26543264238894999 E 12 9 26543264238166977E 612 3 G GGG GGO GO GGO GOGO GO G 147638E 14 1429 3E 14 139363E 14 135537E 14 131926E 14 128268E 14 124466E 14 126882E 14 117138E 14 113919E 14 116593E 14 187672E 14 164774E 14 161661E 14 985463E 15 955720E 15 927607E 15 902031E 15 8 76 798E 15 convergence n in the first column Here we are running ond and fifth columns show numbers which te solution has been found to the problem d 4 columns show a number which approximates the left hand he set of equations to be solved Thes other Clearly in this simulation a accuracy The final energy of the system is hese numbers should be very similar to each has occurred Once the simulation has completed the Status shown returns to READY and the Zinc list file zls becomes viewable Fig 17 This file contains information about Zinc s input reflecting the specifications in the zin file Note that all the components of the Zinc matrices c a f g q are listed in this file including the zero ones so the user can check what was actually set up recall that we only set the c value
14. you should check you get the same answers we did It s probably a good idea to copy the input files to a fresh folder so as not the overwrite the output files we ve supplied That way you can do a comparison easily Example 1 Dielectric Problem Directory examples dielectric Introduction ZINC is a finite element program that is designed to solve problems having the form V CVU aU f 1 which is written explicitly in component form as follows a U Bean A EE ETE 2 Ox Ox where summation over 1 2 N where N is the number of variables is implied by repeated Greek superscripts and where summation over 1 2 3 is implied by Arabic suffices Appendix A gives a more detailed description of the notation and definitions The reader can also refer to the Zinc User Manual Chapter 4 for a more detailed description In this tutorial the problem is to determine the electrical potential distribution and other derived properties such as the electric field and displacement field The system consists of a spherical sample of one type of material that is embedded in a cuboid of another material The system is excited by a potential difference between opposing faces of the cuboid This problem is in fact mathematically eguivalent to both thermal and electrical conduction problems and to diffusion and magnetic permeability problems The field eguations that must be solved is one of Maxwell s eguations in the absence of d
15. zin Meshfile piezo mtf nstep 10000 Constants filepiezo con omega 0 1 i labels ux uy uz V NL file No NL file gt Restartfile None initial condition supplied region 1 elements 144 values C 1 1 0 12600E 12 c11 2 0 00000E 00 c16 3 Q 00000E 00 c15 1 000000E 00 c61 2 023470E411 cb6 3 Q 00000E 00 c65 1 0 00000E 00 c51 2 0 00000E 00 c56 3 0 23000E 11 c55 Output files Outputfile piezo zou List file piezo zls Status READY Open Run Fig 2 Zinc window ready for simulation od EO 0303 NOI nd al Having set up the geometry and output the MTF file we are now ready to run the simulation as shown in Fig 2 Having run the simulation we can take a look at the output using Zmesh Fig 3 Z ZMesh File View Help Showing variable 1 of 4 EP Variable name ux 9 13432E 19 9 13432E 19 Fig 3 Showing the output using Zmesh The Yellow and buttons can be used to cycle between the four variables solved for ux Uy uz V 34 WE now run Zpp to output scans through the data Fig 4 Zinc 2 0 Zpp Path DAMGCyzincprogjexamplesypiezc Inputfiles Zpp piezo zpp Sol n piezo zou Zinc inputfiles linescan 100 Input piezo zin linescan 100 Mesh piezo mtf linescan 100 Constants piezo con linescan 100 Gnuplot piezo gnu 3 Graph files Status READY Open Bun Fig 4 Post Processing with ZPP Let s look at the piezo zpp input file here k
16. 0 0 10000E 01 0 0 005 0 01 scan distance 40 0 015 0 02 Oxx X Scan 0 10000E 01 0 10000E 02 0 40000E 02 to 0 10000E 01 0 10000E 02 0 40000E 02 140 T T T E31INC VZ E31 2 C13INC UZZ 1 C13 2 C12INC UYY 1 C12 2 C11 INC UXX 1 20 L L 0 005 0 01 0 015 0 02 C11 2 o 1 scan distance Oyy X scan 0 10000E 01 0 10000E 02 0 40000E 02 to 0 10000E 01 0 10000E 02 0 40000E 02 150 T T E32INC VZ E32 2 100 fF I C23INC UZZ H1 C23 2 a S T L C22INC UYY 1 C22 2 C12INCPUXX 1 100 1 1 0 0 005 0 01 0 015 0 02 C12 2 scan distance n 41 Comparison graphs To check solution accuracy it is often a good idea to compare simulation results between different FE solvers We have therefore compared u z and u z with Comsol another finite element solver The Comsol results are in files ux txt and uz txt in the Zinc distribution The following gnuplot script file was written to perform the comparison set term emf set output uzcomp emf set xlabel z cm set ylabel uz pm plot piezo09 out u 1 100 5 1e12 t Zinc uz txt i 0 u 1 2 1e12 t Comsol w 1 set output uxcomp emf set xlabel z cm set ylabel ux pm plot piezo08 out u 1 100 5 1e12 t Zinc ux txt i 0 u 1 2 1e12 t Comsol w 1 This results in the foll
17. 12500E 03 0 99000E 03 0 00000E 00 0 1344E Zpp cathodetoken zpp 0 19800E 04 0 12500E 03 0 97020E 03 0 00000E 00 0 13446 0 39600E 04 0 12500E 03 0 95040E 03 0 00000E 00 0 13446 0 59400E 04 0 12500E 03 0 93060E 03 0 00000E 00 0 13446 0 79200E 04 0 12500E 03 0 91080E 03 0 00000E 00 0 13446 0 99000E 04 0 12500E 03 0 89100E 03 0 00000E 00 0 13446 Input cathodetoken zin 0 11880E 03 0 12500E 03 0 87120E 03 0 00000E 00 0 13446 Mesh cathodetoken mtf 0 13860E 03 0 12500E 03 0 85140E 03 0 00000E 00 0 13446 Constants cathodetoken con 0 15840E 03 0 12500E 03 0 83160E 03 0 00000E 00 0 13446 0 17820E 03 0 12500E 03 0 81180E 03 0 00000E 00 0 13446 4 Gnuplot cathodetoken gnu lt gt Graph files cathodetoken01 out lt lt linescan 100 0 125 0 99 0 0 0 125 0 99 0 0 x1 cathodetoken02 out lt lt linescan 100 0 125 0 99 0 0 0 125 0 99 0 0 lt 2 cathodetokenD3 out lt lt linescan 100 0 125 0 99 0 0 0125099 0 0 lt p cathodetoken04 out lt lt planescan 30 02 0 1250 125 0 99 0 99 20 50 x1 cathodetoken05 out lt lt planescan 30 02 0 1250 125 0 99 0 99 20 50 x2 cathodetoken06 out lt lt planescan 30 02 0 1250 125 0 99 0 99 20 50 p Sol n cathodetoken zou Zinc input files Status READY Open Bun Cancel Fig 4 Post processing the fuel cell run Z ZMesh File View Help Showing variable 2 of 3 KIRI Variable name x2 9 28632 9 28533 EA a Fig 5 Fuel cell result viewed with Zmesh 56 Results 0 12500E 0
18. 3 0 99000E 03 0 00000E 00 to 0 12500E 03 0 99000E 03 0 00000E 00 0 1349 0 13485 f 0 1348 T 0 13475 F 0 1347 F 0 13465 F 0 1346 T 0 13455 F 0 1345 F 0 13445 F 0 1344 F 0 13435 0 0002 0 0004 0 0006 0 0008 0 001 0 0012 0 0014 0 0016 0 0018 0 002 0 scan distance z 0 20000E 04 orsrersrers ooooso0o6 LX 0 00049 Fig 6 Variation of xo2 57 0 12500E 03 0 99000E 03 0 00000E 00 to 0 12500E 03 0 99000E 03 0 00000E 00 0 2863 0 2862 r 0 2861 r 0 286 F 0 2859 f 0 2858 f N 0 2857 r 0 2856 f 0 2855 f 0 2854 F 0 2853 0 0002 0 0004 0 0006 0 0008 0 001 0 0012 0 0014 0 0016 0 0018 0 002 0 scan distance 0 20000E 04 z DALONO OOFO 000 COOOOOO CX 0 06098 Fig 7 Variation of xmo 58 0 12500E 03 0 99000E 03 0 00000E 00 to 0 12500E 03 0 99000E 03 0 00000E 00 112000 110000 F 108000 106000 F a 104000 F 102000 100000 0 0002 0 0004 0 0006 0 0008 0 001 0 0012 0 0014 0 0016 0 0018 0 002 0 scan distance 0 20000E 04 z a 0 06098 Fig 8 Variation of pressure 59 Figs 6 8 show the oxygen x water x2 and pressure p distribution in the sample Also in the directory If the fuel cells runs take to long especially the non DLL version cathode zin try uncommenting the lines Try this for a quicker result N
19. View global format ascii region 1 type box fab 222 Fig 4 Window opened by Select File gt Open MIN showing file min Select File Process to generate the mesh Provided no error occurred the resulting mesh is shown on screen and the mesh file is automatically stored as mtf in the folder The main window now shows the geometry and mesh that has been generated Fig 5 To view the geometry from different angles drag with the mouse File View Help Fig 5 Window opened by Select File gt Open MIN showing file min after rotation of model Dragging rotates the model so that its 3D nature can be observed It can be returned to a standard state simply by pressing the buttons x y or z which shows the view along the selected axis NOTE Can rotate by dragging with mouse or use arrow keys Hold down SHIFT while dragging with mouse to shift rather than rotate the figure Press keys A Z to zoom in out and N M to rotate in the plane To observe the sphere alone press the button Elem 1 in the key at the right which switches off the elements in Region 1 and reveals the embedded sphere in Region 2 see Fig 6 When the elements of a region are switched off a cross appears in the key symbol Click on the key symbol again to toggle the elements back on 1 0000 z 1 0000 xstice ystice zstice Fig 6 Spherical region is exposed by switching off the elements of Regi
20. ZINC Tutorial Manual Introduenonu posts anbe Sep 1 Example 1 DEMENS OSE 2 IntroduetiOni Jet mesen E nan 2 Geometry and TSS IN nan Sa BAN 5 Tipu Pile SZIN Ls EA RENA LL 13 Running the simulation a een s s ae a eaa erten an ede 14 Niewing S 16 A quick look at the resulE 205 enamel yaaa 17 PP tates ean node geist ama E E A EEIN 18 IBT DiS ee 24 More advanced graph plonge sss sese 24 Ifyou don t ike TSTS sirsenis 31 Also the directory ASEAN PA 31 Egampl 2 pi zoelectric examples isisisi inest am 31 Introduction en Aan ta e RA NA E R eee 31 SIM AON E E dere 34 GS EE 37 Companson Strap bit ANA AL SEN EN RN 42 Also in the direct9xaexxxx AE ENERE TA E 43 Example 3 eeeeJfJ 43 IntrodustOn super EA Ma SR INA e AL ALAN NA 43 ZING O lem ena se Sess Bee ISA EL AE NAN ode 48 DEE method hron EEE EE 49 Ex pressions metho EERTE 53 Running thesim l tion oreen 54 Also AM LS directory a aa a 60 Appendix A Description of differential equations sese sees eee eee ee eee eee 60 Introduction This tutorial manual shows a series of real life problems being solved with Zinc It is probably the best place to learn how to use the program though it may be helpful to refer back to the User Manual from time to time Each tutorial corresponds to one of the examples in the examples folder that came with your Zinc distribution These folders contain the input and output files for the simulation as generated by us When you run the simulations
21. Zinc 2 0 Zpp Path DAMGC zincprog new dielectric Input files Zpp file dielectric zpp Sol n file dielectric zou Zinc input files Input file dielectric zin Mesh file dielectric mtf Constants dielectric con Gnuplotfile dielectric gnu Graph files Status READY Run Cancel Fig 19 Creating scans through the data using ZPP The main input files used by zppwin are zpp and zou the simulation solution Zpp will complain if any of these are absent Zpp also reads zin mtf con as used by Zinc so these must also be present As before we can examine files using the view button as shown in Fig 20 Zinc 2 0 Zpp Path DAMGC zincprog new dielectric Input files Zpp file dielectric zpp Sol n file dielectric zou linescan 100 0 101 03 1 0 0101 03 1 0 V cart linescan100 0 101 0 3 1 0 0 101 0 3 1 0 V 2 Sue linescan 100 0101 03 1 0 0101 03 1 0 Vz Inputfile dielectric zin linescan100 0 0 00 10 00 00 1 0 vz Meshfile dielectic mtt linescan100 0 0 0 0 1 0 00 00 1 0 41 eps1 2 eps2 Constants dielectric con a a A A A yplane y xl x2 z1 22 Nx Nz Gnuplotfile dielectric gnu Graph files planescan 200 1 0 10 1 0 10 4 planescan 200 1 0 10 1 0 10 4 0 40 11 eps1 2 e Status READY Open Run Cancel Fig 20 Viewing the dielectric zpp input file which specifies line and plane scans If all is well run the post processor using Run Once the po
22. a po po 0 This will not be needed here however since we are modelling a static situation In the Zinc input file piezo zin it is necessary to enter the C matrix in its component form The eguivalence is shown in Appendix A eguation A3 Thus for example 11 Cun Cu C 11 Cun Cr Cia 11 Cus C13 Cia Etc Which is entered using in the ZIN file using region 1 elements 144 values C 1 0 12600E 12 cll 2 0 00000E 00 l e16 3 0 00000E 00 I G15 etc Note that for example that c16 is just a conventional way to specify the stiffness tensor component C1112 both in N m Thus the Zinc C matrix components corresponding to elasticity are simply equal to stiffness tensor components The relation between tensor and matrix components for elastic stiffness is 11 22 33 23 32 13 31 12 21 matrix 1 2 3 4 5 6 which is the standard used in elasticity textbooks The Zinc input file piezo zin is pretty long for this problem and is not all shown here see examples piezo piezo zin In fact we generated it automatically using a computer program Thus as can be seen above we have set all the stiffness components even those which are zero Also we have input the constants directly rather than using constants in the piezo con file However the piezo con file is used for post processing 33 Simulation Path DAMGC zincprog examples piezo Input files Input piezo
23. agging or use keys A Z or N M normally An example is shown in Fig 11 where the slice intersects both Regions 1 and 2 The elements for Region 2 can be switched off in a slice as shown in Fig 12 10 Z iMesh File View Help y 1 09009 1 0900 Eile View Help y 1 8990 1 0000 z 1 0800 m ustice Zstice Fig 10 The result of changing the colours of the elements of the spherical region 11 File View Help xslice File View Help Fig 12 Example of a slice through Region 1 where elements for Region 2 are switched off In this case the mtf file that is generated by Zmesh exe is 12 IMax 10 JMax 10 KMax 10 I J K Xx y Z RegNo RegUp 0 0 0 1 0000E 00 1 0000E 00 1 0000E 00 4 1 0 0 7 9998E 01 1 0000E 00 1 0000E 00 4 2 0 0 5 9995E 01 1 0000E 00 1 0000E 00 4 3 0 0 3 9991E 01 1 0000E 00 1 0000E 00 4 4 0 0 1 9992E 01 1 0000E 00 1 0000E 00 4 5 0 0 2 2495E 05 1 0000E 00 1 0000E 00 4 6 0 0 1 9996E 01 1 0000E 00 1 0000E 00 4 7 0 0 3 9992E 01 1 0000E 00 1 0000E 00 4 8 0 0 5 9993E 01 1 0000E 00 1 0000E 00 4 9 0 0 7 9996E 01 1 0000E 00 1 0000E 00 4 1 10 0 O 1 0000E 00 1 0000E 00 1 0000E 00 4 0 0 1 0 1 0000E 00 7 9998E 01 1 0000E 00 4 BUG Sunset The MTF file is not mysterious it is just a text file showing the locations of nodes their Region and their connectivity to neighbouring regions This f
24. an GI CR alu ui CH Ch Gla CH GU lau Ga BIE Ce GJL a It then follows that V CVU CUR a 12 N A8 The term V CVU is regarded as an 1 x N column vector When N 2 it would be written cpus Vi CVU Se e A9 Gus J The matrix C used by the ZINC code has components of the form Cn which are given by aipj 61 CT Ge 1 je A10 oiBj Boundary conditions In order that the differential equations A2 have a unique solution a set of boundary conditions must be imposed on the external boundary of the system Such boundary conditions are assumed to be of the form generalised Neumann condition n CVU qU g A11 where n is the outward unit vector on the external boundary and where g is a prescribed scalar function and g is a prescribed vector function dimension N In component form this condition is written HC Ur ES MAKE AN A12 62
25. ants dielectric con NL file lt No NL file gt Restartfile None initial condition supplied Output files Outputfile dielectric zou List file dielectric zls Status READY Bun Cancel Fig 14 Zinc window with dielectric zin selected The program checks all input files are present and shows the filename of each To run Zinc at least an input file zin and a mesh file mtf are needed The constants and non linearity NL files are optional A restart file is needed only in the case key rdu 1 which means continue where you left off Since we ve set this to zero s above this is not needed You can view any file but clicking on the view button next to it In Fig 15 we are viewing the zin file Note that the program does not allow you to edit files you should use a text editor to do that You can just use Notepad but a programming text editor is probably better We personally use Emacs but there are many free editors to chose from Zinc 2 0 Path DAMGC zincprog new dielectric l nstep 10000 Input files ooo lomega 0 1 Input file dielectric zin nreg 4 Meshfile dielectic mtf key rdu 0 Constants dielectric con scale 1 EG S 5 key Q 1 NL file lt No NL file gt view labels lt V Restartfile None initial condition supplied 1 elements 3 values C 8 854187818e 12 Output files 8 854187818e 12 Outputfile dielectric zou 8 854187818e 12 L
26. considering 6 13 For example if the B11 token is called for the function implements equation 9 This function calls Nlinit see cathodetoken f90 which simply reads in the cathodetoeken con file in order to obtain the material properties stored there This illustrates how user written functions can call other functions and perform input output operations as required The cfun function in C is written as double cfun char label double x double y double z double ur 3 double dur 3 3 int nvar int istep int length double D11 D12 D13 D21 D22 D23 D31 D32 D33 double x1 x2 p x3 wl w2 w3 m denom NLinit 1 IG r x1 ur x2 ur p ur 2 0 1 l x3 1 x1 x2 m ml x1 m2 x2 m3 x3 wl ml x1 m w2 m2 x2 m w3 m3 x3 m denom x1 D12tilde D13tilde x2 D12tilde D23tilde x3 D13tilde D23tilde D11 w2 w3 w2 w3 x1 D23tilde w2 w2 x2 D13tilde w3 w3 x3 D12tilde denom D22 wl w3 wl w3 x2 D13tilde wl wl x1 D23tilde w3 w3 x3 D12tilde denom D33 wl w2 wl w2 x3 D12tilde wl wl x1 D23tilde w2 w2 x2 D13tilde denom D12 wl1 w2 w3 x1 D23tilde w2 w14tw3 x2 D13tilde w3 w3 x3 D12tilde denom D13 wl1 w2 w3 x1 D23tilde w3 wl w2 x3 D12tilde w2 w2 x2 D13tilde denom D23 w2 wl w3 x2 D13tilde w3 wl w2 x3 D12tilde wl wl x1 D23tilde denom D21 D12 D31 D13 D32 D23 if eq
27. dielectric simulation we have generated 33 files Dielectric zin Dielectric min Dielectric mtf Dielectric con Dielectric mls Dielectric zls Dielectric zou Dielectric zpp Dielectric gnu Dielectricxx out xx 01 to 10 Dielectricxx emf xx 01 to 10 This might seem like rather a lot of files but this structure provides traceability Imagine 6 months down the line you look at file dielectric01 emf and wonder how on earth did I generate this graph By looking at dielectric zpp we can see that it is a line scan between such and such points of the electric potential V But how did we get this simulation in the first place By looking at dielectric zin physics and dielectric min geometry we can see exactly how the simulation was run We can then reproduce the same simulation or a variation of it easily Thus the Zinc text based file based paradigm give clear traceability and is superior to the nested dialog box approach used by most FE packages With these other packages it is common to end up with graph picture files and not know how you generated them More advanced graph plotting We ve seen that Zpp can generate basic graphs which give a good immediate insight into what is happening in the simulation For publication however these graphs are probably not sufficient For a start it will be necessary to re label the y axis of most graphs Labels like 1 epsl 2 eps2 Vz don t mean much to most people Zinc
28. dielectric gt gnuplot comp gnu D MGC z inc prog examples die lectric gt Fig 26 Running gnuplot under batch control If all goes well the graph files are silently created in the same directory The three graphs generated are shown in Figs 27 29 1 T T 0 9 F J 0 8 F J 0 7 F 4 06 F J V Volts o a T 1 03 F J 0 2 F 4 01 F 4 0 0 5 1 1 5 2 Scan distance m 27 Fig 27 Graph plotting with better axis labelling 0 35 x 0 1m x 0 0m eae Sisk A 1 aa 0 45 0 5 Ez V m gt Less 0 5 1 Scan distai nce m Fig 28 Graph with two curves and a legend 1 0 5 0 0 5 1 1 0 5 Dz C m 2 0 x m Fig 29 Colour coded graph of Dz 0 5 1 3e 012 3 5e 012 4e 012 4 5e 012 5e 012 5 5e 012 6e 012 6 5e 012 7e 012 Figs 27 and 28 were prepared in black and white with dashed lines which is more usual for journal papers This is just a taste of the graphs one can prepare using gnuplot The full range of graphs commonly found in papers inset graphs multiple axes histograms vector fields error bars etc can all be generated using gnuplot 28 For completeness Fig 30 shows the same 3 graphs rendered as eps encapsulated postscript and being viewed using ghostview freeware EPS is the normal graph format when using LaTeX it arguably gives the best quality graphs
29. e Ox x dm I mm nRT ox dx mox dy dy m dy dz dz maz Pm P X m 2 Ox x m Ba am 2 MRT dx ox max dy dy may dz daz mor f 0 From 5 the Neumann boundary conditions that will exist on non fixed surfaces of the simulation are n 1 n J with the flux j given in 1 The is the so called convective boundary condition we want to solve in this problem The way we have set things up the Zinc C matrix encodes the div jx term in 1 while the Zinc f vector encodes the pu V Since C appears in the boundary condition 5 n j 0 is the natural boundary condition we get There are two ways to solve non linear problems in Zinc 1 Write a set of subroutines to encode the non linearity compile to a dynamic link library DLL and run with Zinc 2 Write expression in the Zinc input file encoding the non linear behaviour In this tutorial we will show both methods and furthermore we will prepare the DLL using both Fortran and C implementations We will first consider the DLL method DLL method We will implement this method using the name cathodetoken The input file cathodetoken zin looks like this key_sim 1 nvar 3 49 nreg 3 key rdu 0 scale le 3 omega 0 1 nstep 30000 nstride 100 nstrideNL 1000 labels xl x2 p region 1 elements 21 values C Row 1 of B L det SBI 1 2 1 2 B11 31 3 B11 L 2 IS SBIZ 1 2 2 2 B12 3 2 3 SB12 Row 2 of B 21 1 1 SB21 2 2 1 2
30. ey plot 2 ftol le 3 N x1 y z1 x2 y2 z2 linescan 100 1 0 0 1 0 4 10 Od 0 4 ux linescan 100 LO Dil 0 4 TO 0 1 0 4 uy linescan 100 1 0 Og 0 4 1 0 0 1 0 4 uz linescan 100 WhO Ol 0 4 T0 Oial 0 4 V linescan 100 aD Del 0 4 TRO Pal 0 4 Vz linescan 100 1 0 Ol 0 4 120 0 1 0 4 Vx linescan 100 05 O20 1 0 0 5 0 0 1 0 V linescan 100 DaS 00 1 0 0 45 0 0 1 0 ux linescan 100 aS Del LO 0 5 0 0 1 0 uz linescan 100 1 0 0 1 0 4 1 30 0 1 0 4 amp l cll1 2 cllinc uxx 1 c12 2 cl12inc uyy 1 c13 2 1 e31 2 e31inc Vz cl3inc uzz amp linescan 100 LO Dul 0 4 T0 0 1 0 4 amp f1 c12 2 c12inc uxx 1 c22 2 c22inc uyy 1 c23 2 c23inc uzz amp 1 e32 2 e32inc Vz The first batch of scans are along x and plot the raw output and also the electric fields Recall that Vz means dV dz so Vz is the electric field z component Ez The next batch of line scans are along z between the two electrode plates 35 The final two scans are stress plots specifically G1 0xx and 62 oyy Let us consider the constitutive law of piezoelectrics O Cy 2 i E Eal 05 Cor Cn 3 E x E 1 O3 31 C32 C33 E3 33 E za 2 O C44 E Eny E O C55 E 15 OG Can NE Multiplying out gives O CE tp 0363 763 E3 O CE tc 03383 CE which are precisely the formulas used in the last two linescans The expression f1 c12 2 c12inc means take c12 when we are
31. ficients it is possible to enter these as complex expressions This approach is used in cathode zin which does not require linking to a DLL file Thus the full expression for C matrix element C1111 is region I elements 21 values C x1 D11 D13 RT 1 p ml x1 BT x2 m23 m3 m3 x1 2 x1 D23 m2 2 x2 D13 m3 2 1 x1 x2 D12 amp m1 x2 m23 m3 m3 x1 D23 m3 m1 x1l m2 x2 D12 m2 2xx2 D13 amp ml3 x1 m23 x2 m3 2 x1l facl x2 fac2 fac3 p ml 111 53 This can be derived from 6 and 9 and using the D coefficient expressions just after Figure 1 Everything needs to be substituted in so that the C and f coefficients are written purely in terms of 1 The variables to be solved for x1 x2 p in this case a labels command sets this up 2 Constants in the cathode con file We have introduced a few new extra constants including facl fac2 and fac3 here see cathode con Writing out the expressions explicitly like this is ok but it makes them more complicated to read In the program based approach the expressions were built up using intermediate variables which is not possible here Also the DLL file method will work faster since the expressions are compiled into machine code On the other hand expressions like the one shown above need to be repeatedly interpreted by Zinc when it runs In summary the programming DLL method should be used for more complicated non linearity systems and fuel cells certainly
32. hat needs to be implemented 43 The fundamental equation for diffusion of n species is AO V j pu Va R k 1 n 1 1 where p denotes the density kg m OM P P is the mass fraction of species k u is the mass based mean velocity vector m s and R is the reaction rate kg m s Also Ji P U u describes the diffusion driven transport where ux is the mean velocity of the k species The flux jx is given by 5 VT z hi D PO 2 Dyd k 1 n 2 j l where the driving force vector d is given by Vp a vx l so j 1 n 3 where x c c is the molar fraction of species j and p is the pressure Pa and where cj is the concentration moles m of the T species and c is the mean molar concentration defined by n eye k 1 Thus the driving force for diffusion comes from two parts the variation in species concentration x and the variation of pressure p The molar and mass fractions satisfy the relations yx 1 Yo 1 k 1 k l and can be expressed in terms of each other using 44 The total density p in 1 can be written as p n m X P RTS 3 We now make several assumptions to simplify the problem Firstly we assume steady state s0 that d dt 0 Next we assume no temperature variation for that VT 0 Finally we assume reaction rates Rx of zero With these assumptions substituting 3 into 2 and then 2 into 1 gives v dL x 0 H 0 k 1 n 1 4 J In our
33. hrough the sample You can look at the other graphs using the Next or Prev buttons on the graph window or by returning to the listbox in the zppwin window Note that there is no need to save any of these files they are automatically saved when we run zppwin In fact two types of files are saved in this case the dielecticxx out files which contain a table of numbers corresponding to the scan and the dielectricxx emf which contain the actual graphs Graph files are only generated if key_plot is 1 or 2 key_plot 1 generates encapsulated postscript files eps while key_plot 2 as used here generates enhanced metafiles emf Graph previews within zppwin are only available for emf files key_plot 2 To view eps files you ll need to install ghostview a free program It depends which word processor you intend to import graphs to Eps graphs are suitable for LaTeX while emf are suitable for Microsoft products like Word and Powerpoint Both these graph formats are vector based Fig 22 is the second of the two plane scans showing D the z component of the D field 22 Zinc 2 0 Zpp Graph y 0 00000E 00 Nag RER LT PRET 277 LAL CZ E LLIT EFAS N a a UU U A a GD a UU 1 0 00000 000 225333533533 Id MI IND KAKA NI NING Fig 23 Plot of Dz C m through the central plane of the system Let us look at the command which generated this plot planescan 2 0 0 1 0 1 0 51 0 1 0 40 40
34. ile is read in by the program Zincwin see later Normally the user does not have to look at this file Input file zin ZINC requires information that controls the performance of the software materials data for each region and data relating to boundary conditions that are not set by default For the dielectric problem this information is given in the following zin file key sim 1 I for static or steady cyclic states 2 for transient problems nvar 1 No of unknown variables in this case just one electrostatic potential nstep 10000 No of iterations omega 0 1 Successive over relaxation parameter lt 1 Default 0 1 nreg 4 No of Regions in model key_rdu 0 0 start from initial condition in zin 1 continue from previous sim scale 1 scaling factor for lengths found in mtf labels V List of unknown variables nvar in number region 1 elements 3 values C For elements of Region the are 3 non zero Icoefficients of ZINC matrix C L1 11 epsl 21 2 epsl 38 BE epsl region 2 elements 3 values C For elements of Region 2 the are 3 non zero Icoefficients of ZINC matrix C LITT eps2 2 1 2 eps2 L 3 1 3 eps2 region 3 nodes 1 values For nodes of Region 3 there is I value 1 0 for unknown variable L 1 0 13 region 4 nodes 1 values For nodes of Region 4 there is 1 value 0 0 for unknown variable 1 0 0 Init Initial value 0 0 for unknown variable 1 0 0 In this problem as specified i
35. ist file dielectric zls region 2 elements 3 values C 111111 1 770837564e 11 Status READY Open Run Cancel Fig 15 Viewing one of the input files If all seems to be in order it is now time to run Zinc Press the run button In the console window you will see output shown in Fig 16 15 cx zincwin 9826 9836 9846 9856 9866 2870 9886 98908 92900 9916 9926 2930 99408 9956 2960 2970 9986 9996 16666 3812 3E 14 371486E 14 362953E 14 353629E 14 344153E 14 333995E 14 3229 1E 14 312764E 14 4 36365GE 14 294653E 14 286164E 14 278516E 14 27 851E 14 262586E 14 253513E 14 244433E 14 236299E 14 2292 8E 14 222897E 14 Final Energy Final energy from Q Writing final state to exe 6 544522E 83 6 544522E 83 6 544522E 83 6 544522E 83 6 544522E 3 6 544522E 3 6 544522E 3 6 544522E 83 6 544522E 83 6 544522E 83 6 544522E 83 6 544522E 83 6 544522E 83 6 544522E 83 6 544522E 83 6 544522E 63 6 544522E 63 6 544522E 3 6 544522E 83 ZOU Zinc completed successfully OK Fig 16 The numbe should ge The 3 an side and numbers s up to 10 000 time steps Output console window showing simulation r of time steps is show The sec t low when an approxima right hand side of t hould be equal to each solution output using two calculations other If Viewing has been found to good T
36. istributed charge namely V D 0 3 where D is the electric displacement whose value is given by the constitutive relation D ce E 4 where E is the electric field is the permittivity of free space and is the relative permittivity a dimensionless material property The electric field is related to the electric potential V according to the relation E VV 5 The relations 3 5 imply that the electric potential must satisfy the equation V ee VV 0 6 In component form this equation is written 2eov e 7 The relation 7 results from 2 if the following identifications are made See Zinc User Manual Chapter 4 Nel a SO PE VEV EE 8 where 5 is the Kronecker delta symbol having the value 1 if i j 0 otherwise It follows from A10 that the C matrix used by ZINC has the values given by NI Ciaj 5Cy E 9 It should be noted that in ZINC all coefficients appearing in the general differential equations 2 are zero unless non zero values are set in the ZINC input files zin where denotes an input filename introduced by the user The differential equation 7 has a unique solution only if suitable boundary conditions are applied The boundary condition assumed by Zinc is w du i Gr U j where n is the unit normal on the surface It is easy to check that with q U ee q g 0 this becomes n D 0 or ov U on which is the usual open boundary cond
37. ition used in electrical problems It implies the electric field E will be parallel to the outer edges of the simulation In the example to be solved this boundary condition applies everywhere on the external surface of the system except on two opposing planes of the cuboid where differing uniform electric potentials are to be imposed These Dirichlet boundary conditions are imposed simply be fixed the values of the nodes at the appropriate potential which will be simply I V and OV In ZINC any unset components of the C a f q g matrices defaults to zero so the condition n D 0 occurs automatically At the surface of the spherical region the normal component of the electric displacement must be continuous across the boundary and the electric potential must also be continuous across this boundary Provided we leave q and g at their default zero values this condition will also occur automatically The boundary conditions that apply are illustrated in Fig 1 Fig 1 Schematic diagram of geometry and the boundary conditions The default condition on the external boundary of Region I is therefore 0 implyingthat n D 0 10 n where n is the unit outward normal to the free surface On the plane z I V 1 11 On the plane z 1 V 0 12 On the interface between Regions 1 and 2 the following conditions must be satisfied V l l V V amp AN implying that n D is continuous across the interface 13 an
38. lot from its desktop start menu shortcut producing the window in Fig 25 38 gnuplot ex File Plot Expressions Functions General Axes Chart Styles 3D Help Reptot Jl Open Jl Save J ChDir Jl Print Prise Jl Prev Jl next JE GNUPLOT Version 4 4 patchlevel 2 last modified Wed Sep 22 12 10 34 PDT 2010 System MS Windows 32 bit Copyright C 1986 1993 1998 2004 2007 2010 Thomas Williams Colin Kelley and many others gnuplot home http www gnuplot info faq bugs etc type help seeking assistance immediate help type help plot window hit h ID MGC zincprog examples dielectric Imuplot load comp qnu muplot gt _ Fig 25 The gnuplot window 26 First make sure you are in the right directory using the pwd command as shown and change directory using the cd command if necessary Then use the load command to run the script file you ve prepared comp gnu in this case If all goes well the graphs are created in the same directory Tip Note that we did not type in the full file path for the out files in comp gnu This is because gnuplot works in the current directory by default If you are familiar with using the command prompt you can simply run gnuplot at the command line with the script file as an argument fig 26 cx Command Prompt lo x Microsoft Windows XP Version 5 1 2666 CO Copyright 1985 2801 Microsoft Corp D N gt cd mgc zincprog examples dielectric D MGC z inc prog examples
39. n where n is now the unit outward normal to the interface 0 0n denotes the normal derivative and where is the relative permittivity of Region i The conditions 13 ensure that the normal component of the electric displacement vector and the tangential component of the electric potential are continuous across the surface These continuity conditions are true by default for all such internal interfaces and normally the user only needs to think about applying boundary conditions at the outer surfaces The other boundary condition is the Dirichlet boundary in which the unknown variables U are set at fixed values In this tutorial we will use the Dirichlet boundary at the top and bottom faces of Fig 1 and the Neumann condition n D 0 at all other faces Geometry and meshing Having installed ZINC onto your computer it is useful to generate a new folder having name ZINC Runs or similar into which all simulations will be generated and saved To deal with the first example of this tutorial a new folder within ZINC Runs is created having the filename Dielectric This file will now be denoted in the following text by Input file min To represent the geometry shown in Fig 1 by finite elements the following input file must be written and stored in ZINC_Runs Dielectric The first step defines the global mesh of cubic elements which represents the entire region of the problem The remainder of the file defines pa
40. n Eq 8 the C matrix coefficients are just the permittivities In this simple case all regions are isotropic so permittivities 11 22 33 Note that coefficient 1111 is permittivity 11 1212 is permittivity 22 etc The first and third index numbers refer to the variable while the second and forth relate to spatial dimension See User Manual Chapter 4 In this case there is only one variable so all indexes are Ixly It would of course be easy to simulate an anisotropic dielectric rather than an isotropic one as in this example The constants eps and eps2 are in the constants file dielectric con They specify the permittivity of both regions You can easily change these permittivities without altering the dielectric zin file Running the simulation Having prepared the Zinc input file it is time to run the actual simulation using Zinc shown in Fig 13 Zinc also pops up an additional console window As Zinc runs information about convergence is written into the console window Path none Input files Input file none Mesh file none Constants none NL file none Restartfile none Output files Outputfile none List file none Status READY Cancel Fig 13 Zinc core window Now use the Open button to select the file dielectric zin Fig 14 14 Zinc 2 0 Path DAMGC zincprog new dielectric Input files Input file dielectric zin Mesh file dielectric mtf Const
41. nodes red are clamped and fixed at 0 V We consider a very simple geometry as shown in Fig The outer area is a piezoelectric known as PZT 5H while the inner area is a fictitious material with 10 times greater stiffness and the same piezoelectric constant and permittivity The top and bottom row of nodes shown in Fig 1 are fixed in terms of all the variables the displacements are fixed at zero ux Uy uz 0 while the potential is set at 1 V top and 0 V bottom In Zinc it is of course possible to fix only a selection of the variables on each node Here we will fix all variables We will set the Zinc matrices g q 0 the default which results in a zero stress zero charge boundary condition at the x and y minimum and maximum edges The simulation is rather similar to that of Example 1 but with extra elastic interactions In Chapter 4 of the User Manual the full C matrix for a multiferroic material is derived resulting in equation 4 36 We won t need all that matrix since we are only solving for a piezoelectric In fact the C matrix is given by the top left portion of that multiferroic matrix C11 C16 GR C16 C12 C14 C15 C14 C13 E11 21 31 C61 66 65 C66 C62 64 C65 C64 63 16 26 36 C51 C6 C55 C56 C52 C54 C55 GRA C53 15 25 35 C61 66 65 66 C62 C64 C65 C64 C63 Cip 26 36 C21 C26 C25 C26 622 C24 C25 C24 23 12 22 32 CA C46 CaS C46 C42 C44 CAR C44 Ca3 14 24 34 C51 C56 C55 CKG C52 C54
42. of course just labels the axes with the formula you give it it doesn t know the physics nor the right units to use Beyond this it might be necessary to compare to curves on the same graph with a legend to distinguish them Or we might want to compare FE simulation results with theory or experimental data obtained from elsewhere Zinc creates its graphs by running gnuplot an open source graph plotting package bundled with the Zinc distribution The script file used to generate the graphs is shown called gnu and can be viewed by pressing the view button shown in Fig 24 24 Zinc 2 0 Zpp setterm emf Path DAMGC zincprog examples diele set nokey Input files set hidden Zpp dielectric zpp set view 30 30 Sol n dielectric zou set output dAmgc zincprog examples dielectric dielectric01 er settitle 0 10100E 00 0 30000E 00 0 10000E 01 to 0 1010 SUPE LET setxlabel scan distance Input dielectric zin i set ylabel V Mesh dielectric mtf i plot d mgc zincprog examples dielectric dielectricO1 out u 1 5 Constants dielectric con set output dAmgckzincprogfexamplesidielectricidielectricO2 er J Gnuplot dielectric gnu lt gt Graph files dielectric01 out lt lt linescan 100 0 101 0 3 1 0 0101 03 1 0 v dielectric02 out lt lt linescan 100 0 101 0 3 1 0 0 101 0 3 1 0 v 2 dielectric03 out lt lt linescan 100 0 101 0 3 1 0 0 101 0 3 1 0 vz dielectric04 out lt lt linescan100 00 00 10 00 00 1 0 vz
43. on 1 The elements of Region 2 can be switched off by pressing the button Elem 2 and so on Nodes can also be toggled on or off by clicking on the key symbols marked node 1 etc The geometry of the nodes is best seen by dragging Repeated pressing of all buttons simply toggles between two possible states visible or invisible In Fig 7 is shown the result of retaining the elements of Region 2 but showing the nodes Ze 1 9000 y 1 00989 Ea m me Fig 7 Elements and nodes of spherical region shown by a suitable choice of active buttons The colours of the nodes and elements can be changed by right clicking on node or element buttons on top right hand side of the window The colour is selected by specifying the mix of primary colours using the Set Colour dialogue box see Fig 8 The result of changing node colours is shown in Fig 9 while the result of changing element colours is shown in Fig 10 Basic colors BT ERT BAG GE TTT HANNENE HEHHEHE EH Custom colors efo Red 255 ser ae Green 0 F Pustom C ColorlSolid Lum 120 Blue fo Cancel Add to Custom Colors Fig 8 Dialogue box for selecting the colours of elements and nodes Pressing the xslice yslice or zslice buttons shows only a slice of the model in the x y or z directions The location of the slice can be changed using the Up Down arrows on the keyboard While doing this you can continue to rotate the system by dr
44. owing graphs which show good agreement ux pm 42 2 5 T T T Zinc Comsol uz pm F 0 5 F z cm Also in the directory Piezoslap zin Same simulation but using the Incomplete LU factorisation GMRES method Example 3 Fuel cell Directory examples fuelcell Cathode zin self contained simulation without programming Cathodetoken zin simulation with DLL written in Fortran Cathodetokenc zin simulation with DLL written in C Introduction We now consider a multiphysics non linear example the fuel cell Fuel cells are ideal cases for general purpose multiphysics solvers like Zinc since the precise set of equations and boundary conditions is under constant debate It is often necessary to solve a set of equations only recently derived which cannot be done with traditional FE packages that have built in modules corresponding to common problems Zinc allows you to write down a new set of equations and immediately solve them without having to wait for developers to explicitly add the new equation sets in This is precisely the approach we took in a recent fuel cell project at NPL We derived the equations ourselves from first principles and them implemented them in Zinc The full derivation of the fuel cell equations is given in NPL report MAT RES 108 Here we will briefly summarise the derivation and show the working equation set t
45. rts of the system that are needed to define homogeneous regions having specific properties and boundary and interface conditions In the following code comments are included to provide some explanation More details of the min file are given in the Zinc User Manual global format ascii xmesh 1 1 0 2 The x range 1 lt x lt lis divided into elements of length 0 2 ymesh 1 1 0 2 The y range 1 lt x lt 1 is divided into elements of length 0 2 zmesh 1 1 0 2 The z range 1 lt x lt 1 is divided into elements of length 0 2 part 1 Defines a cubic box having sides of length 2 as Region 1 region 1 type box fab 2 2 2 end part 2 Defines a sphere having radius 0 5 as Region 2 inserted into Region 1 region 2 type sphere fab 0 5 surface region 1 end part 3 Defines Region 3 to be the upper plane z 1 region 3 type boundzup end part 4 Defines Region 4 to be the lower plane z 1 region 4 type boundzdn end endfile The next step is to generate a mesh using zmesh exe which appears as two windows shown in Figs 2 and 3 The ZMesh window Fig 2 enables the user to generate the required mesh while the Output window Fig 3 provides information on locations of files being used Fig 2 Window first opened by Zmesh Welcome to 2Mesh 3 6 2011 Fig 3 Output window first opened by Zmesh Select File gt Open MIN and open the file min whose contents are then displayed see Fig 4 File
46. s are still present on disk we can replot these files in different ways as needed The best Way is to create a new gnuplot script file and run it through gnuplot In principle one could just alter dielectric gnu in this example but this file will be overwritten next time zppwin is run so its probably better to use a new file We will now prepare a file to 1 plot the first linescan with better labelled axes 2 plot the two Ez linescans in the same graph with a legend 3 plot the Dz planescan as a colour grid rather than a surface plot Create a text file called comp gnu and enter the following 25 set term emf mono dashed set output vscan emf set xlabel Scan distance m set ylabel V Volts unset key plot dielectricOl out u 1 5 w 1 set output Ezscan emf set ylabel Ez V m set key plot dielectric03 out u 1 5 dielectric04 out u 1 5 t x set output Dzmap emf set title Dz C m 2 set xlabel x m set ylabel z m set size square set pm3d map splot dielectricO7 out u 1 3 4 By now you should have a rough idea of what this is doing but consult the gnuplot manual for further details We are simply plotting the files dielectricO1 out dielectric03 out dielectric04 out and dielectric07 out with slightly more advanced formatting There are several ways to run gnuplot If you are unfamiliar with the command prompt the simplest is to start gnup
47. s in the zin file so all other matrices should be full of zeros in this simple system 16 Zinc 2 0 Path DAMGC zincprog new dielectric Input files Input file dielectric zin Mesh file dielectric mtf Constants dielectric con NL file lt No NL file gt Restartfile None initial condition supplied Output files Outputfile dielectric zou List file dielectric zls Prepared by Zinc 2 0 imax jmax kmax 10 10 10 Note that i 0 imax etc xrange 1 0000000000000000 1 0000000000000000 yrange 1 0000000000000000 1 0000000000000000 zrange 1 0000000000000000 1 0000000000000000 Control Value User set nvar lo nreg 4 T key_rdu oT key_sim iW nstep 10000 T tstep 0 00000E 00 F key db DE omega 0 10000E 00 T scale 0 10000E 01 T Status READY Open Run Cancel Fig 17 Viewing the output file in a completed run As Zinc shows us the output file containing the simulation result in zou This file is not viewable because it is not intended to be human readable There s no mystery involved however and you can look at this file using a text editor Note that the output files are automatically saved to disk with the names zou zls There is no need to save the files Zinc is now complete and you dismiss the window by pressing Cancel A guick look at the result We can get a simple view of the data using Zmesh In Zmesh open the zinc output file dielectric zou using the File
48. scan 100 0 101 0 3 1 0 0 101 0 3 0 30000E 00 0 10100E 00 0 30000E 00 0 70000E 00 0 1656 dielectric03 out lt lt linescan 100 0 101 0 3 1 0 0 101 0 3 0 32000E 00 0 10100E 00 0 30000E 00 0 68000E 00 0 1771 dielectric04 out lt lt linescan 100 00 0 0 1 0 00 00 11 0 34000E 00 0 10100E 00 0 30000E 00 0 66000E 00 0 188 dielectic05 out lt lt linescan 100 00 00 10 00 00 1 p360006 00 0 10100E 00 0 30000E 00 0 64000E 00 0 1996 dielectric06 out lt lt planescan 2 0 0 1 0 10 1 0 1 0 0 38000E 00 0 10100E 00 0 30000E 00 0 62000E 00 0 2106 dielectric07 out lt lt planescan 200 1 0 10 1 0 10 0 40000E 00 0 10100E 00 0 30000E 00 0 60000E 00 0 2221 0 42000E 00 0 10100E 00 0 30000E 00 0 56000E 00 0 2337 Solr fle diglenticzou Zinc input files Status READY Open Bun Cancel Fig 21 Once a run is complete the Graph files list is active and we see the resulting graphs But it s much more fun to view this as a graph To do this click Graph at the bottom right of the window Another window pops up Fig 22 showing the graph In this case it is a graph of voltage versus distance along the z axis The voltage goes from 0 to 1 V as we expect and there is a kink in the middle where the high permittivity sphere is located 21 Zinc 2 0 Zpp Graph L 0 10100E 00 0 30000E 00 0 10000E 01 to 0 10100E 00 0 30000E 00 0 10000E 01 1 scan distance Fig 22 Variation of voltage t
49. st processing is complete the graph files list becomes active This list shows us which datafiles correspond to which of 20 the selected scans For example dielectricO1 out will contain the first linescan requested Clicking on this line as illustrated in Fig 21 allows us to view the file Zinc 2 0 Zpp HFINAL ISTEP 10000 Path DAMGC zincprog new dielectric 0 00000E 00 0 10100E 00 0 30000E 00 0 10000E 01 0 8654 nput files 0 20000E 01 0 10100E 00 0 30000E 00 0 98000E 00 0 1097 REE BREN 0 40000E 01 0 10100E 00 0 30000E 00 0 96000E 00 0 2194 0 60000E 01 0 10100E 00 0 30000E 00 0 94000E 00 0 3291 0 80000E 01 0 10100E 00 0 30000E 00 0 92000E 00 0 4388 0 10000E 00 0 10100E 00 0 30000E 00 0 90000E 00 0 548 i 0 12000E 00 0 10100E 00 0 30000E 00 0 88000E 00 0 658 Inputfile dielectric zin 0 14000E 00 0 10100E 00 0 30000E 00 0 86000E 00 0 767 Mesh file dielectric mtf 0 16000E 00 0 10100E 00 0 30000E 00 0 84000E 00 0 8776 Constante didlecmeuah 0 18000E 00 0 10100E 00 0 30000E 00 0 82000E 00 0 987 0 20000E 00 0 10100E 00 0 30000E 00 0 80000E 00 0 1096 Gnuplotfile dielectric gnu 0 22000E 00 0 10100E 00 0 30000E 00 0 78000E 00 0 1208 S 0 24000E 00 0 10100E 00 0 30000E 00 0 76000E 00 0 1321 i saias 0 26000E 00 0 10100E 00 0 30000E 00 0 74000E 00 0 143 dielectric0 out lt lt linescan 100 0101 03 10 0101 0 3 0 28000E 00 0 10100E 00 0 30000E 00 0 72000E 00 0 1546 dielectric02 out lt lt line
50. step 10000 Key rdu 1 in the ZIN files of this directory This will continue from a previous simulation rather than starting from scratch Appendix A Description of differential equations The differential equations solved by ZINC have the form V CVU L aU f Al which are written explicitly in component form as follows GPO j ta UP it a 12 N A2 where summation over 1 2 N is implied by repeated Greek superscripts and where summation over 1 2 3 is implied by Arabic suffices Here f df 0x The structure of ij the 3N x 3N C matrix when N 2 is Ci Cp G G Cp G Cy Cy Ca i amp Cy Gi Cy Ca Ca Co Ci Ch C3 Ci Ch Ch G Ch C3 Cr Co C3 21 21 21 22 22 22 L C3 Cy C3 C31 C3 C3 A3 Note that 0 Cy Caip in the notation of the User Manual Chapter 4 C components in the input ZIN file need to be entered in the order i B j 60 The general form of the matrix C is an NxN array of 3 x 3 matrices The term VU expands in component form as vu us B 1 Di ea N j 1 2 3 A4 The term VU is regarded as a 1 x 3N column vector When N 2 it would be written EEE ailn sr ag VU A5 c HN c ws 72 KE The term CVU then expands in component form as follows CVU CPUE NELS A6 The term CVU is regarded as a 1 x 3N column vector When N 2 it would be written IB 78 CU Ch Ch oa CR sil Su S amp BA Ce GBP ag
51. ual label length SB11 return p ml x1 D11 D13 RT else if equal label length S B12 52 return p ml x1 D12 D13 RT else if equal label length SB21 return p m2 x2 D21 D23 RT else if equal label length 5 B22 return p m2 x2 D22 D23 RT else if equal label length SQ1 return ml x1 RT x1 D11 1 ml m D13 1 m3 m x2 D12 1 m2 m D13 1 m3 m D13 1 m3 m else if equal label length 5 Q2 return m2 x2 RT x1 D21 1 ml m D23 1 m3 m x2 D22 1 m2 m D23 1 m3 m D23 1 m3 m else if equal label length S E return p m kappa RT eta else printf Error in Cfun token not recognised n exit 1 Note that the token is passed in as a fixed length variable length 1ength rather than being null terminated For this reason we have written a function equal which the user should examine in the cathodetokenc c file We used the free gcc compiler including gfortran to compile the Fortran and C DLLs gfortran c fnounderscoring cathodetoken f90 gfortran s shared mrtd o cathodetoken dll cathodetoken o and gcc c cathodetokenc c gcc shared mrtd o cathodetokenc dll cathodetokenc o The batch files compile bat and ccompile bat included in the fuel cell directory contain these commands Expressions method Instead of entering literal constants for the C f coef
52. wl ml x1 m w2 m2 x2 m w3 m3 x3 m denom x1 D12tilde D13tilde x2 D12tilde D23tilde x3 D13tilde D23tilde tw3xx2 x3 D12tilde denom tw3xx2 x3 D12tilde denom tw2xx2 x2 D13tilde denom x2 D13tilde w3 2 x3 D12tilde denom x3 D12tilde w2 2 x2 D13tilde denom x3 D12tilde wl 2 x1 D23tilde denom D11 w2 w3 2 x1 D23tilde w2 2 x2 D13tilde 4 D22 wl w3 2 x2 D13tilde w1 2 x1 D23tilde 4 D33 wl w2 2 x3 D12tilde w1 2 x1 D23tilde 4 D12 wl w2 w3 x1 D23tilde w2 w1 w3 D13 wl w2 w3 x1 D23tilde w3 wl w2 D23 w2 wl w3 x2 D13tilde w3 wl w2 D21 D12 D31 D13 D32 D23 if label eq B11 then Cfun p m1 x1 D11 D13 RT else if label eg B12 then Cfun p m1 x1 D12 D13 RT else if label eg B21 then Cfun p m2 x2 D21 D23 RT else if label eg B22 then Cfun p m2 x2 D22 D23 RT else if label eq SQ1 then Cfun m1 x1 RT x1 D11 1 m1 m D13 1 m3 m amp 51 x2 D12 1 m2 m D13 1 m3 m amp D13 1 m3 m else if label eq Q2 then Cfun m2 x2 RT x1 D21 1 m1 m D23 1 m3 m amp x2 D22 1 m2 m D23 1 m3 m K D23 1 m3 m else if label eq SE then Cfun p m kappa RT eta else print Error in Cfun token not recognised trim label stop endif end function Cfun The operation of this function should be clear by
53. y s constant for O2 dissolved in H2O Pa m3 mol mO2 32 0x107 Molar mass of 02 kg mol mH2 2 0x10 Molar mass of H2 kg mol mH20 18 0x10 Molar mass of H2O kg mol mN2 28 0x10 Molar mass of N2 kg mol K 1 0x10 Permeability M n 2 1x10 gt Fluid gas viscosity Pa s Pret 1 013x10 Reference pressure Pa 0 4 Macroscopic porosity 46 Table 2 Effective diffusion coefficients Name Value Description Units Pegg 0 0000074 Binary diffusion coefficient for O and N m s TD A 0 0000087 Binary diffusion coefficient for O and H O m s Diana 0 0000285 Binary diffusion coefficient for O and N m s Peng 0 0000080 Binary diffusion coefficient for O and N m s Table 3 Inlet parameters Oi 0 1 Mass fraction of H gt at anode inlet Oi po 0 9 Mass fraction of H gt 0 at anode inlet O55 0 168 Mass fraction of O gt at cathode inlet Oi po 0 2 Mass fraction of H gt 0 at cathode inlet Ok 0 632 Mass fraction of N at cathode inlet SX 0 5 Molar fraction of H at anode inlet ER 0 5 Molar fraction of H2O at anode inlet xd 0 1348486 Molar fraction of O gt at cathode inlet Kn 0 2853939 Molar fraction of H gt 0 at cathode inlet Ks 0 5797574 Molar fraction of N gt at cathode inlet p 111430 Pressure at the anode inlet Pa ps 111430 Pressure at the cathode inlet Pa Table
54. you can change the source file extension When you run zinc the presence of the DLL file is detected and in this case cathodetoken dll is copied into the Zinc directory Zinc then runs linked to that DLL Path DAMGC incprog examples fuelcell Input files Input cathodetoken zin Mesh file cathodetoken mtf Constants file cathodetoken con NL file cathodetoken f90 Restartfile None initial condition supplied Output files Outputfile cathodetoken zou List file cathodetoken zls Status READY Open R Cancel module nlinc private I public Cfun afun ffun gfun gfun BCfun scanfun integer save init 1 double precision save R T m1 m2 m3 kappa eta Pref epsilc double precision save D1 2tilde D1 3tilde D23tilde D21tilde C double precision save x1c x2c x3cpcRT ml13 m23 end module nlinc subroutine NLinit I If the first call read from cathode con use nlinc implicit none lt Fig 3 Zinc solver ready to run the fuel cell problem Here the non linear source code routines have the extension f90 for a Fortran file You can change this to match your language 55 Having run the simulation it is now time to post process with Zpp Fig 4 Here we stick to simple output of the raw quantities the mass fractions and the pressure Fig 5 shows the simulation result viewed with Zmesh Zinc 2 0 Zpp XI Path DAMGC zincprog examples fuelce HFINAL ISTEP 30000 Inputfiles 0 00000E 00 0

ZINC Tutorial Manual

Contents

Download Pdf Manuals

Related Search

Related Contents