Computer Graphics & Geometry

Finite Element Data Processing Based on Implicit Complexes

E. Kartasheva, A.Boldarev, S. Boldyrev, G. Bagdasarov, S. D'yachenko,
V. Gasilov, O. Olkhovskaya, V. Shmirov

Institute for Mathematical Modeling, Russian Academy of Science , Moscow , Russia ,

E-mail: ekart@imamod.ru

Abstract: Implicit complexes (ICs) provide a novel framework that makes it possible to model heterogeneous objects exploiting hybrid representation schemes. The paper is focused on the application of IC based heterogeneous models to designing data processing tools for finite element/finite difference (FE/FD) analysis. We give a brief description for implicit complexes, their structure, basic operations and implementation. Then we represent our approach to designing FE data structures within the IC framework. As an example we consider software architecture of a radiative magneto-hydrodynamics simulation code MARPLE developed in IMM RAS for numerical solution of research problems in the field of gas-dynamics and plasma-dynamics.

Introduction

In this paper, we discuss problems concerning design and implementation of data structures for FE/FD (finite element/finite difference) analysis.

Numerical FE/FD methods are based on discrete representation (surface and volume meshes) of a computational domain, although meshfree analysis and simulation methods are also emerging. Meshes are also actively used now in visualization, animation, computational geometry, image processing, and other areas. However, requirements for discrete models in FE/FE analysis (FEA) are stricter than in other areas. Such requirements are formulated in terms of discrete model topology and metric characteristics, exactness of the geometry approximation, and conformity with initial attributes. This results from the use of meshes in FEA for approximation of systems of equations of mathematical physics. The size and shape of mesh elements and mesh structure seriously influence the stability of numerical simulation procedures and accuracy of obtained solutions.

As the experience shows, using of index coordinate grids is not worthwhile in domains of complex geometry shape (especially non-uniformly scaled), when a significant difference in physical parameters exists or an accurate presentation of near-boundary processes is of interest. That is why it is necessary to use combined meshes with irregular structure. While handling such meshes, a number of geometrical operations are required: monitoring of the mesh structure and various relations of its elements, calculation of geometrical characteristics of mesh elements (i.e. segment length, angles, cell squares etc.), installing correspondence of various objects with the mesh. In addition, a possibility should be provided for modifications of the grid during computations and using in calculations different types of incidences (node-edge, node-cell and others). thereby it is advisable to have a specialized geometrical system as a part of the numerical simulation code. This module should perform data preparation and preprocessing with extended functions including creation of the geometric model of physical simulation domain, mesh generation, optimization and adaptation, as well as the preparation of initial and boundary conditions. Such a geometrical system has to provide tools for supporting data structures used by both the data preprocessor and the solver.

In this paper we consider an approach to designing the required FE data preparation system based on hybrid cellular-functional geometrical models represented by implicit complexes (IC). Object oriented programming is used as the most suitable implementation technique for such an approach.

The structure of the paper as follows. Initially we discuss main FE data preparation problems, review related works and explain the motivation of our approach to development of FE data processing tools. Then we describe main conceptions of IC based geometrical model. After that we describe IC model implementation and discuss our approach to solution of the FE data preparation problems in the frame of IC based models. Finally, a radiative magneto-hydrodynamics (RMHD) simulation code MARPLE is described. MARPLE data preprocessor is considered as an example of realization of the concepts considered in the paper. A few examples of practical problems solved using MARPLE tools are also represented.

1.     FE data preparation

The main goal of FE/FD data preparation is to provide a discrete representation of the computational domain with assigning boundary conditions and initial values distributions. Let us consider a formulation of the discretization problem in terms of initial and resulting objects along with a set of special requirements.

Generally, a computational domain can be represented as a heterogeneous geometrical object having an internal structure with non-uniform distribution of material and other attributes of an arbitrary nature (photometric, physical, statistical, etc.), and consisting of elements of different dimension. Attributes are used for representation of initial data and boundary conditions of the simulated process. Attributes can also provide a formal description of requirements and constraints imposed by FE methods, for example, constraints on the size of mesh elements.

Thus, given the initial object, it is necessary to build a resulting heterogeneous object with discrete  geometry. In conventional terms, such a discrete model is called a mesh. The mesh provides an approximation of the initial object geometry. The resulting mesh depends not only on the initial geometry  but also on the initial attributes (in particular, on the mesh density attribute). Properties of the object after the discretization are described by attributes associated with the mesh. Note that the resulting attributes can differ from the initial ones in terms of their number, set of values, and description. Some of the initial attributes may not have direct counterparts in the discrete model. Other attributes associated with the mesh can have the same meanings and similar values as their initial counterparts but may be described by another representational scheme. Completely new attributes can also appear in the discrete model. For example, such an attribute can describe normals to the initial implicit surface at the nodes of a mesh or can represent values of 3D cell volumes evaluated in advance and useful for speeding up FE-based calculations. So, in general the resulting attributes are created depending  on the initial attributes as well as on the initial geometry 

Mesh generation oriented towards FEA applications  should  take into account the following requirements and constraints upon the resulting discrete structures (surface and volume meshes):

1. The boundary mesh must provide an adequate approximation of the underlying boundary surface of the initial object. To this end, it is necessary to bound the distance between the mesh and the initial surface and to limit the deviation of normal vectors to the mesh from those to the surface. The surface mesh  has to include all the sharp features of the initial object's surface. This means that there should be conformity between initial object ridges and peaks, and edges and nodes described by the boundary mesh;
2. There can be specific constraints upon the size, the shape of mesh elements and upon proportions of adjacent cell sizes in the mesh. These constraints are useful for ensuring a reliable discretization and convergence of numerical simulation procedures.
3. Non-geometrical attributes reflecting some features of FE modelling can also result in specific constraints upon surface and volume meshes. Let us describe the following cases:
1. Given the attribute function, one can set the constraint concerned with the limitation of this function's variations within a mesh . Such a function can describe some estimated characteristics of the process being modeled or just error estimations. Attributes representing material or medium properties can also impose constraints of this kind.
2. Then, while decomposing the object  one should take into account that attribute function features, namely singular points, have to coincide with mesh nodes, and lines/surfaces of discontinuity have to be exactly described by 1D/2D elements of the resulting mesh. For example, such requirements are necessary when a computational domain is decomposed into sub-areas corresponding to different materials or different properties of medium. Singular points and lines can also conform to disposition of sources or to load application areas.
3. If no special constraints are formulated, the number of elements (cells) in meshes should be minimal.

Some of these requirements can be satisfied during mesh generation. For example, the adaptive advancing front method [FG00] allows to constructing meshes compatible to a predefined cell sizes distribution.  But in general special optimization procedures are applied to initially constructed meshes for their adaptation to all the FEA requirements.

Thus we can summarize all the phases of the FE data preparation process: initial model description, attribute assigning, volume mesh generation according to predefined constraints, subsequent remeshing for adaptation to special FEA needs; converting initial attributes for getting new ones associated with the resulting computational mesh.

Analyzing the listed above processes we can conclude that FE data processing tools should support heterogeneous geometrical models of different types and provide functionalities for data conversion and discretization of the models.

Then we give a review of works in the areas that are relevant in the context of our consideration and provide motivation of our approach to FE data processing tools design.

2.     Related works

We briefly discuss some previous works on modeling dimensionally heterogeneous objects, objects with varying distribution of material and other attributes, and approaches to combining various geometrical representations. We also consider the known FE mesh generation methods.

2.1     Heterogeneous objects modeling

