Next: Eigenstrain programs
Up: Actual program operation
Previous: Actual program operation
The following gives some actual details of use for the finite element programs ELAS2D.F, ELAS3D.F, ELECFEM2D.F, and ELECFEM3D.F. The eigenstrain programs have some differences in them, and will be discussed separately. The 2-D and 3-D programs are very similar to each other, for both the elastic and electric cases, so that the discussion below assumes 3-D.
The first task is to assign the system size, by choosing the values of nx, ny, and nz. These can be different from each other. The total number of nodes is, in the program, ns = nx x ny x nz. The comments in the finite element programs show clearly which arrays need to be dimensioned by the system size = total number of nodes. These variables can be dimensioned at the same time with a global replacement, since the dimensions must be changed in each subroutine as well as the main program.
The next task is to generate a digital-image structure in a separate program,
or in the subroutine PPIXEL. Define i, j, and
k as the indices running in the x, y, and z directions in a right-handed
Cartesian coordinate system. Once the structure has been generated, for an image
that is nx x ny x nz, the finite element programs assume it is
stored in a file that
has been generated by the following piece of code:
DO 10 K = 1, NZ
DO 10 J = 1, NY
DO 10 I = 1, NZ
M = NX*NY*(K - 1) + NX*(J - 1) + I
WRITE(8,1) PIX(M)
1 FORMAT(i1)
10 CONTINUE
and reads the file in the same way in subroutine PPIXEL. In the programs, this file is called MICROSTRUCTURE.DAT (the name of which can be set by the user). The two OPEN statements in the beginning of the programs give the file names for the microstructure input file (unit=9) and the results output file (unit=7). The values in pix(m) are phase labels, starting at the value one, i.e., phase one has label 1, phase two has label 2, and so on, up to the maximum value of nphase, which must be set by the user. If the user wishes to generate a simple structure within the program, the place to do it is in subroutine PPIXEL.
The value of gtest is also assigned, and is used for determining when the gradient of the elastic energy is small enough. One usually determines this using solutions of exactly known composite problems (see Section 5), or operating directly on the problem of interest and checking how the various calculated quantites, like effective elastic moduli (conductivity), local stresses (currents), etc. change with relaxation. The value of gtest should be determined by numerical experimentation, for the specific problem of interest. The only way to be sure that gtest is small enough is to print out the values of the quantity desired and see how many conjugate gradient cycles are enough so that this quantity is no longer changing. If gtest is picked to be too small, much CPU time can be wasted in relaxing a system beyond the point at which the desired answer no longer changes significantly.
Material properties have to be assigned to each distinct phase. For the elastic programs, in the isotropic case bulk and shear moduli are assigned in the main program. The variable containing the bulk and shear moduli for each phase is phasemod. In the anisotropic case, the entire elastic moduli tensors must be inputted in subroutine FEMAT. In the electric programs, the entire conductivity tensor must be supplied in the main program, for isotropic or anisotropic cases. In ELAS3D.F, the values of the applied strain are assigned (exx, eyy, ezz, exz, eyz, and exy) by the user, and in ELECFEM3D.F, the components of the applied electric field are assigned (ex, ey, ez) (similarly in 2-D).
Sometimes more than one microstructure will be considered in the same program, or the same microstructure but with several different sets of properties will be considered. In this case, the parameter npoints can be set to the appropriate value of iterations (micro = 1, npoints). One then gives an initial set of displacements (voltages). In ELAS3D.F (ELECFEM3D.F), the initial displacements (voltages) are obtained by assuming that the applied strain (field) is the same everywhere. It is sometimes advantageous, in the case when several sets of phase properties are used on the same microstructure, for the initial displacements (voltages) of one computation to be the final values of the last computation. This can save CPU time, and is easy to implement, by only going through the displacement (voltage) initialization if micro = some fixed value.
The subroutine FEMAT is then run, and sets up: the stiffness matrix dk for each kind of pixel, the elastic moduli tensors cmod, and the various linear and constant terms needed because of the periodic boundary conditions. Subroutine DEMBX then runs the conjugate gradient algorithm to find the answer to the problem that satisfies the gtest criterion. The user must supply the values (or use the default values) of two parameters, kmax and ldemb, that control how dembx works. The value of kmax controls how many times DEMBX can be called in a given computation, and ldemb tells how many conjugate gradient iterations DEMBX will perform at each call. These two parameters can be adjusted so as to be able to see how the relaxation is going. Each time DEMBX returns to the main program, intermediate values of the average stress (current), the energy, and the value of gb x gb are printed out, to serve as a check on how the relaxation is progressing. If the gtest criterion has not been met, DEMBX is then re-entered to perform the next ldemb conjugate gradient cycles. When relaxation has finally met the gtest criterion, subroutine STRESS (CURRENT) is then used again to output the strain and stress (field and current) at every pixel, in order to compute some kind of average or produce a map of local quantities.
Next: Eigenstrain programs Up: Actual program operation Previous: Actual program operation