Particular attention in solid modeling is paid to modeling heterogeneous objects with multiple materials and non-uniform internal material distribution [CAD05].  

Such objects may represent mechanical parts or assemblies, geological and medical models, the results of physical simulations as well as time-dependant models of artistic nature, and are heterogeneous in terms of view of their internal structure and their dimensionality. Such objects may represent mechanical parts or assemblies, geological and medical models, the results of physical simulations as well as time-dependant models of artistic nature, and are heterogeneous in terms of view of their internal structure and their dimensionality. Both topological subdivisions and constructive procedural methods are used for description of heterogeneous objects.

Various topological stratifications and complexes are used to describe spatial subdivisions [R94, MRG00]. Types of stratifications and complexes differ in the constraints imposed on the topology of their elements and on the connections between these elements. Multidimensional simplicial complexes are used in [PBCF93] for dimension-independent geometric modelling for various applications. A Selective Geometric Complex (SGC) [RC90] is a non-regularised non-homogeneous point set represented through enumeration as the union of mutually disjoint connected open subsets of the real algebraic variety. A SGC provides a framework for representing objects of mixed dimensionality possibly having internal structures and incomplete boundaries.

The Djinn API for solid modelling [ABC00] is based on objects partitioned in a cellular way and containing mutually disjoint cells which are manifold point-sets of differing dimensionality in 3D. In [GH95] an extension of B-spline surfaces to surfaces of arbitrary topology is proposed. Polyhedral complexes are used to describe the surface topology. A procedure for designing cellular models based on CW-complexes with an emphasis on the topological validity of the resulting shapes is considered in [K99,OK01]. Selective Nef complexes were proposed in [HK05].

The work on constructive topological representations [RAG06] introduces the stratified structure that is quite different from the topological complexes. This stratification is defined on n-dimensional solids using ‘natural' topology based on the neighborhood concept. It considers only the n-dimensional atoms and ignores the lower dimensional ones as well as the connectivity characteristics of the atoms and of the corresponding solid. Such a model can not be used for describing heterogeneous objects containing components of different dimensionalities.

To specify non-geometric properties of objects, spatial subdivisions are also used in computer graphics and in finite element analysis (FEA) as the underlying structures for piecewise analytical descriptions of attribute functions. Usually a basic topological subdivision is selected, which can be described by a topological stratification [RC90, ABC00], a cell complex [CD*02, KBDH99], or a voxel model [CT00]. Different types of functions can be used to describe attributes [JL99,PCB01, MC01, ST02].

Constructive approaches are applied to description of heterogeneous objects consisting of multiple materials and combining different geometrical representations. In the STC framework [R97], a composite object is defined using a combination of layers each of which is described by a geometric complex, which is homogeneous with respect to the representations of the components. A model for objects with fixed dimensionality and heterogeneous internal structure (multidimensional point sets with multiple attributes or constructive hypervolumes) was proposed in [PAS01]. This model supports uniform constructive modelling of point set geometry and attributes using real functions of point coordinates.

A feature based design methodology is represented in [QD04]. Under this methodology, heterogeneous objects consisting of multiple materials are constructed from engineering significant high-level components called form features and material features. Form features describe the shape of the objects and material features are used for modeling material variation. The relationships between form features and material features in heterogeneous objects are examined in [QD04] along with constructive feature operations. Set-theoretic operations for modeling functionally graded materials associated with a BRep geometry model are discussed in [SD01].

Constructive Volume Geometry (CVG) [CT00] combines geometry and attributes in a systematic manner. The model is presented as algebra of 3D spatial objects utilizing voxel arrays and continuous scalar fields for representing both geometry and photometric attributes (opacity, color, etc.). A distance field based approach for heterogeneous object modeling, in which the space is parametrized by distance to the geometry boundaries is proposed in [BST02].

The HybridTree [AGCA06] is a constructive tree with leaves defined by a number of representations. Both the function evaluation and the surface mesh generation are provided for modelled objects. Depending on the user's query, corresponding conversions between representations are applied. A hybrid constructive tree in [FG*05] has leaves with both implicit and parametric representations. To polygonize the surface of a complex object, surface meshes of primitives are classified against the subtree defining function, trimmed, and then merged into the resulting mesh. However, both of these approaches do not support heterogeneous objects with components of different dimensionalities and do not provide a description of the topological structure of the object being modelled.

A hybrid cellular-functional model based on the notion of an implicit complex, which provides a valid topological description of heterogeneous objects and allows for the flexible combination of cellular and functional representations of the geometry of objects and their attributes, was introduced in [AK02].

In papers [AK02, KAC*05] an IC's general structure and some basic procedures of the IC-based model construction along with suitable rendering methods were represented. Work [KAC*08] provides a systematic description of the IC-based hybrid model, algorithms for construction and discretization of ICs as well as the set-theoretic operations on heterogeneous objects represented by ICs are considered.

In this paper we propose to apply the heterogeneous model based on IC to solving FE data preparation problems.

The IC-based framework provides a unified description of heterogeneous object's geometry, topology, and attributes. An object is described as the union of cells of various dimensionalities, geometric representation types, and non-geometric attributes along with the relations between cells. Hybrid models constructing within the IC framework is not supposed to serve just for the combined usage of separate representations but is genuinely unified.

The IC-based framework allows for representing a heterogeneous object by a union of high-level components that are significant for a given application. For example, such components can describe mechanical parts or parts of an experimental installation. The components can be overlapped but special constraints on the description of the mutual dispositions of the components are introduced. Thus the IC-based framework provides the representation with a distinctive structure that can easily be reduced to the cellular topological subdivision. The intersections of the components are described by the constructive methods which preserve the precision of the representation. The representation can be quite compact if it involves only those entities which are necessary for descriptions of the initial components and their mutual dispositions.

The main relations characterizing mutual locations of cells are the boundary and the "to contain" relations. Non-geometric attributes are independently described by functional or cellular models and are associated with IC's cells by means of attribute relations. IC based heterogeneous representations provide all information necessary for numerical simulation, for contact and friction analysis and for other CAD/CAM applications.

2.2     Finite element mesh generation and refinement

Issues of finite element mesh generation are discussed in detail in [FG00]. Tetrahedrization is one of the widely used methods of 3D discretization. The main approaches to automatic tetrahedral (triangular) mesh generation include spatial decomposition based methods, Delaunay type methods, and advancing-front techniques. Algorithms based on spatial decomposition are relatively easy to implement, but they do not allow for detection of boundary sharp features and cannot distinguish boundary entities which are rather close but not directly connected. The boundary connectivity constraint is not taken into account in Delaunay tetrahedrization. So, local mesh modifications are necessary to fit the boundary.

More accurate boundary representation is supported by the advancing-front method. This method starts from a domain boundary discretization and marches into the region to be processed by adding one element at a time. However, since the method is based on local operations, convergence problems may be encountered. The convergence problem is common for all methods as there is no theoretical result which can guarantee that a polyhedron with the given boundary triangulation can be subdivided into tetrahedrons without adding internal points. In spatial decomposition methods, the convergence problems appear when refining small details and sharp features. For the Delaunay type methods, the convergence of the boundary fitting procedures has not been proven. We use the advancing-front technique as it is applicable to arbitrary solids and allows us to control shapes and sizes of tetrahedrons during the mesh generation process. A modification of this technique for increasing the effectiveness of the tetrahedrization procedure for FRep solids is described in [KAP*03].

Issues of surface mesh optimization for FEA are considered in detail in [FB98, FG00]. The described techniques are based on the consecutive application of different mesh simplification, mesh subdivision, and mesh adaptation procedures. Details of such operations are discussed also in [GH97, S01] and other works. Note that in these works the surface models are not defined in terms of analytical functions but rather by means of triangulation (resulted, for example, from measurements, CAD, biomedical engineering). During the mesh refinement, the exact definition of underlying surfaces is unknown. When optimizing polygonized implicit surfaces, we can use both approximate and precise functional surface models, which provide for more precise calculations of surface characteristics and corrections of the node positions in respect to the underlying surface for remeshing. Mesh optimizations methods for implicit surfaces discretizations are described in [KAP*03]. Iterative simplification and refinement of the implicit surface preserving sharp features are described in detail in [KAP*03]. 3D mesh optimization procedures are described in [RL92, LJ95, FO97].

3.     Implicit complexes.

In this section we provide a brief description of our approach to modeling heterogeneous objects based on the Implicit Complex (IC) notion.

We consider a hybrid model defined in the Euclidian modelling space  as follows. Let g_i be a closed point set called a cell. Then, a geometric object D is defined as the union of cells g_i  of different dimensions under the following conditions:

1. The boundary of each cell  is the union of a finite number of cells of lower dimensions;
2. Cells can overlap each other but the intersection of any two cells is either the union of a finite number of cells or is empty. (Note that we call the cells satisfying 1 and 2 conditions as properly joined cells.)
3. Each cell is unambiguously described by some known geometric representation which provides a set of tools for geometrically and topologically correct discretization of the cell. So, a variety of representations can be used for the description of the cell shapes. However, all of these representations should guarantee a conversion into a mesh described by a polyhedral  complex.

A collection K of cells satisfying the above conditions is called an implicit complex. The dimension of IC is the maximal dimension of its cells. Correspondingly the point set union of all cells of the IC K is denoted by |K| and called a carrier of K. Fig. 1 shows an example of an IC.

The above conditions actually ensure the ability to convert an arbitrary IC K into a polyhedral complex, which approximates, geometrically and topologically, the object D being modeled. In fact, the reducibility of each of the used representations to a polyhedral one guarantees a correct execution of any operation on objects described by various representations. However, for exploiting advantages of different types of representations the initial representations of components are kept and meshes are used only for the implementation of various numerical procedures. The support of overlapping cells allows for inserting components of a composite object into its IC model without subdivision. The IC definition conditions are satisfied by adding extra cells describing mutual intersections of the components. This allows for the preservation of the initial representations of components, which is useful for heterogeneous object modeling.

An implicit complex provides a consistent description of both the geometry and the topology of a modeling object. The geometry is represented by the geometry of the individual cells and the topology is described by means of the relations between cells.

3.1     The IC Topology

The main relations defining the topology of an IC are the boundary relation and the relation "to contain". According to the first two conditions of the IC definition the mutual disposition of any of the IC cells can be evaluated through queries to its main relations.

A boundary relation between p-dimensional cells and s-dimensional cells of an IC is denoted as Rb^ps, s<p. By definition the boundary relation Rb^ps consists of pairs of cells (g^p,g^s) where s-dimensional cell g^s  belongs to the boundary of p-dimensional cell g^p and does not lie in the interior of any other boundary cell of g^p. The relation "to contain" between p-dimensional cells and s-dimensional cells of an IC is denoted by Rc^ps, s ≤ p. A pair cells (g^p,g^s)  belongs to Rc^ps if  g^s  belongs to the interior of cell g^p.

For an unambiguous definition of a 3D IC, it is necessary to describe three boundary relations Rb={Rb¹⁰, Rb²¹, Rb³²} and nine "to contain" relations Rc = {Rc^ps}, (s ≤ p, p=1,2,3, s=0,1,2,3) for all cells of dimensions from 0 to 3. Other boundary relations can be calculated using the composition operation (denoted by a symbol '*'). For example, Rb²⁰=Rb²¹*Rb¹⁰. Thus the description of an implicit complex K consists of the collection  of cells and the sets of boundary Rb relations and "to contain" Rc relations.

Some additional relations can be useful for implementing those operations on the ICs that require a faster access to the information about the mutual disposition of cells. The most often used additional relations are the co-boundary, the "to be contained", the incidence and the adjacency relations. These relations can be derived from the boundary and the "to contain" ones using various operations on relations.

Note that in principle it is not necessary to have all the possible relations explicitly described and stored. As to queries to the relations, they can be evaluated using the described dependences between relations. If needed, the additional relations can be dynamically computed on the base of the main ones and then discarded again (if not needed afterwards).

One can define not only relations within the IC cells but also between different complexes. They describe mutual disposition of cells belonging to different ICs. The equivalence relation and the relations "to contain" , "to be contained" between implicit complexes are introduce in [KAC*08]. The equivalence relation Rq^s between s-dimensional cells of  ICs A and B consists of pairs equivalent s-dimensional cells (a^s,b^s), where a^s belongs to A and b^s belongs to B. The relations "to contain" and "to be contained" between cells of different ICs are introduced by analogy with the corresponding relations between the cells of a single complex.

3.2     The IC geometry

In [AK*02] a basic IC structure have been proposed. Work [KAC*08] introduces more elaborated structural classification in the form of five types of IC cells that differ in their geometric representations:

- The P-cell, which is an explicit cell representing a simple polyhedron.

- The B-cell, which is a cell representing a manifold defined by its boundary. B-cells describe segments of parametric curves, patches of parametric surfaces and boundaries defining 3D solids. A 1D (2D) B-cell is defined by its supporting curve (surface) and by its oriented boundary. A 3D B-cell is defined by its oriented boundary only. In the general case the boundary of a B-cell can consist of cells of all the other types supported in the IC framework.

- The F-cell, which is an implicit cell described by the FRep that is a constructive representation by real-valued functions [PAS95]. A valid variety of 2D and 1D FRep objects are restricted by s-dimensional F-cells (s<3) to those which are represented as subsets of the boundaries of 3D manifolds.

- The C-cell, which is a composite cell aggregating cells of various types. Each C-cell is defined as a carrier of an implicit complex T differing from the complex K containing this C-cell. The complex T can consist of the cells of all types supported in the IC framework. In particular case T can be a simplicial or a polyhedral mesh. The complex T is not a subcomplex of K. Its cells are not properly joined with respect of the cells of K.

- The T-cell which is a cell described by a constructive tree. Its leaves represent objects described by cells of all the other types. The tree nodes represent operations admissible for the IC - in particular, some bijective geometric transformations, non-regularized set-theoretic operations and trimming by 3D manifolds. The T-cells allow for the description of the result of applying set-theoretic operations to cells of different types without the need for converting between representations.

Various types of cells are illustrated by an example described in the following section.

In accordance with the IC definition, polyhedral and cellular complexes are particular cases of implicit complexes. There are additional constraints both on the geometry of the complex's cell and on the relations between them, namely all the cells of the polyhedral complex are polyhedra by definition, and the intersection of any two cells is either empty or equal to the intersection of the boundaries of these cells. So all the "to contain" relations of a polyhedral complex are empty and the topology of a polyhedral complex is described by the boundary relations only. Polyhedral complexes are very important for FE applications, as they actually represent meshes As to the cellular complexes, they differ from polyhedral ones only by the shape of the cells. Cellular complexes includes curvilinear edges, faces and volumes.

3.3     IC model example

Let us consider an example illustrating geometric and topological features of an IC based model with cells of different types and with basic relations between them. Fig. 1 shows a 2D object consisting of two 2D components (rectangle DEQS and disk LHKF) and one 1D component (segment OM ). We assume that the components are represented by different methods and have different non-geometrical properties (attributes). Rectangle DEQS  is described by the boundary representation, disk LHKF is defined functionally and segment OM is described as straight line segment. One of the possible IC representations of this heterogeneous model is described by the complex K consisting of the following cells:

·         2D cells: rectangle DEQS, disk LHKF , half-disk LKF;

·         1D cells: closed polylines DEQS, circle LHKF, arc KFL, segments KL, OM , ON ;

·         0D cell: points K, L, O, N, M.

H

Fig. 1 A 2D object consisting of two 2D components (rectangles DEQS and disk LHKF) and one 1D component (segment OM ).

According to the IC definition, the complex K includes the cells representing the initial components (DEQS, LHKF, OM ), their boundaries (DEQS, LHKF,O,M), and the cells (LKF; KFL, KL, OM , ON, K, L, N) describing the mutual intersections between the listed initial ones . We assume that initially the rectangle DEQS is defined by the boundary representation, the disk LHKF is described functionally and the segment OM is specified explicitly by its end points. Consequently we represent these components by the B-cell DEQS, the F-cell LHKF and the P-cell OM , correspondingly. Other cells are added to satisfy the constraints of the IC definition. For example, the 2D cell LKF represents the intersection of the initial 2D components. So it is represented constructively by a T-cell. 

Fig. 2. Graph representing the boundary relations for an IC of the object shown in Fig. 2

Fig. 2 shows a graph representing the basic boundary relations in the IC K. Here the graph nodes represent the cells of the IC K, and the graph edges represent connections between cells in the boundary relations.

The "to contain" relations for the IC K are represented by a graph shown in the Fig. 3. Here once again the graph nodes represent the cells of the IC K, and the graph edges represent the connections between the cells in the relations "to contain".

Fig. 3. Graph representing the relations "to contain" for the IC of the object shown in Fig. 1.

3.4     The Attribute Model

In this section we consider a cellular-functional representation of attributes associated with an IC according to work [KAC*08]. The attributes of the object and its geometry are independently described. Each attribute Ai is described by a set Ni of its values embedded into a multidimensional real number space R^mi of a proper dimension mi. The interpretation of attribute values depends on its nature and specifics. For the sake of uniformity, the value set Ni of each attribute is supplemented with a special "empty" value  and so all the attributes are defined at each point of the modeling space. An attribute associated with  an IC K is represented by a collection of attribute functions Si={Sij}_j=1,J and a set of attribute relations  Ri={Ri³,Ri²,Ri¹,Ri⁰}. The relations (Ri)^p  (p=0,1,2,3) associate functions of an attribute Ai with the cells of the IC K. The relation (Ri)^p  consists of pairs (g^p,Sij)  where g^p is p-dimensional cell of the complex K and Sij is an attribute function..

Each function Sij of an attribute Ai maps the modeling space into the attribute value set Ni. Attribute functions can be analytic, piecewise analytic or be defined by an interpolation methods [KAC*08]. Note that for defining interpolated attribute functions, one can use a space partition different from the one associated with the object. This means in the general case introducing another complex different from K. Note that it is possible to introduce various operations on attributes in the model. As attribute value sets belong to mathematical spaces R^m, all these operations are reduced to operations on real numbers and vectors. The attribute functions are defined on the entire modeling space. So we can construct new attribute functions on the base of th e k nown ones using various mixing functions.

4.     Operations in the IC framework

4.1     The basic operations in the IC framework

Here we listed basic operations over cells and entire implicit complexes. These operations are especially important in the context of constructing and manipulating implicit complexes. For each operation, there are constraints on input data that should be checked before their actual evaluating. More detailed description of the operations introduced within IC framework is given in [KAC*08].

Cell_adding adds a new cell to IC. Cell g can be added to IC K if g is properly joined to all the cells of K and the boundary of g  is represented as the union of cells of K. The input data of the procedure include the cell geometry description and the list of cells of K related to the cell being added.

Attribute_adding defines an attribute on K. The input data include a collection of attribute functions and the list of cells of K being associated with this attribute.

Cell_removing deletes a cell from IC. It is follows from the IC definition that cell g  can not be deleted if one of the following conditions or both of them are satisfied

1) Cell g  has co-boundary cells in K.

2) Cell g  represents the intersection of some other cells of  K. These conditions are checked automatically using IC relations.

Cell_cutting removes a cell with all its boundary cells from IC. It is implemented under the same restrictions as the Cell_removing operation.

IC_adding implements sum operations over properly joined complexes. . Two implicit complexes are called properly joined if their cells altogether satisfy conditions 1 and 2 of the IC definition. Given any two properly joined implicit complexes C and T, the sum of C and T is a complex M consisting of all the cells of the complexes  and , this is denoted as M=C+T. The boundary relations and the relations "to contain" of M are automatically formed on the base of the same relations defined on the complexes T, C. Attributes are established on the IC M by combining those ones defined on the initial complexes C and T using an appropriate mixing function.

IC_transform implements bijectivegeometric transformations over ICs including affine transformation and nonlinear bijective transformations.

4.2     Set-theoretic operations over objects represented by IC.

Set-theoretic operations are the main mechanism for constructing composite geometric objects starting from more primitive ones. In the IC framework, the union and intersection operations over ICs with attributes are introduced. Given two implicit complexes A and B, the intersection of A and B is an implicit complex Ñ whose carrier is equal to the set-theoretic intersection of the carriers of A and B, C=AÇB. The union of two implicit complexes is a complex whose carrier is equal to the set-theoretic union of the initial complexes. Thus the non-regularized set-theoretic operations are considered. The difference operation is more problematic as a mere set-theoretic difference of the carriers results in non-closed objects. So the restricted version of the difference operation are introduced, namely, trimming with a 3D manifold. Let the IC B represent a 3D manifold; then the result of trimming the IC A by the complex B is a complex C whose carrier is equal to the set-theoretic intersection between the carrier A and the point set described as the inversed carrier of B (that is the cavity in the whole solid space). This operation is denoted as C=A-B. Realization of set-theoretic operations over ICs is considered in detail in work [KAC*08].

Four procedures for the implementation of the set-theoretic operations are introduced in [KAC*08]. Each of them takes two input ICs A and B and returns the IC C. The intersection and the trimming operations are implemented as the procedures IC_intersection and IC_trimming, correspondingly. There are two procedures implementing the union operation, namely, IC_union and IC_subtractive_union. The procedure IC_union returns the IC C involving all the cells of the input ICs A and B and is realized as. C=A+(AÇB)+B. The procedure IC_subtractive_union calculates the union of the ICs A and B according to the formula C=A+(B-A). Thus the resulting IC C does not include those cells of IC B which lie inside the carrier of A. 

4.3     Discretization of IC based models.

The IC_meshing procedure implements the conversion of the IC representation into the simplicial one. The discretization of IC models are guaranteed by the IC definition according to which the appropriate mesh generation methods have to be available for all types of the IC cells. The discretization of an implicit complex K is implemented as an iterative process. The mesh generation of each kind of cell is implemented using specific meshing algorithms. An extensive survey of discretization methods are presented in the book [FG00]. 

We subdivide IC cells in the order of their increasing dimensionality. Among the cells of the same dimensionality we first subdivide those ones which do not contain other yet unprocessed cells. Thus, at the moment of the meshing of a cell we already know the discretization of its boundary and the subdivision of all the cells lying inside the considered one. Then we subdivide the cell into mesh elements such that they are compatible with other meshes which already belong to it. The corresponding incremental mesh generation approaches that allow for preserving existing mesh elements can be found in works [AGCA06], [KAP*03] and in the references contained in these works. The most incremental mesh generation  algorithms are based on the advancing front technique. This technique allows to subdivide multi-component domains with holes according to a predefined mesh density function. For illustration let us consider the application of IC-meshing operation to the IC shown in Fig1. The corresponding meshes are shown in Fig.4. According to the described meshing algorithm initially all the nodes K,L,O,N,M  are included into the mesh and then the other cells are subdivided in the following order:

1. 1D cells KL, ON, KFL  (where KFL is subdivided taking into account contained point O)
2. 1D cells OM (compatible to the subdivision of contained 1D cell ON);
3. 1D cells DEQS (compatible to the subdivision of 1D cell KL and to point O);
4. 1D cells LHKF ( compatible to the subdivision of contained 1D cell KFL);
5. 2D cell LKF;
6. 2D cell LHKF (compatible to contained 2D cell LKF);
7. 2D cell DEQS (taking into account the subdivision of the contained 2D cell LKF and 1D cell ON).

The three last steps of the meshing procedure are shown in Fig.4

   

Fig. 4. The sequence of the 2D cells discretization steps for the IC containing overlapped cells

5.     IC based data processing

Here we describe our implementation of IC based models and consider in detail all phases of the FE data preparation process from initial object description to final mesh generation and its adaptation to special FEA needs.

5.1     IC based data model implementation

We develop software for IC based modeling of heterogeneous objects within an object-oriented framework. Let us outline here only the principal classes which are directly derived from the presented theoretical description.

Data structures used for descriptions of topological subdivisions of objects or their boundaries are typically represented by the adjacency graphs whose nodes correspond to vertices, edges, faces, connected volume regions or combinations of these, and whose links capture information related to adjacency, orientation, and ordering [R94]. Instead of the frozen combination of the adjacency graphs optimized for a narrow range of applications we have elaborated more flexible data structures based on the concept of relations developed in discrete mathematics and widely used in computer science. Relations provide a unified description of various connections between individual cells, collection of cells, and entire complexes. They are also used for assigning attributes. The use of well known operations on relations makes it possible to change the data structures dynamically thus adapting them to specific applications. Relations based tools allow us to realize various operations on ICs in an flexible and compact form. A flexible mechanism of dynamically creating and removing relations provides an IC program implementation that is effective for design of complicated assemblies as well as for working with cellular complexes and meshes.

We suppose that ICs can change in the process of calculations. Data assigned to IC's elements are not predefined  and can vary from one application to another so we introduced  the conception of a numerated set. We use such a set for representation of a collection of cells of complexes.  We imply that all the cells are enumerated and each of them has a unique number. So, having the number of a cell, it is possible to find all the data assigned to this cell. Usage of the known event handling technique allows us to implement changes in all the data structures when some cells are added to or deleted from a complex. 

The basic SetNumbers class represents a numerated set. It supports information about the size of the set of numbers, contains references to objects dependent on the numeration and initializes the corresponding event updating process for any changes of the numeration.  All classes dependent on numerated sets are inherited from the Dependent class. The following classes: Relation, ComplexRelations, AttributeArray and ArrayOfNumbers are dependent on numerated sets of ICs' cells. The Relation class contains all the pairs of numbers of related cells. The operations of the Relation class allow us to get the indices of all the related cells as well as to add and delete pairs of cells. The Relation class can describe connections between cells of the same complex as well as between cells of different complexes. The AttributeArray class is used for description of attributes associated with cells. This class is parametrized by the type of the attribute and depends on a numerated set of cells. The AttributeArray class assigns values  of the predefined attribute type to all the cells of the tied numerated set. The ArrayOfNumbers class  contains an arbitrary collection of  numbers of the numerated set which this class is linked to.

Object-oriented programming technology allows independent realization of the main components of the IC based heterogeneous model (IC's topology, geometry and non-geometrical attributes)

The basic IComplex class represents the topology of an implicit complex. It contains four SetNumbers objects for description of four numerated sets of cells of different dimensions and Relation objects for representation of connections between cells. The Relation objects are linked to the corresponding SetNumbers objects. Only the base relations (namely, three boundary relations and nine "to contain" relations) are permanently supported in the frame of the IComplex class, while other relations are calculated automatically by means of private methods of the ICcomplex class. The methods of the IComplex class realize the operations on cells of the complex as well as operations on ICs. The IComplex class includes operations for modifying the relations as well as inquiry operations on the relations between the IC's cells.

The geometry of 0-dimensional cells is described by the GVector class representing point coordinates. The geometry of 1D, 2D and 3D cells is specified by classes inherited from the base abstract classes EdgeShape, FaceShape and VolumeShape correspondingly. Each of these abstract classes describes virtual operations for defining the point membership as well as for rendering and discretization of the cells. Classes inherited from the base classes are designed for representation of particular shapes.

The Geometry class represents the geometry of a complex. It contains four objects of the AttributeArray class template instantiated for the references to the GVector, EdgeShape, FaceShape and VolumeShape data types correspondingly.

The GeomIC class is designed for description of a geometrical implicit complex. It is inherited from the IComplex class and contains the Geometry object for the cells' shape representation.

The ICattribute class represents an attribute data structure. It contains attribute relations represented by objects of the AttributeArray class template instantiated for appropriate attribute data types. Attribute values are described by classes inherited from the abstract class ICFunction.

We implemented a few different attribute types including the integer type (marker), space distribution types representing functions dependent on space coordinates and grid space distribution types describing functions as an interpolation of the values defined at nodes of a grid.

Heterogeneous objects are represented by the ICModel class inherited from the GeomIC class.  The ICModel class contains a list of ICattribute objects for description of non-geometrical attributes. The number of attributes associated with the model is unlimited. Additional AttributeArray objects can be defined outside the ICModel for representation of various physical data calculated during the process of numerical simulation. The class template is instantiated for application data types. The AttributeArray dependence on the numerated cells sets provides connection between the model and the associated attributes.

The ComplexRelation class is implemented for description of the relationship between different complexes. The ComplexRelation class contains two references to the tied IComplex objects and Relation objects for representation of the equivalence relations and the relations "to contain", "to be contained" between complexes.

IC based models can be used for representations of geometrical complexes of different types. But for more effective implementation of C++ code we realized separate classes for polyhedral and cellular complexes which are particular cases of ICs. The TComplex class is an analog of the IComplex class but the TComplex describes the topology of complexes consisting of non-overlapping cells so it does not include objects for storing the relation "to contain". The Grid class and the GeomComplex class describe geometrical complexes consisting of non-overlapping cells. Each of them is inherited from the TComplex class and contains an object of Geometry class. The Grid class represents polyhedral complex whose geometry is unambiguously defined by nodes coordinates. So another shape attributes of the Geometry class are not used inside the Grid class. Heterogeneous models based on cellular complexes consisting of non-overlapping cells are described by the GridModel and ComplexModel classes inherited from the Grid and the GeomComplex classes, correspondingly. The GridModel and the ComplexModel classes additionally include ICattribute objects for description of non-geometrical attributes.

5.2     FE data preparation within the IC framework

ICs provide the background for representation of heterogeneous geometrical models combining various representations and describing objects with arbitrary distribution of non-geometrical properties. Particular cases of ICs, polyhedral and cellular complexes, are suitable for representation of structured and unstructured grids of polyhedral or curvilinear elements with attributes. So FE data preparation problems can be solved by means of model converting tools inside the unified IC based software environment.

Now we can formulate the discretization problem as the problem of conversion of the initial IC based heterogeneous object D into the object G with grid based geometry. D is described by an object of the ICModel class and G is represented by an object of GridModel class. The geometry of the initial object can be constructed using the functionality of the GridModel class or can be imported from some other 3D modelers (for example, from CAD systems). Attributes of D object are used for description of initial and boundary conditions of a solving application problem.

Given the initial object D of the ICModel type, we are going to build the resulting heterogeneous object G of the GridModel type. Here, the geometry component of G is an approximate discrete representation of the initial geometry G. Conversion of the initial model into a grid is implemented using mesh generation methods. Constraints on the size of mesh elements are expressed through the mesh density attribute A_r whose map function defines a proper element size at each point of the modeling space. Such an attribute can be given by the user or can be calculated on the basis of other given attributes. There is a promising way of defining A_r  through functional representation on the basis of so called "sources" [L97, AK*02].

In the process of mesh generation we also form the relations between the topology complexes describing the initial model and the resulting grid. These relations are represented by an object of GridRelation class. Cells of the resulting grid inherit attributes from their preimages in the initial model. Thus, using the relations between the objects we can assign values of the initial attributes to the cells of the mesh. Additional attributes  can be also introduced on the resulting grid model.

In principle, there are two main strategies to discretize a heterogeneous object with generation of a volume mesh. The first strategy involves decomposing an initial 3D object into 3D elements (tetrahedrons, blocks, prisms, some combined volume mesh) that are optimised under FE requirements and constraints described above. Here, the boundary mesh  appears as a side effect of the 3D initial object decomposition. Another approach implies that first we decompose the surface of the initial object thus yielding the surface mesh. Then, this surface mesh is subjected to optimisation and refinement to make sure that it satisfies all the possible requirements following which one can build a volume mesh conformable to the refined surface mesh.

Our IC based data structures are suitable for realization both the approaches to 3D mesh generation. But, we follow the second approach, because most of the constraints and requirements deal with the boundary mesh  whose quality is crucial in the context of FEA. As some of the constraints contradict each other, it is important to ensure that all the accessible iterative optimisation procedures are performed to provide the best possible result. In addition, it is known [FG00] that some effective methods of volume mesh generation are actually based on boundary descriptions of computational domains.

So, we decompose the discretization problem into two relatively independent sub-tasks:

·         generation of a mesh of the object surface along with its refinement;

·         generation of a volume mesh conformable to the boundary one.

This approach agrees with the general IC discretization  technique described in the previous sections. So initially we subdivide 1D and 2D cells of the complex describing the initial object D and generate a boundary mesh B represented by an object of Grid class. Boundary subdivision is implemented taking into account the initial geometry and the mesh requirements described by the mesh density attribute. Altogether with the boundary mesh we create an object of GridRelations class describing connections between the boundary mesh and the initial object. Then we can apply mesh optimization procedures to the boundary mesh with simultaneous updating the grid relations. Note, that using the relations between the grid and the initial model we can get information about underlying surfaces during  the mesh refinement  process and introduce new nodes lying exactly on the boundary surfaces.

After the boundary mesh generation we construct volume grid that conforms with the boundary one and satisfy mesh density requirements. We form the relations between volume and boundary grids and between the volume grid and the initial model. They can be used during volume mesh optimization processes.

6.     Software architecture of a radiative magneto-hydrodynamics simulation code MARPLE

6.1     Scientific software.

The theoretical level of the power engineering problems investigations may be improved drastically at present time due to burst out progress of high performance computing systems enabling more realistic simulations. Real-world mathematical models are very complicated. They include the description of numerous physical and other processes and many variable parameters. In general the multiphysics models are so sophisticated that should be implemented only on the basis of parallel or distributed TFLOPS+ performance computations. This is true for both the industrial power engineering problems and the forward looking basic research, particularly in the field of plasma technologies. The team in the Institute for Mathematical Modeling, Russian Ac. Sci. (IMM RAS) in collaboration with colleagues from Kurchatov Institute and Troitsk Institute For Innovation & Fusion Research and other research centers proposed new computer models of high-speed radiative magnetohydrodynamic processes and developed novel problem statement data representations especially efficient for implementation of complex models. These theoretical and practical studies resulted in creation of application scientific high-temperature hydrodynamics and magnetohydrodynamics software. A number of urgent problems in pulsed power and controlled fusion were investigated using these codes. The suggested techniques allowed numerical validation of some innovative schemes for high-temperature plasma production with the help of super high electric power generators.

Modern problems in pulsed-power energetic issue a real challenge to the computer simulation theory and practice. High-performance computing is a promising technology for modeling complex multiscale nonlinear processes such as transient flows of strongly radiative multicharged plasmas. An essential part of such numerical investigations is devoted to computer simulation of self-constricted discharges, or pinches, resulted from electric explosion of cold matter, e.g. gas-puff jets, foam strings, or metallic wire arrays. These investigations were significantly stimulated by impressive results in the soft X-ray yield produced by wire cascades or double arrays obtained in Sandia National Laboratories [H*99]. The goal of numerical research in pulsed-power is to study the evolution of very intensive transient electric discharges and to perform a multiparametric optimization of future experimental schemes.

The studies of plasma shells (liners) magnetic implosion by multi-megaamper electric pulses and initiation of superintensive Z-pinches with the characteristic time of the order of 100 ns is a part of controlled fusion research. Plasma liner as well as classical pinch and plasma focus may be used as an intensive neutron source with 10⁹ to 3×10¹² neutrons pulse. The kinetic energy of the plasma shell may be used for compression and heating of thermonuclear target. The future prospect is application of the fast pinches technique for fusion ignition and creation of a neutron source with 10¹⁶ to 10¹⁸ neutrons per pulse. But the most wide application of plasma liners is for powerful soft X-rays sources (radiation with effective temperature 2-3 million Kelvin). In typical experimental schemes the energy is stored during the implosion (explosive compression) of plasma accelerated by the magnetic field to 300-500 km/sec. Then it is converted to thermal and subsequently to X-rays emission energy as a result of magnetic compression of the liner or collision of two liners. The characteristic energy range of X-rays sources based on plasma liners is from one to ten million Kelvin with the power from 10¹⁰ to 10¹⁴ W/cm. It allows planning investigations in the fields of X-ray laser, powerful X-rays emission impact on materials (laboratory research) X-ray lithography, materials hardening by radiation treatment, development of new diagnostics and tools in nanotechnology.

In practice numerical simulation is an essential part of large-scale experiments of this sort. Computer simulation in high-temperature plasma research projects saves both time and money. Numerical experiments cover considerable part of laborious work in planning natural experiments, prediction of experimental schemes efficiency, and experimental data analysis. In order to solve these problems the specialized scientific codes are created. One of them is MARPLE (Magnetically Accelerated Radiative PLasma Explorer), which is developed in IMM RAS. The code is based on the single-fluid two-temperature magnetohydrodynamic model [B63] taking into account radiative energy transfer [C51], supplemented with the prolonged plasma ablation model [AG*02], the electrical equation for the full circuit (generator, leading-in systems and the discharge chamber with the plasma in it) and computed/experimental data base (equations of state, transport and kinetic coefficients, opacity and emissivity coefficients) in the form of tables and analytical dependencies [NNU00]. Different types of 2D symmetry and full 3D problem statement are provided.

Up-to-date numerical methods and software engineering were applied for implementation of the above models and computer simulations. Complex geometry of experimental and industrial facilities as well as spatially non-uniformly scaled physical processes is handled by use of unstructured meshes: triangular/ quadrilateral/ blocked. The MHD system is solved by the generalized TVD Lax-Friedrichs scheme which was developed for the unstructured mesh applications. For the solution of parabolic equations describing the magnetic diffusion and conductive heat transfer, we developed the new finite-volume schemes constructed by analogy with mixed finite-element method. Diverse physical models processing and modification is supported by overall splitting scheme when the physical processes are included consequently.

The inhomogeneity of magnetically-driven plasmas is a known challenge for one who intends to calculate energy balance taking into account radiative processes. To overcome this problem we solve the radiative transport equation by means of semi-analytical characteristic algorithm. The analytical solution along the characteristic direction is constructed by means of the backward-forward angular approximation to the photon distribution function known as the Schwarzschield-Schuster technique. The two-group angular splitting gives an analytical expression for radiation intensity dependent on opacity and emissivity coefficients. The energy exchange between radiation field and plasma is taken into account via a radiative flux divergence, which is incorporated into the electron plasma component energy balance as some source function. A special class was developed for handling data tables including parameters of ionization, electric conductivity, thermodynamic and optical properties.

Flexible software architecture based on C++ programming provides updating possibility and readability of the code. The developed application software is compatible with different operating systems (MS Windows, Linux, UNIX) and different platforms (PC, powerful workstation, mainframe including parallel processing).

In order to ensure team work of all the project developers at all the stages - programming, simulations, further accompaniment of the code and its further development - the code construction is based on bloc principle. The program kernel includes:

• structured data base, which must contain all the information concerning the problem statement (simulation domain description, properties of media and materials, the "keys" of algorithms to be used) and the current state of the process (the computed physical values and auxiliary functions).
• a set of functional programs interacting with this data base (including preprocessor, computational modules and probably some tools for calculations control and primary results analysis).

These both parts are open, i.e. permitting extensions and revisions, including that by an expert user (with necessary limitations). Interaction between these structures is organized by supervisor (control module) receiving information from control files.

The object-oriented programming based on C++ language allows:

1. To develop unified descriptions of the program package components (data structures, their interrelations and handling procedures) thus facilitating and regulating the work of both developers and users.
2. To provide hierarchical access to the above components and specific access procedures at each level.
3. To organize cooperative implementation of the project including replacement/adding of personnel if necessary.
4. To organize step-by-step and combined verification, testing and adequacy estimation of the developed software.
5. To provide the possibilities of further development, enhancement and modification of both the separated elements and completed modules of the program package.
6. To provide software portability for different operating systems (MS Windows, Linux, Unix) and different platforms (PC, powerful workstation, mainframe including those with parallel processing).

 

6.2     MARPLE data preparation tools

MARPLE data preparation tools are developed according to the described above technique based on IC models.

3D MARPLE preprocessor is under contracture now. To illustrate the IC based modeling potential we represent an example of heterogeneous computational domain shown in Fig.5. This model initially was described in work [KAC*08].

               

                                     a)                                                                          b)

   

c)                                                                     d)

Fig. 5. The multicomponent mechanical assembly consisting of a gear, a shaft, and 2 end-round key: a) general view; b) shaft  c) , d) sectiones.

2D MARPLE preprocessor is fully completed. It is actively used for solving practical problems. 2D MARPLE preprocessor comprises all activities concerning data preparation, namely geometry definition, mesh generation, problem description (boundary conditions, material properties, and so on).

The preprocessor deals with several IC based models. Initial computational domain is described by an object of the ICModel consisting of cells of polyhedral and boundary types (P-cell, B-cell). The result of the computational domain discretization is represented by an object of the ICModel consisting of composite C-cells. The shape of C-cells is described by objects of the Grid type. Each of the composite cells represents subdivision of the corresponding cell of the initial ICModel. Then some of the neighboring cells are merged to form the resulting computational mesh that is described by ICModel consisting of C-cells or by GridModel, depending on homogeneity of the considered physical model.

We assume that the domain is a two-dimensional area defined by their boundary contours. It can consist of several components and may have holes. So the geometrical model can contain a few subdomains. Each 2D cell describes a subdomain which is unambiguously defined by their boundary contours that should be oriented. We use counterclockwise orientation for external contours and clockwise orientation for internal ones. Contours consist of segments. Cells of the initial IC model are formed in terms of the problem statement taking into account materials, physical properties distributions, requirements of applied numerical methods. For description of non-geometric attributes such as boundary conditions, material properties we use markers. Markers can be assigned to segments and subdomains of the initial geometric model as well as to elements, edges and nodes of the mesh.

MARPLE preprocessor allows to introducing periodical boundary conditions. Such conditions are interpreted as a consequence of the problem formulation symmetry. So definition of the periodical boundary conditions requires description of the corresponding plane transformations with listing all boundary elements that are identical in relation to these transformations.

MARPLE meshing tools provide automatic generation of triangle meshes according to the given elements size distribution with following optimization improving the mesh quality.

We assume that mesh density is described by a real-valued function S(X) which defined everywhere in the space E² and the value S(X_j) is interpreted as the desired mesh element size at the point X_j. The appropriate element size distribution may be prescribed by the user if he has knowledge of the physical situation a priori. For these purposes we use functional description based on so-called sources . Given a set of m sources, we calculate the overall element size distribution (mesh density) S(X) depending on all m sources Si(X) , as follows

 

 

Each source Si(X) describes the element size for a location X in the domain as a function of the distance r_i(X) to the source. In many cases the element size should be constant near the source (r_i<r_imin) and far from it (r_i>r_imax) and have to change smoothly between r_imin and r_imax . So we use the following function Si(r(X)) which provides the geometrical progression law of the element size increase/decrease.

Here h_imin, and h_imax are the limit values of the element size near the source and far from it correspondingly. And k_i is the coefficient of the progression . If k_i>1 then h_imin should be less than h_imax otherwise h_imin > h_imax.

Calculation of the distance r_i(X) from a point X to the i^th source depends on the shape of the source. We use point and line sources. The point source is defined by coordinates of its center. The line source is described by coordinates of the ends of the corresponding straight segment.

Mesh generation is illustrated by the Fig.6. We use the advancing front technique which allows to subdivide multi-component domains with holes according to the predefined mesh density function. Figure 6a illustrates the example of mesh density function defined by five sources (one of the point and four of the segment types). Colors are used for visualizing the distribution of the element size values. The corresponding adaptive mesh is shown in the Fig.6b.

 

a                                                         b

Fig. 6. Adaptive mesh  generation. a)  mesh density distribution b) adaptive mesh created by the advancing front method

6.3    Examples of simulations.

Since its origin, MARPLE was applied to many problems in the field of plasma dynamics, especially to those related with modern pulsed-power systems. Let's consider some of them.

6.3.1     Nozzle modeling.

Fig.7-9  illustrate data preparation process for modeling of a nozzle used in experiments on ultrafast magnetic implosion of plasma liners formed by gas-puff jets.An initial CAD model of an experimental installation is shown in Fig.7. One can see the annular nozzle-like flow path which surrounds central nozzle. Both nozzles, external and internal, are of Laval (converging-diverging) type. The corresponding axially-symmetric domain designed by MARPLE preprocessing tools for numerical calculations is represented in Fig.8. An adaptive mesh generated by MARPLE preprocessor is shown in Fig.9

Fig. 7. Initial CAD model of  an experimental installation.

Fig. 8. Computational domain designed by MARPLE preprocessor tools

Fig. 9. Fragments of adapted computational mesh

Another example of nozzle modeling is illustrated with Fig.10-11. Here Fig.10 represents a computational domain created for numerical simulation of a gas flow in a double ring nozzle, which was used in plasma physics experiments in attempt to obtain a double gas liner.

Fig. 10. Computational domain

Fig. 11 Steady density distribution patterns resulted from the flow simulation for the nozzle shown in Fig.10

6.3.2     Modelling of the electrodynamical implosion of the conical plasma liners

MHD numerical simulations were applied for multiparameter studies of 3D effects in magnetic compression of plasma liners created by conical multiwire arrays electrical explosion at «ANGARA-5-1» facility (TRINITI) with the discharge current 2 to 3  ÌÀ and the pulse rise time about 100 ns [GGA*08].

The description of the plasma dynamics at different stages of implosion was reproduced as a result of simulation. numerical and experimental time profiles of voltage drop at the load and soft X-ray yield power were compared. The effect of the geometry changes upon the implosion process was studied.

Numerical simulation is based on 2D RMHD code MARPLE using unstructured triangular grids. The model of prolonged plasma ablation [AG*02] was introduced to simulate plasma source. The governing MHD system was completed by electrical equation for the full circuit including the generator itself, leading-in systems and the discharge chamber with the plasma in it. Equations of state, transport and kinetic coefficients, opacity and emissivity coefficients are taken from the data tables [NNU00].

The behavior of the discharge is satisfactorily described in general by the above RMHD model. The plasma ablation model appeared to have a significant effect on both the entire scheme of plasma dynamics and such values as voltage drop at the load and soft X-ray yield power. The MARPLE code calibrated against the conical liners simulations proved to be a useful tool for computations aimed to optimization of the experimental setup for 3D implosion of plasma.

The results of the conical plasma liner implosion simulation for two different types of the discharge chamber are illustrated with Fig.12 and Fig.13. Contour plot of plasma density shown in Fig.12  allows to observe the bottleneck of the plasma flow at the bottom of the Z-axis. This result agrees with the experiment that was performed at ANGARA-5-1 (TRINITY, Russia ). Numerical calculations of the conical plasma liner implosion for the variant of the discharge chamber with specially fitted conical cathode insertion (bottom of the Z axis) were implemented for confirming a hypothesis about prevention of plasma flow bottleneck formation. Contour plot of plasma density for this type of the discharge chamber shown in Fig.13 demonstrates that the cathode insertion successfully prevents formation of the bottleneck at the bottom of the Z-axis. This MARPLE prediction is in good agreement with the corresponding experiment at ANGARA-5-1 (TRINITY, Russia ). It is a good demonstration of the predictive power of our numerical code.

Fig. 12. Conical plasma liner implosion simulation using the MARPLE code. Variant of the discharge chamber without conical cathode insertion (bottom of the Z axis). Plasma density patterns are shown at the moment close to that of maximum current. Strong peculiarity of the plasma flow (the bottleneck) is observed at the bottom of the axis. This result is in good agreement with the experiment that was performed at ANGARA-5-1 (TRINITY, Russia).

Fig. 13. Conical plasma liner implosion simulation using the MARPLE code. Plasma density patterns representing the result of the conical plasma liner implosion simulation for the variant of the discharge chamber with specially fitted conical cathode insertion (bottom of the Z axis). The cathode insertion successfully prevented formation of the bottleneck (see Fig. 12) at the bottom of the axis. This MARPLE prediction is in good agreement with the corresponding experiment at ANGARA-5-1 (TRINITY, Russia).

6.3.3     Magnetic flux compression by a plazma liner.

Magnetic Flux Compression (MFC) represents one of the modern promising schemes realizing magnetic energy storage with further energy compression in space and time. In particular, this scheme implies that an accelerated conductor can compress the flux of initially seeded azimuthal magnetic field. The accelerating force can be the Lorenz force of current created by a pulse-power generator. The generator produces the primary magnetic field accelerating a plasma shell (plasma MFC) which is supposed to compress and amplify the inner (secondary) magnetic field associated with the load current. This technique implies initiation of a hollow plasma cylinder from gas puffs, dielectric films or conducting foils. Therefore, in difference with the previous extensive MFC application at magneto-explosive generators the maximum shell velocity is determined by the Alfven velocity rather than by thermal velocity of explosion products. This opens possibilities for shortening of the load current pulse rise-time to ~100 ns and below.

Numerical, rather than experimental evaluation of the scheme feasibility becomes of particular importance when the experiments require energies of many mega joules. Numerical modeling could thus reveal physical limitations on the scheme operation efficiency and could play a determining role in engineering development of future thermonuclear facilities [GC03].

Despite limited experimental information available for the plasma MFC scheme, its theoretical efficiency was already investigated analytically via ideal MHD model and in simplified one-dimensional (1D) numerical modeling taking into account diffusive flux losses. One of the identified theoretical limitations on the scheme efficiency is related to diffusion of the compressed magnetic flux into the plasma of accelerated shell. In this case, a part of the flux remains frozen in the plasma during the short time of the load current increase and, therefore, can penetrate the load only with the plasma. The results represented in Fig14, Fig.15 illustrates our study of this limitation in two-dimensional (2D) geometry close to the experimental one. Magnetic flux compression by a plasma liner was simulated using the MARPLE code. The late stage of flux compression is illustrated with contour plot of plasma density shown in Fig.14.The Rayleigh-Taylor instability at both outer and inner surfaces of the plasma shell (the compressor) is well-developed. The instability considerably reduces efficiency of the magnetic flux compression scheme. Fig.15 represents contour plot of plasma density for the latest stage of flux compression. The Rayleigh-Taylor instability at the outer surface of the plasma shell is very strong. The instability at the inner surface of the shell is suppressed by this moment of time. The instability considerably reduces efficiency of the magnetic flux compression scheme. Observing Figs. 14 and 15 one can conclude that the optimal value of the initially seeded current should be taken close to the first of the considered variants, i.e. somewhat higher than 0.3 MA.

Fig. 14. Simulation using the MARPLE code: magnetic flux compression by a plasma liner at 1 MA initially seeded inner current. Density patterns are shown at the moment of maximal compression. The Rayleigh-Taylor instability at both outer and inner surfaces of the plasma shell (the compressor) is well-developed. The plasma of the liner-compressor does not penetrate into the load space due to the deceleration action of the compressed magnetic flux. However the flux compression is not quite effective because hydrodynamic instability strongly distorts the inner liner surface.

Fig. 15. Simulation using the MARPLE code: magnetic flux compression by a plasma liner at 0.3 MA initially seeded inner current. Density patterns are shown at the moment of maximal compression. c The Rayleigh-Taylor instability at the outer surface of the plasma shell is very strong and therefore reduces the compression efficiency due to the intensive mixing between the plasma and outer magnetic flux. The instability at the inner surface of the shell is suppressed by this moment of time. For this case the flux compression is satisfactory, however the plasma nearly penetrates into the load space.

7.     Conclusion

In this paper we considered an approach to designing the FE data preparation tools based on hybrid geometrical models represented by implicit complexes (IC).

ICs provide the background for description of heterogeneous geometrical models combining various representations and describing objects with arbitrary distribution of non-geometrical properties. Particular cases of ICs, polyhedral and cellular complexes, are suitable for representation of structured and unstructured grids of polyhedral or curvilinear elements with attributes. Application of the IC framework to FE data structures design and implementation allows us to develop unified software tools that meet the requirements of all the FE data preparation problems. In this paper we described our approach to realization of FE data processing tools based on the IC concepts and considered the software architecture of a radiative magneto-hydrodynamics simulation code MARPLE exploited this approach. We also represented examples of applied problems solved using MARPLE tools.

References