CN102663153B - Finite element modeling method for heterotype honeycomb structure - Google Patents

Abstract

The invention provides a finite element modeling method for a heterotype honeycomb structure, comprising thirteen steps: 1. establishing a cross section geometry model of the heterotype honeycomb structure; 2. generating a central point of the Kth ring honeycomb structure; 3. screening out an intra-domain honeycomb of the Kth ring honeycomb structure; 4. carrying out a step 11 if there is nonexistence of the intra-domain honeycomb in the Kth ring, otherwise carrying out a step 5; 5. generating nodes of each intra-domain honeycomb in the honeycomb structure on the Kth ring; 6. dealing with the nodes of the intra-domain honeycomb; 7. carrying out a step 8 if number of the nodes in the intra-domain honeycomb exceeds 2, otherwise carrying out a step 9; 8. dividing linear elements of the intra-domain honeycomb; 9. returning to carry out the step 5. if the finite element model of the intra-domain honeycomb on the Kth ring honeycomb structure is not completely generated, otherwise carrying out a step 10;10. returning to carry out the step 2, if the finite element model of the intra-domain honeycomb of thehoneycomb structures on all rings is not completely generated, otherwise carrying out a step 11; 11. outputting shell element information; 12.outputting node information; and 13.outputting element set information of the honeycomb structure sticking-joint.

Description

A kind of finite element modeling method of special-shaped honeycomb

(1) technical field

The invention provides a kind of finite element modeling method of special-shaped honeycomb, belong to aero propulsion technical field.

(2) background technology

Honeycomb comes from bionics, and last century five, the sixties, high performance honeycomb manufacturing technology is ripe, and honeycomb is applied to aerospace industry at large.Compare with traditional packing material, cellular structural material has the advantages such as large, the good energy absorption characteristics of low-density, lower thermal conductivity, high specific stiffness, high specific strength, compressive strain, high designability, is a kind of important compound substance that is applicable to aerospace field.

In general, traditional honeycomb mechanical characteristic analysis method has two kinds: method of equal effects and experimental method.Method of equal effects adopts even structureization theoretical, can derive the equivalent mechanical parameter of honeycomb, and its shortcoming is to concentrate the complex stress strain and the honeycomb cut-boundary that cause to be difficult to analyze on the impact of total mechanical characteristic for stress; And experimental method can only obtain the empirical value of honeycomb mechanical characteristic, be difficult to the design criteria of formation and development honeycomb.

In current engineering, applying more is the special-shaped honeycomb of excision forming, its structure stress is complicated, and the mechanical characteristic of boundary has very important impact to integrally-built intensity, rigidity and vibration characteristics, obviously, traditional method of equal effects and experimental method are difficult to meet the analysis needs of its structural mechanics characteristic, therefore should tend to use finite element method to the mechanical characteristic analysis of special-shaped honeycomb.Yet the difficult point of using finite element method to carry out mechanical characteristic analysis to honeycomb is that honeycomb is difficult to set up finite element model, use existing finite element pre-processing software to carry out finite element modeling to honeycomb and can only adopt manual mode, this will expend a lot of time and resource, and for some large scale structures, manual modeling or even impossible realization.

Therefore, need to develop a kind of finite element modeling method of special-shaped honeycomb, realize the quick division of special-shaped honeycomb finite element grid, to meet requirement of engineering.

(3) summary of the invention

The finite element modeling method that the object of this invention is to provide a kind of special-shaped honeycomb, it has solved the deficiencies in the prior art.

The present invention is intended to special-shaped honeycomb to divide finite element grid, and this honeycomb is the hexagon honeycomb the most often using in engineering, and structure boundary adopts excision forming, and cut place exists multiple imperfect honeycomb; Between honeycomb and honeycomb, by bonding way, be connected, have a bonding direction, and in bonding direction, have the bonding thickness between honeycomb; For obtaining different structural mechanics characteristics, honeycomb can be arranged along any direction.

The fundamental design idea of a kind of special-shaped honeycomb finite element modeling method of the present invention is: first, extract a cross section of special-shaped honeycomb, and according to the geometrical property of hexagon honeycomb, on this cross section, carry out arranging of node, its way is on honeycomb cross section, to get a honeycomb central point as reference point, by reference point, to outer ring-like, generated the node of the honeycomb in this cross section, and carry out special processing to being positioned at the node of the imperfect honeycomb being truncated in structure boundary; Then, the node of take on cross section is basis, utilizes the shell unit of the method generating three-dimensional honeycomb of thickness direction interpolation; Finally, special-shaped honeycomb finite element model information is output as to the file that can be read by common finite element process software.

The finite element modeling method of a kind of special-shaped honeycomb of the present invention, its concrete steps are as follows:

Step 1: the cross section geometric model of setting up special-shaped honeycomb;

First, according to the border characteristic of given special-shaped honeycomb, select a cross section perpendicular to honeycomb thickness direction, set up the cross section geometric model of this honeycomb, in the present invention, the cross section geometric model of honeycomb is based upon on the x-y face of cartesian coordinate system, and the thickness direction of honeycomb is positioned at z direction of principal axis;

Then, cross section geometric model partition is become to the quadrilateral subregion of a plurality of regular shape, in the present invention, need repeatedly put the judgement that whether is positioned at plane geometry model area, and the interpolation calculation of honeycomb thickness direction, these two kinds of calculating are all based on quadrangular plan geometric areas, therefore need suitably divide according to model actual conditions pair cross-section geometric model;

Finally, because the present invention carries out finite element modeling according to the geometrical property of hexagon honeycomb, the bonding direction of method hypothesis honeycomb is positioned at y direction of principal axis, therefore if desired the bonding direction of honeycomb is positioned at other direction, the cross section geometric model of honeycomb to be rotated to respective angles in coordinate system and carry out finite element modeling, after modeling completes, finite element model be rotated back to former angle, wherein, with 1 P ₀(x ₀, y ₀) be rotary middle point, will put P (x, y) and be rotated counterclockwise θ angle and obtain a P _t(x _t, y _t) the computing of coordinate figure be expressed in matrix as:

$\left[\begin{array}{cc}{x}_t& y_t\end{array}\right]=\left[\begin{array}{cc}{x}_0& y_0\end{array}\right]+\left(\begin{array}{c}\left[\begin{array}{cc}x& y\end{array}\right]\end{array}-\left[\begin{array}{cc}{x}_0& y_0\end{array}\right]\right)\left[\begin{array}{cc}\cos \mathrm{\θ}& \sin \mathrm{\θ}\\ -\sin \mathrm{\θ}& \cos \mathrm{\θ}\end{array}\right]---\left(1\right)$

Step 2: the honeycomb central point that generates K (0≤K≤Kp) ring honeycomb;

Select a honeycomb central point C (x on cross section geometric model _c, y _c) as reference point, outwards generate the honeycomb central point of K (0≤K≤Kp) ring honeycomb, centered by C point, the honeycomb of point is considered as the 0th ring here, is outwards followed successively by the 1st, 2 ..., Kp ring, Kp is the honeycomb number of rings upper limit default in the present invention;

When K=0, the honeycomb central point on K ring is C (x _c, y _c), when K ≠ 0, generate as follows K ring honeycomb central point:

First, on generation K (1≤K≤Kp) ring, honeycomb central point and C are ordered line and x axle clamp angle are the honeycomb central point of (i=1,2..., 6), so every ring of honeycomb central point has six, and its coordinate is:

$\left\{\begin{array}{c}{x}_{k,i}=\sqrt{3}\mathrm{Kh}\cos [\frac{\mathrm{\π}}{6}+(i-1)\frac{\mathrm{\π}}{3}]+x_C\\ y_{k,i}=\sqrt{3}\mathrm{Kh}\sin [\frac{\mathrm{\π}}{6}+(i-1)\frac{\mathrm{\π}}{3}]+y_C\end{array}\right.i=\mathrm{1,2},..,6---\left(2\right)$

Wherein, the length of side that h is honeycomb;

Then, consider that honeycomb exists bonding thickness on y direction of principal axis, by the above-mentioned central point generating with the y direction of principal axis between reference point C apart from expanding Ap times, i.e. y ' _{k, i}=y _c+ Ap (y _{k, i}-y _c), wherein, Ap is the honeycomb width ratio of the bonding thickness of consideration, t is the thickness of two honeycomb abutting edges, and the actual (real) thickness of two honeycomb abutting edges and honeycomb cell-wall thickness is poor;

Finally, generate remaining honeycomb central point on K (2≤K≤Kp) ring, these honeycomb central points are positioned on the Along ent of every pair of adjacent line of generating center point, choose the two adjacent central point C that generated _{k, i}(x _{k, i}, y ' _{k, i}) and C _{k, i+1}(x _{k, i+1}, y ' _{k, i+1}), by these 2 the honeycomb central point C that can generate between the two _{k, j}coordinate figure be:

$\left\{\begin{array}{c}{x}_{k,j}=x_{k,i}(1-q)+x_{k,i+1}q\\ y_{k,j}={y^{\′}}_{k,i}(1-q)+{y^{\′}}_{k,i+1}q\end{array}\right.j=1,...,K-1---\left(3\right)$

Wherein, q=j/K, such operation is carried out six times altogether, remaining honeycomb central point on K ring can be generated;

Step 3: filter out the honeycomb (hereinafter to be referred as territory honeycomb) that angle point is positioned at cross section geometric model area in K (0≤K≤Kp) ring honeycomb, in note K ring territory, the quantity of honeycomb is nf;

First, each honeycomb central point C ' (x to K (0≤K≤Kp) ring _i, y _i), write out its corresponding honeycomb j (j=1,2 ..., 6) coordinate of individual angle point:

$\left\{\begin{array}{c}{x}_{i,j}=x_i+h\cos \frac{(j-1)\mathrm{\π}}{3}\\ y_{i,j}=y_i+h\sin \frac{(j-1)\mathrm{\π}}{3}\end{array}\right.j=\mathrm{1,2},...,6---\left(4\right)$

Wherein, the length of side that h is honeycomb;

Then, to the j of this honeycomb (j=1,2 ..., 6) individual angle point (x _{i, j}, y _{i, j}), judge whether it is positioned at cross section geometric model area, if this honeycomb has angle point to be positioned at cross section geometric model area, the center point coordinate of this honeycomb is recorded in a matrix cmf, finally obtain honeycomb number, i.e. the line number nf of matrix cmf in the territory of K ring;

Judge 1 P (x _p, y _p) whether be positioned at geometric model region and can use numerical inverse isoparametric mapping interpolation or triangle area method, two kinds of applicable geometric models of method are all plane quadrilateral regions, therefore be a plurality of quadrilateral areas by the cross section geometric model partition of special-shaped honeycomb in step 1 of the present invention, for these a plurality of quadrilateral areas, put respectively the judgement that whether is positioned at its territory, if exist such region to make a little within it in the model of cross section, explanation point is positioned at cross section geometric model area;

Wherein, numerical inverse isoparametric mapping interpolation is according to the analyticity confrontation quadrilateral area of quadrilateral area and puts P (x _p, y _p) carry out coordinate transform, after coordinate transform, quadrilateral area becomes D{ (ξ, η) |-1≤ξ≤1 ,-1≤η≤1}, and the coordinate of point becomes P ' (ξ _p, η _p), if some P is (ξ _p, η _p) be located at region D{ (ξ, η) |-1≤ξ≤1, within-1≤η≤1}, an explanation point P (x _p, y _p) be positioned at former quadrilateral area; Triangle area method is with a P (x _p, y _p) be four triangles of four edges structure of summit and quadrilateral area Ω to be judged, if the area that four leg-of-mutton area sums are less than or equal to quadrilateral area Ω illustrates a P (x _p, y _p) be positioned at region Ω;

Step 4: if nf < 1, K (0≤K≤Kp) encircles the not interior honeycomb of existence domain, illustrate that the honeycomb central point having generated has covered the cross section geometric model area of special-shaped honeycomb completely, and then execution step 11, otherwise continue execution step five, the effect of this step is the waste of avoiding computational resource;

Step 5: generate K (0≤K≤Kp) ring honeycomb i (i=1 ..., the nf) node of honeycomb on cross section geometric model in individual territory, and remember that the node number that this honeycomb is positioned at geometric model region is nin1;

This step generates the i (i=1 of K (0≤K≤Kp) ring honeycomb, ..., nf) node of honeycomb on cross section geometric model in individual territory, and each nodal information of honeycomb in this territory is write to a matrix N Inf, NInf is the matrix of 6nd * 5, nd is the cell density on honeycomb limit, and the information that the first row to the of NInf five row will record is respectively: whether angle point, node that whether node ID, node are positioned at honeycomb are positioned at the horizontal ordinate of geometric model region, node, the ordinate of node;

Generate by the following method the node of this honeycomb:

First, for center point coordinate, be (x _i, y _i) honeycomb, each node coordinate value of its correspondence is:

$\left\{\begin{array}{c}{x}_{i,j}=\frac{\sin \frac{\mathrm{\&pi;}}{3}h}{\sin (\frac{2\mathrm{\&pi;}}{3}-\mathrm{\&alpha;})}\cos \mathrm{\&theta;}+x_i\\ y_{i,j}=\frac{\sin \frac{\mathrm{\&pi;}}{3}h}{\sin (\frac{2\mathrm{\&pi;}}{3}-\mathrm{\&alpha;})}\sin \mathrm{\&theta;}+y_i\end{array}\right.j=\mathrm{1,2},...,6\mathrm{nd}---\left(5\right)$

Wherein, the length of side that h is honeycomb, (j=1,2 ..., 6nd),

And, because honeycomb exists bonding thickness on y direction of principal axis, therefore need be by each node and central point (x _i, y _i) distance on y direction of principal axis expands as original Ap doubly, i.e. y ' _{i, j}=y _i+ (y _{i, j}-y _i) Ap, wherein, Ap is the honeycomb width ratio of the bonding thickness of consideration, h is the length of side of honeycomb, and t is the bonding thickness of two honeycomb abutting edges;

Then, to the information that writes each node in matrix N Inf:

The coordinate figure of the node of this honeycomb is write to the 4th row and the 5th row of matrix N Inf, the node ID of this honeycomb (above-mentioned j value) is write to the first row of matrix N Inf, position attribution by the node of this honeycomb in honeycomb writes the secondary series of matrix N Inf, if node be positioned at honeycomb angle point element value be 1, if it is 0 that node is positioned on honeycomb limit element value, can adopt formula m=mod (j-1, nd) carry out the judgement of the position attribution of node in honeycomb, position attribution by the node of this honeycomb in geometric model writes the 3rd row of matrix N Inf, if node be positioned at geometric model region element value be 1, if it is 0 that node is positioned at outside geometric model region element value,

I (i=1 ..., nf), after the node of honeycomb all generates in individual territory, find out the node number that is positioned at geometric model region in this honeycomb, the line number that in matrix N Inf, the 3rd column element value is 1, is denoted as nin1;

Step 6: to the i of K (0≤K≤Kp) ring honeycomb (i=1 ..., nf), honeycomb is processed in individual the territory in, and remembers that the node number that this honeycomb is after treatment positioned at geometric model region is nin2;

This step is intended to processing and is positioned at the incomplete honeycomb on cross section geometric model boundary, it is taked to special node arrangement method, to reach mating of finite element model and geometric model, i (i=1 to K (0≤K≤Kp) ring honeycomb, ..., nf) after in individual territory, honeycomb is processed, the information matrix NInf of this honeycomb has corresponding renewal, the method that the node of honeycomb is processed will after describe in detail;

After processing, i (i=1 ..., nf) in individual territory will there is certain variation in the nodal community of honeycomb, find out the node number that this honeycomb is now positioned at geometric model region, the line number that in matrix N Inf, the 3rd column element value is 1, is denoted as nin2;

Step 7: judge whether nin2 is greater than 2, if be greater than 2, continue execution step eight, otherwise directly perform step nine;

Nin2≤2 explanation i (i=1, ..., nf) in individual territory, honeycomb can be regarded as approx and only has a honeycomb limit to be positioned on cross section geometric model boundary, in actual conditions, this honeycomb limit does not exist, therefore abandon this honeycomb, and then honeycomb in i+1 territory of processing K (0≤K≤Kp) ring;

Step 8: the i of division K (0≤K≤Kp) ring honeycomb (i=1 ..., nf) the line unit of honeycomb in individual territory;

The territory interior nodes of honeycomb is write to matrix nodes in order successively, from matrix nodes, choose suitable node and construct the line unit of honeycomb, the node number that forms line unit is write to a matrix elems, for honeycomb in complete territory, in its honeycomb, every two adjacent nodes generate a line unit, for the incomplete honeycomb being positioned on cross section geometric model boundary, every two adjacent territory interior nodes generate a line unit, but do not generate line unit between first territory interior nodes and last territory interior nodes, be that honeycomb is not had cell wall structures by model boundary cut place,

Because honeycomb abutting edge exists bonding thickness, while carrying out finite element analysis, abutting edge shell unit need to be assigned to other one-tenth-value thickness 1/10, therefore the shell unit of abutting edge need be classified as to a set, with matrix sets, record these unit, this modeling method hypothesis y direction of principal axis is bonding direction, will in each honeycomb, perpendicular to the line unit number on the honeycomb limit of y axle, write matrix sets, concerning complete honeycomb, the order generating by its node, the condition that in honeycomb, n unit is unit, abutting edge is: by (n-1)/nd, to zero value of rounding, be 1 or 4, wherein, nd is the cell density on honeycomb limit,

After line unit on cross section geometric model all generates, matrix nodes stores the whole nodal informations on cross section geometric model, matrix elems stores the whole line unit informations on cross section geometric model, and matrix sets stores the unit number of honeycomb abutting edge on cross section geometric model;

Step 9: i=i+1, and judge whether i value is greater than nf, if be less than or equal to nf, return to execution step five, if be greater than nf, perform step ten;

I > nf illustrates that the finite element model of honeycomb in the territory of K (0≤K≤Kp) ring set up, and then processes K+1 and encircle honeycomb;

Step 10: K=K+1, and judge whether K value is greater than Kp, if be less than or equal to Kp, return to execution step two, if be greater than Kp, illustrate that the finite element model of the maximum number of rings honeycomb of this method regulation all generates, continue execution step 11;

Step 11: output shell unit information;

Obtain cross section geometric model reach the standard grade unit sum ElemNum and node sum NodeNum, suppose has NLay layer shell unit on z direction of principal axis;

For iLay (iLay=1 ..., NLay) layer on j (j=1, ..., ElemNum) individual shell unit, its node number is n1, n2, n1+NodeNum and n2+NodeNum, obviously, during iLay=1, n1 and n2 are the corresponding line cell node number in matrix elems, during iLay > 1, n1 and n2 i.e. the node number of two higher values of iLay-1 layer shell unit;

According to output requirement, the node number that unit is comprised is by row writing unit information output file;

Step 12: output node information;

For the node being positioned on cross section geometric model, its ordinate is 0, the ordinate of other node utilizes the interpolation method of thickness direction to obtain according to special-shaped honeycomb geometric model characteristic, it for ordinate value, is not 0 node, for node coordinate vector coord3 of its structure, the first two element of coord3 is horizontal ordinate and the ordinate value that this node correspondence records in matrix nodes, the 3rd element of coord3, the i.e. ordinate value of this node, the method that four angular coordinates of employing based on geometric model quadrilateral subregion carry out planar interpolation obtains, the final finite element model forming will be the plane of burst at the end face of honeycomb, or more accurate, the quadrilateral subregion at the projection place based on node on x-y face selects corresponding interpolating function to carry out surface interpolation, end face at honeycomb is obtained and the realistic model patch surface of exact matching more,

Coordinate transform by the coordinate of each node through over-angle rotation, obtain the node coordinate under former angle position, by output requirement output, here should be noted that the output order of nodal information, the node being positioned on cross section geometric model should be according to the order of the storage of matrix nodes, other every node layer all should, corresponding to the storage order of matrix nodes, so just can make nodal information and unit information defeated corresponding;

Step 13: the unit set information of output honeycomb abutting edge;

According to output requirement, the unit number of honeycomb abutting edge is write to aggregate information output file.

This self geometrical property according to hexagon honeycomb is carried out the method for finite element modeling, can generate fast the finite element model of a special-shaped honeycomb, and, the generation of finite element model does not depend on entity, do not need to set up the three-dimensional entity model of special-shaped honeycomb, also abandoned traditional solid model of take is basic manual modeling method simultaneously, has greatly improved modeling efficiency.The present invention has considered that honeycomb exists bonding thickness between honeycomb and the angle problem of arranging of honeycomb, and the imperfect honeycomb being positioned on model boundary is processed, and makes finite element model reach the effect matching with solid model.

Wherein, in the finite element modeling method step 6 of a kind of special-shaped honeycomb of the present invention, the method that the node of honeycomb is processed, its concrete implementation step is as follows:

Step 1: find out the node number that this honeycomb is positioned at cross section geometric model area, the line number that in the nodal information matrix N Inf of this honeycomb, the 3rd column element value is greater than zero, is denoted as nin;

Step 2: judge whether nin equals 6rd, if nin=6nd illustrates that this honeycomb is honeycomb in a complete territory, directly perform step 13, if nin ≠ 6nd illustrates that honeycomb is blocked by cross section geometric model boundary, lost integrality, execution step three;

Step 3: the node to this honeycomb sorts;

The object of this step is: make in the nodal information matrix N Inf of this honeycomb the node of each row representative from top to bottom in honeycomb, become order counterclockwise to arrange, and the first row representative in matrix N Inf is by the first point in counterclockwise tactic cross section geometric model area;

The specific implementation method of this step is: the first row of matrix N Inf is moved to last column, and all the other each row upwards move respectively a line, and this operation is carried out for several times until matrix N Inf meets the requirement of this step;

Step 4: in a row vector pend, obviously, pend represents will enter by counterclockwise order the nodal information of cross section geometric model area by the information storage of last column of matrix N Inf;

Step 5: stipulate an error range herr;

Step 6: process this honeycomb and be positioned at cross section geometric model area by counterclockwise tactic last Nodes;

First, last node in cross section geometric model area is denoted as to P ₁(nin of homography NInf is capable), P ₁next node be denoted as P ₂(nin+1 of homography NInf is capable), tries to achieve P ₁and P ₂the intersection point pin1 of line and cross section geometric model boundary;

Then, investigate P ₁and the distance between pin1, if distance is greater than herr, extra-regional node is moved on on cross section geometric model boundary, the coordinate figure that is about to some pin1 is assigned to the 4th capable row and the 5th row of nin+1 of matrix N Inf, simultaneously, the attribute of the node of the capable representative of nin+1 of matrix N Inf changes " not being positioned at honeycomb angle point " and " being positioned at cross section geometric model area " into, and the capable secondary series element value of nin+1 that is about to matrix N Inf changes 0, the three column element value into and changes 1 into;

If P ₁and the distance between pin1 is less than or equal to herr, consider the situation that nin is 1, this illustrates that this honeycomb only has a node in region, and some P in boundary ₁with the distance of intersection point pin1 within error range, this honeycomb is given up, turn back to the step 12 of this method main flow, if nin is not 1, last node in region is moved on on cross section geometric model boundary, the coordinate figure that is about to some pin1 is assigned to the 4th capable row and the 5th row of nin of matrix N Inf, simultaneously, the attribute of the node of the capable representative of nin of matrix N Inf changes " not being positioned at honeycomb angle point " and " being positioned at cross section geometric model area " into, the capable secondary series element value of nin that is about to matrix N Inf changes 0, the three column element value into and changes 1 into;

Step 7: find out the line number that the 3rd column element value of matrix N Inf is greater than 0, and be again assigned to nin, this represents that this honeycomb is through the operation of above step, and the node number nin that this honeycomb is positioned at cross section geometric model area upgrades;

Step 8: process this honeycomb and be positioned at cross section geometric model area by counterclockwise tactic first Nodes;

First, node row vector pend being recorded is denoted as P ₁, by cross section geometric model area first be denoted as P ₂(the 1st row of homography NInf), the intersection point of trying to achieve honeycomb and cross section geometric model boundary is pin2;

Then, investigate P ₂and the distance between pin2, if distance is greater than herr, consider that nin equals the situation of the line number of matrix N Inf, now this honeycomb all becomes point in territory through all nodes of operation of step 6, need to increase 6nd+1 node in last interpolation a line of matrix N Inf, the content of this line can temporarily copy the content of last column of matrix N Inf, if nin is not equal to the line number of matrix N Inf, do not need to do the operation that increases row;

Last column of matrix N Inf is moved to the first row, all the other each row move respectively a line backward, the first row of matrix N Inf is now represented to the information of first node in territory, the coordinate figure of pin2 is assigned to the 4th row and the 5th row of matrix N Inf the first row, simultaneously, the nodal community of matrix N Inf the first row representative changes " not being positioned at honeycomb angle point " and " being positioned at cross section geometric model area " into, and the secondary series element value that is about to the first row of matrix N Inf changes 0, the three column element value into and changes 1 into;

If P ₂and the distance between pin2 is less than or equal to herr, first node in region is moved on on cross section geometric model boundary, the coordinate figure that is about to pin2 is assigned to the 4th row and the 5th row of matrix N Inf the first row, simultaneously, the nodal community of matrix N Inf the first row representative changes " not being positioned at honeycomb angle point " and " being positioned at cross section geometric model area " into, the secondary series element value that is about to the first row of matrix N Inf changes 0, the three column element value into and changes 1 into;

Step 9: find out the line number that the 3rd column element value of matrix N Inf is greater than 0, and be again assigned to nin, this represents that this honeycomb is through the operation of above-mentioned steps, and the node number nin that this honeycomb is positioned at cross section geometric model area upgrades;

Step 10: be less than herr if honeycomb is positioned at first node and last nodal distance of cross section geometric model area, carry out node merging;

Investigate the represented node of the 1st row of matrix N Inf and the distance value between the capable represented node of nin, if distance is less than herr, the represented node of the 1st row of matrix N Inf is moved on to the midpoint of two nodes, soon the 4th of the first row of matrix N Inf the row and the 5th column element value change the mid point coordinate figure of two nodes into, simultaneously, the attribute of nin node is changed into " being positioned at cross section geometric model area outer ", be about to the capable tertial element value of nin and change 0 into;

Step 11: find out and be positioned at cross section geometric model area in this honeycomb and be positioned at the node number on honeycomb angle point, be that in matrix N Inf, secondary series and the 3rd column element value are all 1 line number, be denoted as ncd, if ncd < 1, illustrate that this honeycomb does not have angle point to be positioned at cross section geometric model area, returns to the step 12 in this method main flow, if ncd >=1, find out the line number that the 3rd column element value of matrix N Inf is greater than 0, and be again assigned to nin;

Step 12: judge the pattern of this honeycomb and carry out cellular node sequence;

First, remember that this honeycomb first node being arranged on honeycomb angle point in cross section geometric model area is positioned at pos at matrix N Inf ₁oK, and remember d ₁=pos ₁-1, remember that the node that last in cross section geometric model area of this honeycomb is arranged on honeycomb angle point is positioned at pos at matrix N Inf ₂oK, and remember d ₂=nin-pos ₂;

Then, if d ₁> d ₂, need the node of this honeycomb to become following order: the last point in former territory is the first point (the first row of homography NInf) in territory, and each node is by order arrangement clockwise;

Step 13: find out the line number that the 3rd column element value of matrix N Inf is greater than 0, and be again assigned to nin, this represents that this honeycomb is through the operation of above-mentioned steps, and the node number nin that this honeycomb is positioned at cross section geometric model area upgrades.

When honeycomb is processed in this method, need to ask the intersecting point coordinate of two line segments of 4 formations, concrete grammar is as follows:

First, suppose two line segments with end points coordinate be respectively P ₁(x ₁, y ₁), P ₂(x ₂, y ₂), P ₃(x ₃, y ₃), P ₄(x ₄, y ₄);

Then, structure following parameters:

$\left\{\begin{array}{c}{d}_1=-(x_3-x_4)(y_1-y_2)+(x_1-x_2)(y_3-y_4)\\ d_2=-(x_3-x_4)(y_1-y_2)+(x_1-x_2)(y_3-y_4)\end{array}\right.---\left(6\right)$

If d ₁d ₂=0, two line segments are parallel, without intersection point, otherwise, the intersection point of two line segments (or intersection point of extended line) P ₀the coordinate figure of (x, y) is:

$\left\{\begin{array}{c}x=\frac{x_1[x_4(-y_2+y_3)+x_3(y_2-y_4)]+x_2[x_4(y_1-y_3)+x_3(-y_1+y_4)]}{d_1}\\ y=\frac{x_4(y_1-y_2)y_3+x_1{y}_2(y_3-y_4)+x_3(-y_1+y_2)y_4+x_2{y}_1(-y_3+y_4)}{d_2}\end{array}\right.---\left(7\right)$

In above Overall Steps, the main symbol of using is as shown in table 1, and all symbols have all provided definition when occurring for the first time.

Table 1

The invention has the advantages that:

(1) can utilize the geometrical property of hexagon honeycomb self to generate fast the finite element grid of special-shaped honeycomb, do not rely on the solid model of honeycomb, the modeling process of having avoided the 3D solid of honeycomb, has improved modeling efficiency, is applicable to engineering reality;

(2) on the cut-boundary of special-shaped honeycomb, there are a large amount of imperfect honeycombs, and borderline imperfect honeycomb has very important impact to mechanical characteristics such as the intensity of total, rigidity, vibrations, if wish, special-shaped honeycomb is carried out to effective mechanics property analysis, inevitable requirement finite element model has higher accuracy, the present invention has realized the accurate foundation of finite element model on honeycomb border, significant to the mechanics property analysis of honeycomb;

(3) the present invention is applicable to have the honeycomb of complex boundary, needs only the boundary Control point to fixed structure, can generate fast the finite element model of hexagon honeycomb, has reached the effect of the intelligent border cutting of finite element model of honeycomb;

(4) finite element modeling method of a kind of special-shaped honeycomb provided by the invention, have efficient, generate fast the feature of honeycomb grid, for given honeycomb, only need to control several parameters, as the cell density on honeycomb limit, the shell unit number of plies on honeycomb thickness etc., just can generate fast satisfactory honeycomb finite element grid, same, in the design process of honeycomb, as long as the given honeycomb length of side, bonding thickness, the parameters such as bonding direction, get final product the finite element grid of the multiple honeycomb form of a member of Rapid Establishment, facilitated designing and calculating.

(4) accompanying drawing explanation

Fig. 1 is the process flow diagram of the finite element modeling method of a kind of special-shaped honeycomb of the present invention;

Fig. 2 carries out the process flow diagram of disposal route to the node of honeycomb in the present invention;

Fig. 3 is the geometric parameter of the selected example of the present invention;

The meaning of label: A in Fig. 3---the part of selected example geometric model, corresponding with Fig. 4;

Fig. 4 is the local geometric model of the selected example of the present invention and the position of reference point while setting up finite element model;

Fig. 5 is the bottom surface geometric model of honeycomb in example of the present invention;

Fig. 6 is that the line unit on the bottom surface geometric model of example of the present invention generates situation;

Fig. 7 is the honeycomb finite element model of example of the present invention generation situation in Hypermesh.

(5) embodiment

Below in conjunction with accompanying drawing and example, the present invention is described in further detail.

The finite element modeling method that the object of this invention is to provide a kind of special-shaped honeycomb, it has solved the deficiencies in the prior art.

The present invention is intended to special-shaped honeycomb to divide finite element grid, and this honeycomb is the hexagon honeycomb the most often using in engineering, and structure boundary adopts excision forming, and cut place exists multiple imperfect honeycomb; Between honeycomb and honeycomb, by bonding way, be connected, have a bonding direction, and in bonding direction, have the bonding thickness between honeycomb; For obtaining different structural mechanics characteristics, honeycomb can be arranged along any direction.

The fundamental design idea of a kind of special-shaped honeycomb finite element modeling method of the present invention is: first, extract a cross section of special-shaped honeycomb, and according to the geometrical property of hexagon honeycomb, on this cross section, carry out arranging of node, its way is on honeycomb cross section, to get a honeycomb central point as reference point, by reference point, to outer ring-like, generated the node of the honeycomb in this cross section, and carry out special processing to being positioned at the node of the imperfect honeycomb being truncated in structure boundary; Then, the node of take on cross section is basis, utilizes the shell unit of the method generating three-dimensional honeycomb of thickness direction interpolation; Finally, special-shaped honeycomb finite element model information is output as to the file that can be read by common finite element process software.

As shown in Figure 1, its concrete steps are as follows for the process flow diagram of the finite element modeling method of a kind of special-shaped honeycomb of the present invention:

Step 1: the cross section geometric model of setting up special-shaped honeycomb;

First, according to the border characteristic of given special-shaped honeycomb, select a cross section perpendicular to honeycomb thickness direction, set up the cross section geometric model of this honeycomb, in the present invention, the cross section geometric model of honeycomb is based upon on the x-y face of cartesian coordinate system, and the thickness direction of honeycomb is positioned at z direction of principal axis;

Then, cross section geometric model partition is become to the quadrilateral subregion of a plurality of regular shape, in the present invention, need repeatedly put the judgement that whether is positioned at plane geometry model area, and the interpolation calculation of honeycomb thickness direction, these two kinds of calculating are all based on quadrangular plan geometric areas, therefore need suitably divide according to model actual conditions pair cross-section geometric model;

Finally, because the present invention carries out finite element modeling according to the geometrical property of hexagon honeycomb, the bonding direction of method hypothesis honeycomb is positioned at y direction of principal axis, therefore if desired the bonding direction of honeycomb is positioned at other direction, the cross section geometric model of honeycomb to be rotated to respective angles in coordinate system and carry out finite element modeling, after modeling completes, finite element model be rotated back to former angle, wherein, with 1 P ₀(x ₀, y ₀) be rotary middle point, will put P (x, y) and be rotated counterclockwise θ angle and obtain a P _t(x _t, y _t) the computing of coordinate figure be expressed in matrix as:

$\left[\begin{array}{cc}{x}_t& y_t\end{array}\right]=\left[\begin{array}{cc}{x}_0& y_0\end{array}\right]+\left(\begin{array}{c}\left[\begin{array}{cc}x& y\end{array}\right]\end{array}-\left[\begin{array}{cc}{x}_0& y_0\end{array}\right]\right)\left[\begin{array}{cc}\cos \mathrm{\&theta;}& \sin \mathrm{\&theta;}\\ -\sin \mathrm{\&theta;}& \cos \mathrm{\&theta;}\end{array}\right]---\left(8\right)$

Step 2: the honeycomb central point that generates K (0≤K≤Kp) ring honeycomb;

Select a honeycomb central point C (x on cross section geometric model _c, y _c) as reference point, outwards generate the honeycomb central point of K (0≤K≤Kp) ring honeycomb, centered by C point, the honeycomb of point is considered as the 0th ring here, is outwards followed successively by the 1st, 2 ..., Kp ring, Kp is the honeycomb number of rings upper limit default in the present invention;

When K=0, the honeycomb central point on K ring is C (x _c, y _c), when K ≠ 0, generate as follows K ring honeycomb central point:

First, on generation K (1≤K≤Kp) ring, honeycomb central point and C are ordered line and x axle clamp angle are the honeycomb central point of (i=1,2..., 6), so every ring of honeycomb central point has six, and its coordinate is:

$\left\{\begin{array}{c}{x}_{k,i}=\sqrt{3}\mathrm{Kh}\cos [\frac{\mathrm{\&pi;}}{6}+(i-1)\frac{\mathrm{\&pi;}}{3}]+x_C\\ y_{k,i}=\sqrt{3}\mathrm{Kh}\sin [\frac{\mathrm{\&pi;}}{6}+(i-1)\frac{\mathrm{\&pi;}}{3}]+y_C\end{array}\right.i=\mathrm{1,2},..,6---\left(9\right)$

Wherein, the length of side that h is honeycomb;

Then, consider that honeycomb exists bonding thickness on y direction of principal axis, by the above-mentioned central point generating with the y direction of principal axis between reference point C apart from expanding Ap times, i.e. y ' _{k, i}=y _c+ Ap (y _{k, i}-y _c), wherein, Ap is the honeycomb width ratio of the bonding thickness of consideration, t is the thickness of two honeycomb abutting edges, and the actual (real) thickness of two honeycomb abutting edges and honeycomb cell-wall thickness is poor;

Finally, generate remaining honeycomb central point on K (2≤K≤Kp) ring, these honeycomb central points are positioned on the Along ent of every pair of adjacent line of generating center point, choose the two adjacent central point C that generated _{k, i}(x _{k, i}, y ' _{k, i}) and C _{k, i+1}(x _{k, i+1}, y ' _{k, i+1}), by these 2 the honeycomb central point C that can generate between the two _{k, j}coordinate figure be:

$\left\{\begin{array}{c}{x}_{k,j}=x_{k,i}(1-q)+x_{k,i+1}q\\ y_{k,j}={y^{\&prime;}}_{k,i}(1-q)+{y^{\&prime;}}_{k,i+1}q\end{array}\right.j=1,...,K-1---\left(10\right)$

Wherein, q=j/K, such operation is carried out six times altogether, remaining honeycomb central point on K ring can be generated;

Step 3: filter out the honeycomb (hereinafter to be referred as territory honeycomb) that angle point is positioned at cross section geometric model area in K (0≤K≤Kp) ring honeycomb, in note K ring territory, the quantity of honeycomb is nf;

First, each honeycomb central point C ' (x to K (0≤K≤Kp) ring _i, y _i), write out its corresponding honeycomb j (j=1,2 ..., 6) coordinate of individual angle point:

$\left\{\begin{array}{c}{x}_{i,j}=x_i+h\cos \frac{(j-1)\mathrm{\&pi;}}{3}\\ y_{i,j}=y_i+h\sin \frac{(j-1)\mathrm{\&pi;}}{3}\end{array}\right.j=\mathrm{1,2},...,6---\left(11\right)$

Wherein, the length of side that h is honeycomb;

Then, to the j of this honeycomb (j=1,2 ..., 6) individual angle point (x _{i, j}, y _{i, j}), judge whether it is positioned at cross section geometric model area, if this honeycomb has angle point to be positioned at cross section geometric model area, the center point coordinate of this honeycomb is recorded in a matrix cmf, finally obtain honeycomb number, i.e. the line number nf of matrix cmf in the territory of K ring;

Judge 1 P (x _p, y _p) whether be positioned at geometric model region and can use numerical inverse isoparametric mapping interpolation or triangle area method, two kinds of applicable geometric models of method are all plane quadrilateral regions, therefore be a plurality of quadrilateral areas by the cross section geometric model partition of special-shaped honeycomb in step 1 of the present invention, for these a plurality of quadrilateral areas, put respectively the judgement that whether is positioned at its territory, if exist such region to make a little within it in the model of cross section, explanation point is positioned at cross section geometric model area;

Wherein, numerical inverse isoparametric mapping interpolation is according to the analyticity confrontation quadrilateral area of quadrilateral area and puts P (x _p, y _p) carry out coordinate transform, after coordinate transform, quadrilateral area becomes D{ (ξ, η) |-1≤ξ≤1 ,-1≤η≤1}, and the coordinate of point becomes P ' (ξ _p, η _p), if some P is (ξ _p, η _p) be located at region D{ (ξ, η) |-1≤ξ≤1, within-1≤η≤1}, an explanation point P (x _p, y _p) be positioned at former quadrilateral area; Triangle area method is with a P (x _p, y _p) be four triangles of four edges structure of summit and quadrilateral area Ω to be judged, if the area that four leg-of-mutton area sums are less than or equal to quadrilateral area Ω illustrates a P (x _p, y _p) be positioned at region Ω;

Step 4: if nf < 1, K (0≤K≤Kp) encircles the not interior honeycomb of existence domain, illustrate that the honeycomb central point having generated has covered the cross section geometric model area of special-shaped honeycomb completely, and then execution step 11, otherwise continue execution step five, the effect of this step is the waste of avoiding computational resource;

Step 5: generate K (0≤K≤Kp) ring honeycomb i (i=1 ..., the nf) node of honeycomb on cross section geometric model in individual territory, and remember that the node number that this honeycomb is positioned at geometric model region is nin1;

This step generates the i (i=1 of K (0≤K≤Kp) ring honeycomb, ..., nf) node of honeycomb on cross section geometric model in individual territory, and each nodal information of honeycomb in this territory is write to a matrix N Inf, NInf is the matrix of 6nd * 5, nd is the cell density on honeycomb limit, and the information that the first row to the of NInf five row will record is respectively: whether angle point, node that whether node ID, node are positioned at honeycomb are positioned at the horizontal ordinate of geometric model region, node, the ordinate of node;

Generate by the following method the node of this honeycomb:

First, for center point coordinate, be (x _i, y _i) honeycomb, each node coordinate value of its correspondence is:

$\left\{\begin{array}{c}{x}_{i,j}=\frac{\sin \frac{\mathrm{\&pi;}}{3}h}{\sin (\frac{2\mathrm{\&pi;}}{3}-\mathrm{\&alpha;})}\cos \mathrm{\&theta;}+x_i\\ y_{i,j}=\frac{\sin \frac{\mathrm{\&pi;}}{3}h}{\sin (\frac{2\mathrm{\&pi;}}{3}-\mathrm{\&alpha;})}\sin \mathrm{\&theta;}+y_i\end{array}\right.j=\mathrm{1,2},...,6\mathrm{nd}---\left(12\right)$

Wherein, the length of side that h is honeycomb, (j=1,2 ..., 6nd),

And, because honeycomb exists bonding thickness on y direction of principal axis, therefore need be by each node and central point (x _i, y _i) distance on y direction of principal axis expands as original Ap doubly, i.e. y ' _{i, j}=y _i+ (y _{i, j}-y _i) Ap, wherein, Ap is the honeycomb width ratio of the bonding thickness of consideration, h is the length of side of honeycomb, and t is the bonding thickness of two honeycomb abutting edges;

Then, to the information that writes each node in matrix N Inf:

The coordinate figure of the node of this honeycomb is write to the 4th row and the 5th row of matrix N Inf, the node ID of this honeycomb (above-mentioned j value) is write to the first row of matrix N Inf, position attribution by the node of this honeycomb in honeycomb writes the secondary series of matrix N Inf, if node be positioned at honeycomb angle point element value be 1, if it is 0 that node is positioned on honeycomb limit element value, can adopt formula m=mod (j-1, nd) carry out the judgement of the position attribution of node in honeycomb, position attribution by the node of this honeycomb in geometric model writes the 3rd row of matrix N Inf, if node be positioned at geometric model region element value be 1, if it is 0 that node is positioned at outside geometric model region element value,

I (i=1 ..., nf), after the node of honeycomb all generates in individual territory, find out the node number that is positioned at geometric model region in this honeycomb, the line number that in matrix N Inf, the 3rd column element value is 1, is denoted as nin1;

Step 6: to the i of K (0≤K≤Kp) ring honeycomb (i=1 ..., nf), honeycomb is processed in individual the territory in, and remembers that the node number that this honeycomb is after treatment positioned at geometric model region is nin2;

This step is intended to processing and is positioned at the incomplete honeycomb on cross section geometric model boundary, it is taked to special node arrangement method, to reach mating of finite element model and geometric model, i (i=1 to K (0≤K≤Kp) ring honeycomb, ..., nf) after in individual territory, honeycomb is processed, the information matrix NInf of this honeycomb has corresponding renewal, the method that the node of honeycomb is processed will after describe in detail;

After processing, i (i=1 ..., nf) in individual territory will there is certain variation in the nodal community of honeycomb, find out the node number that this honeycomb is now positioned at geometric model region, the line number that in matrix N Inf, the 3rd column element value is 1, is denoted as nin2;

Step 7: judge whether nin2 is greater than 2, if be greater than 2, continue execution step eight, otherwise directly perform step nine;

Nin2≤2 explanation i (i=1, ..., nf) in individual territory, honeycomb can be regarded as approx and only has a honeycomb limit to be positioned on cross section geometric model boundary, in actual conditions, this honeycomb limit does not exist, therefore abandon this honeycomb, and then honeycomb in i+1 territory of processing K (0≤K≤Kp) ring;

Step 8: the i of division K (0≤K≤Kp) ring honeycomb (i=1 ..., nf) the line unit of honeycomb in individual territory;

The territory interior nodes of honeycomb is write to matrix nodes in order successively, from matrix nodes, choose suitable node and construct the line unit of honeycomb, the node number that forms line unit is write to a matrix elems, for honeycomb in complete territory, in its honeycomb, every two adjacent nodes generate a line unit, for the incomplete honeycomb being positioned on cross section geometric model boundary, every two adjacent territory interior nodes generate a line unit, but do not generate line unit between first territory interior nodes and last territory interior nodes, be that honeycomb is not had cell wall structures by model boundary cut place,

Because honeycomb abutting edge exists bonding thickness, while carrying out finite element analysis, abutting edge shell unit need to be assigned to other one-tenth-value thickness 1/10, therefore the shell unit of abutting edge need be classified as to a set, with matrix sets, record these unit, this modeling method hypothesis y direction of principal axis is bonding direction, will in each honeycomb, perpendicular to the line unit number on the honeycomb limit of y axle, write matrix sets, concerning complete honeycomb, the order generating by its node, the condition that in honeycomb, n unit is unit, abutting edge is: by (n-1)/nd, to zero value of rounding, be 1 or 4, wherein, nd is the cell density on honeycomb limit,

After line unit on cross section geometric model all generates, matrix nodes stores the whole nodal informations on cross section geometric model, matrix elems stores the whole line unit informations on cross section geometric model, and matrix sets stores the unit number of honeycomb abutting edge on cross section geometric model;

Step 9: i=i+1, and judge whether i value is greater than nf, if be less than or equal to nf, return to execution step five, if be greater than nf, perform step ten;

I > nf illustrates that the finite element model of honeycomb in the territory of K (0≤K≤Kp) ring set up, and then processes K+1 and encircle honeycomb;

Step 10: K=K+1, and judge whether K value is greater than Kp, if be less than or equal to Kp, return to execution step two, if be greater than Kp, illustrate that the finite element model of the maximum number of rings honeycomb of this method regulation all generates, continue execution step 11;

Step 11: output shell unit information;

Obtain cross section geometric model reach the standard grade unit sum ElemNum and node sum NodeNum, suppose has NLay layer shell unit on z direction of principal axis;

For iLay (iLay=1 ..., NLay) layer on j (j=1, ..., ElemNum) individual shell unit, its node number is n1, n2, n1+NodeNum and n2+NodeNum, obviously, during iLay=1, n1 and n2 are the corresponding line cell node number in matrix elems, during iLay > 1, n1 and n2 i.e. the node number of two higher values of iLay-1 layer shell unit;

According to output requirement, the node number that unit is comprised is by row writing unit information output file;

Step 12: output node information;

For the node being positioned on cross section geometric model, its ordinate is 0, the ordinate of other node utilizes the interpolation method of thickness direction to obtain according to special-shaped honeycomb geometric model characteristic, it for ordinate value, is not 0 node, for node coordinate vector coord3 of its structure, the first two element of coord3 is horizontal ordinate and the ordinate value that this node correspondence records in matrix nodes, the 3rd element of coord3, the i.e. ordinate value of this node, the method that four angular coordinates of employing based on geometric model quadrilateral subregion carry out planar interpolation obtains, the final finite element model forming will be the plane of burst at the end face of honeycomb, or more accurate, the quadrilateral subregion at the projection place based on node on x-y face selects corresponding interpolating function to carry out surface interpolation, end face at honeycomb is obtained and the realistic model patch surface of exact matching more,

Coordinate transform by the coordinate of each node through over-angle rotation, obtain the node coordinate under former angle position, by output requirement output, here should be noted that the output order of nodal information, the node being positioned on cross section geometric model should be according to the order of the storage of matrix nodes, other every node layer all should, corresponding to the storage order of matrix nodes, so just can make nodal information and unit information defeated corresponding;

Step 13: the unit set information of output honeycomb abutting edge;

According to output requirement, the unit number of honeycomb abutting edge is write to aggregate information output file.

This self geometrical property according to hexagon honeycomb is carried out the method for finite element modeling, can generate fast the finite element model of a special-shaped honeycomb, and, the generation of finite element model does not depend on entity, do not need to set up the three-dimensional entity model of special-shaped honeycomb, also abandoned traditional solid model of take is basic manual modeling method simultaneously, has greatly improved modeling efficiency.The present invention has considered that honeycomb exists bonding thickness between honeycomb and the angle problem of arranging of honeycomb, and the imperfect honeycomb being positioned on model boundary is processed, and makes finite element model reach the effect matching with solid model.

Introduce the method for in the finite element modeling method step 6 of a kind of special-shaped honeycomb of the present invention, the node of honeycomb being processed below, as shown in Figure 2, concrete implementation step is as follows for its process flow diagram:

Step 1: find out the node number that this honeycomb is positioned at cross section geometric model area, the line number that in the nodal information matrix Inf of this honeycomb, the 3rd column element value is greater than zero, is denoted as nin;

Step 2: judge whether nin equals 6nd, if nin=6nd illustrates that this honeycomb is honeycomb in a complete territory, directly perform step 13, if nin ≠ 6nd illustrates that honeycomb is blocked by cross section geometric model boundary, lost integrality, execution step three;

Step 3: the node to this honeycomb sorts;

The object of this step is: make in the nodal information matrix N Inf of this honeycomb the node of each row representative from top to bottom in honeycomb, become order counterclockwise to arrange, and the first row representative in matrix N Inf is by the first point in counterclockwise tactic cross section geometric model area;

The specific implementation method of this step is: the first row of matrix N Inf is moved to last column, and all the other each row upwards move respectively a line, and this operation is carried out for several times until matrix N Inf meets the requirement of this step;

Step 4: in a row vector pend, obviously, pend represents will enter by counterclockwise order the nodal information of cross section geometric model area by the information storage of last column of matrix N Inf;

Step 5: stipulate an error range herr;

Step 6: process this honeycomb and be positioned at cross section geometric model area by counterclockwise tactic last Nodes;

First, last node in cross section geometric model area is denoted as to P ₁(nin of homography NInf is capable), P ₁next node be denoted as P ₂(nin+1 of homography NInf is capable), tries to achieve P ₁and P ₂the intersection point pin1 of line and cross section geometric model boundary;

Then, investigate P ₁and the distance between pin1, if distance is greater than herr, extra-regional node is moved on on cross section geometric model boundary, the coordinate figure that is about to some pin1 is assigned to the 4th capable row and the 5th row of nin+1 of matrix N Inf, simultaneously, the attribute of the node of the capable representative of nin+1 of matrix N Inf changes " not being positioned at honeycomb angle point " and " being positioned at cross section geometric model area " into, and the capable secondary series element value of nin+1 that is about to matrix N Inf changes 0, the three column element value into and changes 1 into;

If P ₁and the distance between pin1 is less than or equal to herr, consider the situation that nin is 1, this illustrates that this honeycomb only has a node in region, and some P in boundary ₁with the distance of intersection point pin1 within error range, this honeycomb is given up, turn back to the step 12 of this method main flow, if nin is not 1, last node in region is moved on on cross section geometric model boundary, the coordinate figure that is about to some pin1 is assigned to the 4th capable row and the 5th row of nin of matrix N If, simultaneously, the attribute of the node of the capable representative of nin of matrix N If changes " not being positioned at honeycomb angle point " and " being positioned at cross section geometric model area " into, the capable secondary series element value of nin that is about to matrix N Inf changes 0, the three column element value into and changes 1 into;

Step 7: find out the line number that the 3rd column element value of matrix N Inf is greater than 0, and be again assigned to nin, this represents that this honeycomb is through the operation of above step, and the node number nin that this honeycomb is positioned at cross section geometric model area upgrades;

Step 8: process this honeycomb and be positioned at cross section geometric model area by counterclockwise tactic first Nodes;

First, node row vector pend being recorded is denoted as P ₁, by cross section geometric model area first be denoted as P ₂(the 1st row of homography NInf), the intersection point of trying to achieve honeycomb and cross section geometric model boundary is pin2;

Then, investigate P ₂and the distance between pin2, if distance is greater than herr, consider that nin equals the situation of the line number of matrix N Inf, now this honeycomb all becomes point in territory through all nodes of operation of step 6, need to increase 6nd+1 node in last interpolation a line of matrix N Inf, the content of this line can temporarily copy the content of last column of matrix N Inf, if nin is not equal to the line number of matrix N Inf, do not need to do the operation that increases row;

Last column of matrix N Inf is moved to the first row, all the other each row move respectively a line backward, the first row of matrix N Inf is now represented to the information of first node in territory, the coordinate figure of pin2 is assigned to the 4th row and the 5th row of matrix N Inf the first row, simultaneously, the nodal community of matrix N Inf the first row representative changes " not being positioned at honeycomb angle point " and " being positioned at cross section geometric model area " into, and the secondary series element value that is about to the first row of matrix N Inf changes 0, the three column element value into and changes 1 into;

If P ₂and the distance between pin2 is less than or equal to herr, first node in region is moved on on cross section geometric model boundary, the coordinate figure that is about to pin2 is assigned to the 4th row and the 5th row of matrix N Inf the first row, simultaneously, the nodal community of matrix N Inf the first row representative changes " not being positioned at honeycomb angle point " and " being positioned at cross section geometric model area " into, the secondary series element value that is about to the first row of matrix N Inf changes 0, the three column element value into and changes 1 into;

Step 9: find out the line number that the 3rd column element value of matrix N Inf is greater than 0, and be again assigned to nin, this represents that this honeycomb is through the operation of above-mentioned steps, and the node number nin that this honeycomb is positioned at cross section geometric model area upgrades;

Step 10: be less than herr if honeycomb is positioned at first node and last nodal distance of cross section geometric model area, carry out node merging;

Investigate the represented node of the 1st row of matrix N Inf and the distance value between the capable represented node of nin, if distance is less than herr, the represented node of the 1st row of matrix N Inf is moved on to the midpoint of two nodes, soon the 4th of the first row of matrix N Inf the row and the 5th column element value change the mid point coordinate figure of two nodes into, simultaneously, the attribute of nin node is changed into " being positioned at cross section geometric model area outer ", be about to the capable tertial element value of nin and change 0 into;

Step 11: find out and be positioned at cross section geometric model area in this honeycomb and be positioned at the node number on honeycomb angle point, be that in matrix N Inf, secondary series and the 3rd column element value are all 1 line number, be denoted as ncd, if ncd < 1, illustrate that this honeycomb does not have angle point to be positioned at cross section geometric model area, returns to the step 12 in this method main flow, if ncd >=1, find out the line number that the 3rd column element value of matrix N Inf is greater than 0, and be again assigned to nin;

Step 12: judge the pattern of this honeycomb and carry out cellular node sequence;

First, remember that this honeycomb first node being arranged on honeycomb angle point in cross section geometric model area is positioned at pos at matrix N Inf ₁oK, and remember d ₁=pos ₁-1, remember that the node that last in cross section geometric model area of this honeycomb is arranged on honeycomb angle point is positioned at pos at matrix N Inf ₂oK, and remember d ₂=nin-pos ₂;

Then, if d ₁> d ₂, need the node of this honeycomb to become following order: the last point in former territory is the first point (the first row of homography NInf) in territory, and each node is by order arrangement clockwise;

Step 13: find out the line number that the 3rd column element value of matrix N Inf is greater than 0, and be again assigned to nin, this represents that this honeycomb is through the operation of above-mentioned steps, and the node number nin that this honeycomb is positioned at cross section geometric model area upgrades.

When honeycomb is processed in this method, need to ask the intersecting point coordinate of two line segments of 4 formations, concrete grammar is as follows:

First, suppose two line segments with end points coordinate be respectively P ₁(x ₁, y ₁), P ₂(x ₂, y ₂), P ₃(x ₃, y ₃), P ₄(x ₄, y ₄);

Then, structure following parameters:

$\left\{\begin{array}{c}{d}_1=-(x_3-x_4)(y_1-y_2)+(x_1-x_2)(y_3-y_4)\\ d_2=-(x_3-x_4)(y_1-y_2)+(x_1-x_2)(y_3-y_4)\end{array}\right.---\left(13\right)$

If d ₁d ₂=0, two line segments are parallel, without intersection point, otherwise, the intersection point of two line segments (or intersection point of extended line) P ₀the coordinate figure of (x, y) is:

$\left\{\begin{array}{c}x=\frac{x_1[x_4(-y_2+y_3)+x_3(y_2-y_4)]+x_2[x_4(y_1-y_3)+x_3(-y_1+y_4)]}{d_1}\\ y=\frac{x_4(y_1-y_2)y_3+x_1{y}_2(y_3-y_4)+x_3(-y_1+y_2)y_4+x_2{y}_1(-y_3+y_4)}{d_2}\end{array}\right.---\left(14\right)$

The finite element modeling method of a kind of special-shaped honeycomb of the present invention is described below by an example.

The honeycomb of certain excision forming, sets up its finite element model, the geometric parameter of honeycomb as shown in Figure 3, bonding direction and y direction of principal axis angle at 45 ° between honeycomb in this structure, between honeycomb, bonding thickness is 0.12mm, honeycomb width is 15mm;

This routine finite element modeling development platform is mathematical software MatLab, in platform, set up the geometric model of this excision forming honeycomb bottom surface, according to the geometrical property of this structure, bottom surface geometric model is divided into two subregions, these two subregions are all quadrilaterals, and end face corresponding to subregion is all smooth flat;

The bonding direction of honeycomb and y direction of principal axis angle at 45 ° in this example, therefore by the geometric model 45° angle that turns clockwise, select the reference point of finite element model as rotary middle point, as shown in Figure 4, the bottom surface geometric model that this example is set up in platform Matlab as shown in Figure 5 for reference point location;

Point centered by reference point, to outer ring-like, generate the line unit of honeycomb, to each honeycomb, first generate its honeycomb central point, and judge whether this honeycomb has angle point to be positioned at geometric model region, bottom surface, this example adopts numerical inverse isoparametric mapping interpolation to judge whether honeycomb angle point is positioned at geometric model region, bottom surface, to angle point P ₀(x _p, y _p) and one of them quadrilateral subregion P of this example ₁p ₂p ₃p ₄, calculate as follows:

The coordinate of four angle points of note quadrilateral subregion is P ₁(x ₁, y ₁), P ₂(x ₂, y ₂), P ₃(x ₃, y ₃), P ₄(x ₄, y ₄), constructing variable

$\left\{\begin{array}{c}A=(x_1-x_4)(y_2-y_3)+x_3(y_1-y_4)+x_2(-y_1+y_4)\\ B=x_4{y}_2+x_1{y}_3-{2x}_4{y}_3+x_p(y_1-y_2+y_3-y_4)-x_2{y}_4-x_1{y}_p+x_2{y}_p+x_4{y}_p-x_3(y_1-2y_4+y_p)\\ D=x_p(-y_3+y_4)+x_4(y_3-y_p)+x_3(-y_4+y_p)\end{array}\right.$

(15)

If A=0, note np=1-2r ₁if, A ≠ 0 and B ²-4AD>=0, has

$\left\{\begin{array}{c}{r}_{11}=-\frac{B+\sqrt{B^2-4\mathrm{AD}}}{2A}\\ r_{12}=-\frac{B-\sqrt{B^2-4\mathrm{AD}}}{2A}\end{array}\right.---\left(16\right)$

$\left\{\begin{array}{c}{\mathrm{np}}_1=1-{2r}_{11}\\ {\mathrm{np}}_2=1-{2r}_{12}\end{array}\right.---\left(17\right)$

$\left\{\begin{array}{c}{E}_1=|y_p-(\frac{1-{\mathrm{np}}_1}{4}{y}_1+\frac{1-{\mathrm{np}}_1}{4}{y}_2+\frac{1+{\mathrm{np}}_1}{4}{y}_3+\frac{1+{\mathrm{np}}_1}{4}{y}_4)|\\ E_2=|y_p-(\frac{1-{\mathrm{np}}_2}{4}{y}_1+\frac{1-{\mathrm{np}}_2}{4}{y}_2+\frac{1+{\mathrm{np}}_2}{4}{y}_3+\frac{1+{\mathrm{np}}_2}{4}{y}_4)|\end{array}\right.---\left(18\right)$

If E ₁< E ₂, by np=np ₁, otherwise, by np=np ₂, constructing variable

$\left\{\begin{array}{c}{A}_1=x_4(y_1-y_2)+x_3(-y_1+y_2)+(x_1-x_2)(y_3-y_4)\\ B_1=x_4{y}_2-x_1{y}_3+2x_2{y}_3-x_2{y}_4+x_p(-y_1+y_2-y_3+y_4)+x_1{y}_p-x_2{y}_p-x_4{y}_p+x_3(y_1-2y_2+y_p)\\ D_1=x_p(-y_2+y_3)+x_3(y_2-y_p)+x_2(-y_3+y_p)\end{array}\right.$

(19)

If A ₁=0, note ep=1-2r ₂if, A ₁≠ 0 and have

$\left\{\begin{array}{c}{r}_{21}=-\frac{B_1+\sqrt{B_1^2-4A_1{D}_1}}{2A_1}\\ r_{22}=-\frac{B_1-\sqrt{B_1^2-4A_1{D}_1}}{2A_1}\end{array}\right.---\left(20\right)$

$\left\{\begin{array}{c}{\mathrm{ep}}_1=1-2r_{21}\\ {\mathrm{ep}}_2=1-{2r}_{22}\end{array}\right.---\left(21\right)$

$\left\{\begin{array}{c}{E}_1=|x_p-[\frac{(1-{\mathrm{ep}}_1)(1-\mathrm{np})}{4}{x}_1+\frac{(1+{\mathrm{ep}}_1)(1-\mathrm{np})}{4}{x}_2+\frac{(1+{\mathrm{ep}}_1)(1+\mathrm{np})}{4}{x}_3+\frac{(1-{\mathrm{ep}}_1)(1+\mathrm{np})}{4}{x}_4\left]\right|\\ E_2={|x}_p-[\frac{(1-{\mathrm{ep}}_2)(1-\mathrm{np})}{4}{x}_1+\frac{(1+{\mathrm{ep}}_2)(1-\mathrm{np})}{4}{x}_2+\frac{(1+{\mathrm{ep}}_2)(1+\mathrm{np})}{4}{x}_3+\frac{(1-{\mathrm{ep}}_2)(1+\mathrm{np})}{4}{x}_4]|\end{array}\right.$

(22)

If E ₁< E ₂, by ep=ep ₁, otherwise, by ep=ep ₂, some P ₀point coordinate after coordinate transform is P ₀₁(ep, np), if-1≤ep≤1, and-1≤np≤1, explanation point P ₀be positioned at a P ₁, P ₂, P ₃, P ₄in quadrilateral area for angle point;

If this honeycomb has angle point to be positioned at geometric model region, bottom surface, this honeycomb is honeycomb in a territory, and it is carried out to the operation of line dividing elements;

Before dividing line unit, need process honeycomb, the node of determining this honeycomb position of arranging, in this example, the cell density of setting on honeycomb limit is nd=2, and the nodal information of this honeycomb is recorded in matrix N Inf in order;

The territory interior nodes of honeycomb is write to matrix nodes in order successively, soon the 1st row of the nodal information matrix N Inf of each honeycomb is capable to nin2, the 4th row write to the 5th row, obtain the most at last storing the matrix nodes of all node coordinates on this routine bottom surface geometric model;

The node number that each line unit is comprised writes matrix elems, in honeycomb, every two adjacent territory interior nodes will generate a line unit, incomplete honeycomb is tectonic line unit not in cut place, obtain the most at last storing on this routine bottom surface geometric model the matrix elems of wired unit information;

The line unit number of abutting edge is write among matrix sets, in this example, the method of finding abutting edge line unit number is: remember ie (ie=1 in this honeycomb, ..., ne) node ID of individual node (being the value of the capable first row of ie in matrix N Inf) is n1, the more lower node ID of remembering ie node is n2 (if ie node has been last node in territory, the sequence number of next node is 1), constructing variable s ₁=n1/nd and s ₂=n2/nd, if parameter s ₁and s ₂fix (s satisfies condition simultaneously ₁-2)=0 or fix (s ₁-5)=0 and condition fix (s ₂-2)=0 or fix (s ₂-5)=0, the ie in this honeycomb (ie=1 ..., ne) individual line unit is positioned at abutting edge, and this unit number is write to matrix sets, obtains the most at last storing the matrix sets of all abutting edges line unit number on this routine bottom surface geometric model;

The line unit, bottom surface of this routine excision forming honeycomb generates situation as shown in Figure 6;

In this example, on this routine excision forming honeycomb thickness direction, generate three layers of shell unit, by the planar interpolation of thickness direction, obtain the ordinate of space nodes, for a node P (x, y) who is stored in matrix nodes, judge that it is positioned at a quadrilateral subregion P ₁p ₂p ₃p ₄within, the coordinate of four angle points of note quadrilateral subregion is P ₁(x ₁, y ₁, z ₁), P ₂(x ₂, y ₂, z ₂), P ₃(x ₃, y ₃, z ₃), P ₄(x ₄, y ₄, z ₄), according to this quadrilateral subregion, adopting numerical inverse isoparametric mapping interpolation is P ' (ep, np) by the coordinate transform of putting P (x, y), constructing variable

$\left\{\begin{array}{c}{n}_1=\frac{(1-\mathrm{ep})(1-\mathrm{np})}{4}\\ N_2=\frac{(1+\mathrm{ep})(1-\mathrm{np})}{4}\\ N_3=\frac{(1+\mathrm{ep})(1+\mathrm{np})}{4}\\ N_4=\frac{(1-\mathrm{ep})(1+\mathrm{np})}{4}\end{array}\right.---\left(23\right)$

Then, the ordinate value that draws the corresponding model end face of this node is Ptmp=N ₁z ₁+ N ₂z ₂+ N ₃z ₃+ N ₄z ₄, iLay (iLay=1 ..., NLay) layer node ordinate be

The message file that this routine output can be read by finite element pre-processing software Hypermesh, first export the node coordinate value information of this routine honeycomb, each node is done to the rotation contrary with geometric model, obtain the node coordinate value under former angle, write output file, the message file form of Hypermesh can be " .cmf " form, each node coordinate takies a line, information format is " * createnode (x, y, z) ";

Then export the node number information that the shell unit of this routine honeycomb comprises, for iLay (iLay=1, ..., NLay) j shell unit on layer, its node number is n1, n2, n1+NodeNum and n2+NodeNum, during iLay=1, n1 and n2 are the corresponding line cell node number in matrix elems, during iLay > 1, n1 and n2 i.e. the node number of two higher values of iLay-1 layer shell unit, wherein NodeNum is the node number on the geometric model of bottom surface, the node number information of shell unit is write to output file, similarly, file is " .cmf " form, in unit information output file, to each unit, should write two row information, this example generates the shell unit of four nodes, the first row information is " * createlist (nodes, 1) i j k m ", wherein, i, j, k and m are the node numbers that this unit comprises, the second row information is " * createelement (104, 1, 1, 1) ",

Finally, export the shell unit aggregate information of this routine honeycomb abutting edge, the file of output " .cmf " form, file header is written as " * createmark (elements; 1) ", obtains the line number k of matrix sets, to the i (i=1 in matrix sets, ..., k) the unit number sets of row representative _i, obtain its correspondence j (j=1 ..., NLay) floor shell unit NLay (sets _i-1)+j, and writing information file, finally write the input that " * entitysetcreate (E_HH, elements, 1) " completes shell unit aggregate information;

By nodal information, shell unit information and shell unit aggregate information are read in software Hypermesh, can obtain the finite element model of this excision forming honeycomb, the generation situation of the finite element model that Fig. 7 is this routine honeycomb in Hypermesh.

For the special-shaped honeycomb of any given geometric parameter, utilize method provided by the invention can carry out fast grid division, by setting different parameters, obtain satisfactory honeycomb finite element model, there is rapidity, the feature of universality.

Claims (2)

1. a finite element modeling method for special-shaped honeycomb, is characterized in that: the method concrete steps are as follows:

Step 1: the cross section geometric model of setting up special-shaped honeycomb;

First, according to the border characteristic of given special-shaped honeycomb, select a cross section perpendicular to honeycomb thickness direction, set up the cross section geometric model of this honeycomb, the cross section geometric model of honeycomb is based upon on the x-y face of cartesian coordinate system, and the thickness direction of honeycomb is positioned at z direction of principal axis;

Then, cross section geometric model partition is become to the quadrilateral subregion of a plurality of regular shape, here need repeatedly to put the judgement that whether is positioned at plane geometry model area, and the interpolation calculation of honeycomb thickness direction, these two kinds of calculating are all based on quadrangular plan geometric areas, therefore divide according to model actual conditions pair cross-section geometric model;

Finally, according to the geometrical property of hexagon honeycomb, carry out finite element modeling, suppose that the bonding direction of honeycomb is positioned at y direction of principal axis, if desired the bonding direction of honeycomb is positioned at other direction, the cross section geometric model of honeycomb to be rotated to respective angles in coordinate system and carry out finite element modeling, after modeling completes, finite element model be rotated back to former angle, wherein, with 1 P ₀(x ₀, y ₀) be rotary middle point, will put P (x, y) and be rotated counterclockwise θ angle and obtain a P _t(x _t, y _t) the computing of coordinate figure be expressed in matrix as:

$\left[\begin{array}{cc}{x}_t& y_t\end{array}\right]=\left[\begin{array}{cc}{x}_0& y_0\end{array}\right]+\left(\begin{array}{c}\left[\begin{array}{cc}x& y\end{array}\right]-\left[\begin{array}{cc}{x}_0& y_0\end{array}\right]\end{array}\right)\left[\begin{array}{cc}\cos \mathrm{\&theta;}& \sin \mathrm{\&theta;}\\ -\sin \mathrm{\&theta;}& \cos \mathrm{\&theta;}\end{array}\right]---\left(1\right)$

Step 2: generate the honeycomb central point of K ring honeycomb, K value is 0≤K≤Kp herein;

Select a honeycomb central point C (x on cross section geometric model _c, y _c) as reference point, outwards generating the honeycomb central point of K ring honeycomb, K value is 0≤K≤Kp herein, centered by C point, the honeycomb of point is considered as the 0th ring here, is outwards followed successively by the 1st, 2 ..., Kp ring, Kp is the default honeycomb number of rings upper limit;

When K=0, the honeycomb central point on K ring is C (x _c, y _c), when K ≠ 0, generate as follows K ring honeycomb central point:

First, on generation K ring, honeycomb central point and C are ordered line and x axle clamp angle are the honeycomb central point of (i=1,2..., 6), K value is 1≤K≤Kp herein, and so every ring of honeycomb central point has six, and its coordinate is:

$\left\{\begin{array}{c}{x}_{k,i}=\sqrt{3}\mathrm{Kh}\cos [\frac{\mathrm{\&pi;}}{6}+(i-1)\frac{\mathrm{\&pi;}}{3}]+x_C\\ y_{k,i}=\sqrt{3}\mathrm{Kh}\sin [\frac{\mathrm{\&pi;}}{6}+(i-1)\frac{\mathrm{\&pi;}}{3}]+y_C\end{array}\right.i=\mathrm{1,2},...,6---\left(2\right)$

Wherein, the length of side that h is honeycomb;

Then, consider that honeycomb exists bonding thickness on y direction of principal axis, by the above-mentioned central point generating with the y direction of principal axis between reference point C apart from expanding Ap times, i.e. y' _k,i=y _c+ Ap (y _k,i-y _c), wherein, Ap is the honeycomb width ratio of the bonding thickness of consideration, t is the thickness of two honeycomb abutting edges, and the actual (real) thickness of two honeycomb abutting edges and honeycomb cell-wall thickness is poor;

Finally, generate remaining honeycomb central point on K ring, K value is 2≤K≤Kp herein, and these honeycomb central points are positioned on the Along ent of every pair of adjacent line of generating center point, choose the two adjacent central point C that generated _k,i(x _k,i, y' _k,i) and C _{k, i+1}(x _{k, i+1}, y' _{k, i+1}), the honeycomb central point C by this two dot generation between the two _k,jcoordinate figure be:

$\left\{\begin{array}{c}{x}_{k,j}=x_{k,i}(1-q)+x_{k,i+1}q\\ y_{k,j}={y^{\&prime;}}_{k,i}(1-q)+{y^{\&prime;}}_{k,i+1}q\end{array}\right.j=1,...,K-1---\left(3\right)$

Wherein, q=j/K, such operation is carried out six times altogether, and remaining honeycomb central point on K ring is generated;

Step 3: filtering out the honeycomb that angle point is positioned at cross section geometric model area in K ring honeycomb is territory honeycomb, and K value is 0≤K≤Kp herein, and in note K ring territory, the quantity of honeycomb is nf;

First, each honeycomb central point C'(x to K ring _i, y _i), K value is 0≤K≤Kp herein, writes out the coordinate of j angle point of its corresponding honeycomb:

$\left\{\begin{array}{c}{x}_{i,j}=x_i+h\cos \frac{(j-1)\mathrm{\&pi;}}{3}\\ y_{i,j}=y_i+h\sin \frac{(j-1)\mathrm{\&pi;}}{3}\end{array}\right.j=\mathrm{1,2},...,6---\left(4\right)$

Wherein, the length of side that h is honeycomb;

Then, j the angle point (x to this honeycomb _i,j, y _i,j), judge whether it is positioned at cross section geometric model area, if this honeycomb has angle point to be positioned at cross section geometric model area, the center point coordinate of this honeycomb is recorded in a matrix cmf, finally obtain honeycomb number, i.e. the line number nf of matrix cmf in the territory of K ring;

Judge 1 P (x _p, y _p) whether be positioned at geometric model region and use numerical inverse isoparametric mapping interpolation or triangle area method, two kinds of applicable geometric models of method are all plane quadrilateral regions, therefore be a plurality of quadrilateral areas by the cross section geometric model partition of special-shaped honeycomb in step 1, for these a plurality of quadrilateral areas, put respectively the judgement that whether is positioned at its territory, if exist such region to make a little within it in the model of cross section, explanation point is positioned at cross section geometric model area;

Wherein, numerical inverse isoparametric mapping interpolation is according to the analyticity confrontation quadrilateral area of quadrilateral area and puts P (x _p, y _p) carry out coordinate transform, after coordinate transform, quadrilateral area becomes D{ (ξ, η) |-1≤ξ≤1 ,-1≤η≤1}, and the coordinate of point becomes P'(ξ _p, η _p), if some P is (ξ _p, η _p) be located at region D{ (ξ, η) |-1≤ξ≤1, within-1≤η≤1}, an explanation point P (x _p, y _p) be positioned at former quadrilateral area; Triangle area method is with a P (x _p, y _p) be four triangles of four edges structure of summit and quadrilateral area Ω to be judged, if the area that four leg-of-mutton area sums are less than or equal to quadrilateral area Ω illustrates a P (x _p, y _p) be positioned at region Ω;

Step 4: if nf<1, K encircles the not interior honeycomb of existence domain, K value is 0≤K≤Kp herein, illustrate that the honeycomb central point having generated has covered the cross section geometric model area of special-shaped honeycomb completely, and then execution step 11, otherwise continue execution step five, the effect of this step is the waste of avoiding computational resource;

Step 5: generate the node of honeycomb on cross section geometric model in i territory of K ring honeycomb, K value is 0≤K≤Kp herein, and i value is i=1 herein ..., nf, and remember that the node number that this honeycomb is positioned at geometric model region is nin1;

This step generates the node of honeycomb on cross section geometric model in i territory of K ring honeycomb, K value is 0≤K≤Kp herein, i value is i=1 herein, ..., nf, and each nodal information of honeycomb in this territory is write to a matrix N Inf, NInf is the matrix of 6nd * 5, nd is the cell density on honeycomb limit, and the information that the first row to the of NInf five row will record is respectively: whether angle point, node that whether node ID, node are positioned at honeycomb are positioned at the horizontal ordinate of geometric model region, node, the ordinate of node;

Generate by the following method the node of this honeycomb:

First, for center point coordinate, be (x _i, y _i) honeycomb, each node coordinate value of its correspondence is:

$\left\{\begin{array}{c}{x}_{i,j}=\frac{\sin \frac{\mathrm{\&pi;}}{3}h}{\sin (\frac{2\mathrm{\&pi;}}{3}-\mathrm{\&alpha;})}\cos \mathrm{\&theta;}+x_i\\ y_{i,j}=\frac{\sin \frac{\mathrm{\&pi;}}{3}h}{\sin (\frac{2\mathrm{\&pi;}}{3}-\mathrm{\&alpha;})}\sin \mathrm{\&theta;}+y_i\end{array}\right.j=\mathrm{1,2},...,6\mathrm{nd}---\left(5\right)$

Wherein, the length of side that h is honeycomb, $\mathrm{\&theta;}=\frac{(j-1)\mathrm{\&pi;}}{3\mathrm{nd}}(j=\mathrm{1,2},...,6\mathrm{nd}),\mathrm{\&alpha;}=\mathrm{mod}(\mathrm{\&theta;},\frac{\mathrm{\&pi;}}{3});$

And, because honeycomb exists bonding thickness on y direction of principal axis, therefore need be by each node and central point (x _i, y _i) distance on y direction of principal axis expands as original Ap doubly, i.e. y' _i,j=y _i+ (y _i,j-y _i) Ap, wherein, Ap is the honeycomb width ratio of the bonding thickness of consideration, h is the length of side of honeycomb, and t is the bonding thickness of two honeycomb abutting edges;

Then, to the information that writes each node in matrix N Inf:

The coordinate figure of the node of this honeycomb is write to the 4th row and the 5th row of matrix N Inf, by the node ID of this honeycomb, it is the first row that above-mentioned j value writes matrix N Inf, position attribution by the node of this honeycomb in honeycomb writes the secondary series of matrix N Inf, if node be positioned at honeycomb angle point element value be 1, if it is 0 that node is positioned on honeycomb limit element value, employing formula m=mod (j-1, nd) carry out the judgement of the position attribution of node in honeycomb, position attribution by the node of this honeycomb in geometric model writes the 3rd row of matrix N Inf, if node be positioned at geometric model region element value be 1, if it is 0 that node is positioned at outside geometric model region element value,

After in i territory, the node of honeycomb all generates, i value is i=1 herein ..., nf, finds out the node number that is positioned at geometric model region in this honeycomb, and the line number that in matrix N Inf, the 3rd column element value is 1, is denoted as nin1;

Step 6: honeycomb in i territory of K ring honeycomb is processed, and K value is 0≤K≤Kp herein, and i value is i=1 herein ..., nf, and remember that the node number that this honeycomb is after treatment positioned at geometric model region is nin2;

This step is intended to processing and is positioned at the incomplete honeycomb on cross section geometric model boundary, it is taked to node arrangement method, to reach mating of finite element model and geometric model, after in i territory of K ring honeycomb, honeycomb is processed, K value is 0≤K≤Kp herein, and i value is i=1 herein, ..., nf, the information matrix NInf of this honeycomb has corresponding renewal, and the method that the node of honeycomb is processed will describe in detail in addition;

After processing, in i territory, the nodal community of honeycomb will change, and i value is i=1 herein, ..., nf, finds out the node number that this honeycomb is now positioned at geometric model region, be the line number that in matrix N Inf, the 3rd column element value is 1, be denoted as nin2;

Step 7: judge whether nin2 is greater than 2, if be greater than 2, continue execution step eight, otherwise directly perform step nine;

In i territory of nin2≤2 explanations, honeycomb is regarded as approx and is only had a honeycomb limit to be positioned on cross section geometric model boundary, i value is i=1 herein, ..., nf, in actual conditions, this honeycomb limit does not exist, therefore abandon this honeycomb, and then honeycomb in i+1 territory of processing K ring, K value is 0≤K≤Kp herein;

Step 8: in i territory of division K ring honeycomb the line unit of honeycomb, and K value is 0≤K≤Kp herein, and i value is i=1 herein ..., nf;

The territory interior nodes of honeycomb is write to matrix nodes in order successively, from matrix nodes, choose suitable node and construct the line unit of honeycomb, the node number that forms line unit is write to a matrix elems, for honeycomb in complete territory, in its honeycomb, every two adjacent nodes generate a line unit, for the incomplete honeycomb being positioned on cross section geometric model boundary, every two adjacent territory interior nodes generate a line unit, but do not generate line unit between first territory interior nodes and last territory interior nodes, be that honeycomb is not had cell wall structures by model boundary cut place,

Because honeycomb abutting edge exists bonding thickness, while carrying out finite element analysis, abutting edge shell unit need to be assigned to other one-tenth-value thickness 1/10, therefore the shell unit of abutting edge need be classified as to a set, with matrix sets, record these unit, suppose that y direction of principal axis is bonding direction, will in each honeycomb, perpendicular to the line unit number on the honeycomb limit of y axle, write matrix sets, concerning complete honeycomb, the order generating by its node, the condition that in honeycomb, n unit is unit, abutting edge is: by (n-1)/nd, to zero value of rounding, be 1 or 4, wherein, nd is the cell density on honeycomb limit,

After line unit on cross section geometric model all generates, matrix nodes stores the whole nodal informations on cross section geometric model, matrix elems stores the whole line unit informations on cross section geometric model, and matrix sets stores the unit number of honeycomb abutting edge on cross section geometric model;

Step 9: i=i+1, and judge whether i value is greater than nf, if be less than or equal to nf, return to execution step five, if be greater than nf, perform step ten;

I>nf illustrates that the finite element model of honeycomb in the territory of K ring set up, and K value is 0≤K≤Kp herein, and then processes K+1 and encircle honeycomb;

Step 10: K=K+1, and judge whether K value is greater than Kp, if be less than or equal to Kp, return to execution step two, if be greater than Kp, the finite element model of the maximum number of rings honeycomb of regulation all generates, and continues execution step 11;

Step 11: output shell unit information;

Obtain cross section geometric model reach the standard grade unit sum ElemNum and node sum NodeNum, suppose has NLay layer shell unit on z direction of principal axis;

For j shell unit on iLay layer, iLay value is iLay=1 herein ..., NLay, j value is j=1 herein ..., ElemNum, its node number is n1, n2, n1+NodeNum and n2+NodeNum, obviously, during iLay=1, n1 and n2 are the corresponding line cell node number in matrix elems, during iLay>1, n1 and n2 i.e. the node number of two higher values of iLay-1 layer shell unit;

According to output requirement, the node number that unit is comprised is by row writing unit information output file;

Step 12: output node information;

For the node being positioned on cross section geometric model, its ordinate is 0, the ordinate of other node utilizes the interpolation method of thickness direction to obtain according to special-shaped honeycomb geometric model characteristic, it for ordinate value, is not 0 node, for node coordinate vector coord3 of its structure, the first two element of coord3 is horizontal ordinate and the ordinate value that this node correspondence records in matrix nodes, the 3rd element of coord3, the i.e. ordinate value of this node, the method that four angular coordinates of employing based on geometric model quadrilateral subregion carry out planar interpolation obtains, the final finite element model forming will be the plane of burst at the end face of honeycomb, or the quadrilateral subregion that is more accurately the projection place on x-y face based on node selects corresponding interpolating function to carry out surface interpolation, end face at honeycomb is obtained and the realistic model patch surface of exact matching more,

Coordinate transform by the coordinate of each node through over-angle rotation, obtain the node coordinate under former angle position, by output requirement output, here should be noted that the output order of nodal information, the node being positioned on cross section geometric model should be according to the order of the storage of matrix nodes, other every node layer all should, corresponding to the storage order of matrix nodes, so just can make nodal information and unit information defeated corresponding;

Step 13: the unit set information of output honeycomb abutting edge;

According to output requirement, the unit number of honeycomb abutting edge is write to aggregate information output file.

2. the finite element modeling method of a kind of special-shaped honeycomb according to claim 1, is characterized in that: the method that the node to honeycomb described in step 6 is processed, and its concrete implementation step is as follows:

Step 1: find out the node number that this honeycomb is positioned at cross section geometric model area, the line number that in the nodal information matrix N Inf of this honeycomb, the 3rd column element value is greater than zero, is denoted as nin;

Step 2: judge whether nin equals 6nd, if nin=6nd illustrates that this honeycomb is honeycomb in a complete territory, directly perform step 13, if nin ≠ 6nd illustrates that honeycomb is blocked by cross section geometric model boundary, lost integrality, execution step three;

Step 3: the node to this honeycomb sorts;

This step is to make in the nodal information matrix N Inf of this honeycomb the node of each row representative from top to bottom in honeycomb, become order counterclockwise to arrange, and the first row representative in matrix N Inf is by the first point in counterclockwise tactic cross section geometric model area; The first row of matrix N Inf is moved to last column, and all the other each row upwards move respectively a line, and this operation is carried out for several times until matrix N Inf meets the requirement of this step;

Step 4: in a row vector pend, obviously, pend represents will enter by counterclockwise order the nodal information of cross section geometric model area by the information storage of last column of matrix N Inf;

Step 5: stipulate an error range herr;

Step 6: process this honeycomb and be positioned at cross section geometric model area by counterclockwise tactic last Nodes;

First, last node in cross section geometric model area is denoted as to P ₁the nin that is homography NInf is capable, P ₁next node be denoted as P ₂the nin+1 that is homography NInf is capable, tries to achieve P ₁and P ₂the intersection point pin1 of line and cross section geometric model boundary;

Then, investigate P ₁and the distance between pin1, if distance is greater than herr, extra-regional node is moved on on cross section geometric model boundary, the coordinate figure that is about to some pin1 is assigned to the 4th capable row and the 5th row of nin+1 of matrix N Inf, simultaneously, the attribute of the node of the capable representative of nin+1 of matrix N Inf changes " not being positioned at honeycomb angle point " and " being positioned at cross section geometric model area " into, and the capable secondary series element value of nin+1 that is about to matrix N Inf changes 0, the three column element value into and changes 1 into;

If P ₁and the distance between pin1 is less than or equal to herr, consider the situation that nin is 1, this illustrates that this honeycomb only has a node in region, and some P in boundary ₁with the distance of intersection point pin1 within error range, this honeycomb is given up, turn back to the step 12 in claim 1, if nin is not 1, last node in region is moved on on cross section geometric model boundary, the coordinate figure that is about to some pin1 is assigned to the 4th capable row and the 5th row of nin of matrix N Inf, simultaneously, the attribute of the node of the capable representative of nin of matrix N Inf changes " not being positioned at honeycomb angle point " and " being positioned at cross section geometric model area " into, the capable secondary series element value of nin that is about to matrix N Inf changes 0, the three column element value into and changes 1 into;

Step 7: find out the line number that the 3rd column element value of matrix N Inf is greater than 0, and be again assigned to nin, this represents that this honeycomb is through the operation of above step, and the node number nin that this honeycomb is positioned at cross section geometric model area upgrades;

Step 8: process this honeycomb and be positioned at cross section geometric model area by counterclockwise tactic first Nodes;

First, node row vector pend being recorded is denoted as P ₁, by cross section geometric model area first be denoted as P ₂be the 1st row of homography NInf, the intersection point of trying to achieve honeycomb and cross section geometric model boundary is pin2;

Then, investigate P ₂and the distance between pin2, if distance is greater than herr, consider that nin equals the situation of the line number of matrix N Inf, now this honeycomb all becomes point in territory through all nodes of operation of step 6, need to increase 6nd+1 node in last interpolation a line of matrix N Inf, the content of this line temporarily copies the content of last column of matrix N Inf, if nin is not equal to the line number of matrix N Inf, do not need to do the operation that increases row;

Last column of matrix N Inf is moved to the first row, all the other each row move respectively a line backward, the first row of matrix N Inf is now represented to the information of first node in territory, the coordinate figure of pin2 is assigned to the 4th row and the 5th row of matrix N Inf the first row, simultaneously, the nodal community of matrix N Inf the first row representative changes " not being positioned at honeycomb angle point " and " being positioned at cross section geometric model area " into, and the secondary series element value that is about to the first row of matrix N Inf changes 0, the three column element value into and changes 1 into;

If P ₂and the distance between pin2 is less than or equal to herr, first node in region is moved on on cross section geometric model boundary, the coordinate figure that is about to pin2 is assigned to the 4th row and the 5th row of matrix N Inf the first row, simultaneously, the nodal community of matrix N Inf the first row representative changes " not being positioned at honeycomb angle point " and " being positioned at cross section geometric model area " into, the secondary series element value that is about to the first row of matrix N Inf changes 0, the three column element value into and changes 1 into;

Step 9: find out the line number that the 3rd column element value of matrix N Inf is greater than 0, and be again assigned to nin, this represents that this honeycomb is through the operation of above-mentioned steps, and the node number nin that this honeycomb is positioned at cross section geometric model area upgrades;

Step 10: be less than herr if honeycomb is positioned at first node and last nodal distance of cross section geometric model area, carry out node merging;

Investigate the represented node of the 1st row of matrix N Inf and the distance value between the capable represented node of nin, if distance is less than herr, the represented node of the 1st row of matrix N Inf is moved on to the midpoint of two nodes, soon the 4th of the first row of matrix N Inf the row and the 5th column element value change the mid point coordinate figure of two nodes into, simultaneously, the attribute of nin node is changed into " being positioned at cross section geometric model area outer ", be about to the capable tertial element value of nin and change 0 into;

Step 11: find out and be positioned at cross section geometric model area in this honeycomb and be positioned at the node number on honeycomb angle point, be that in matrix N Inf, secondary series and the 3rd column element value are all 1 line number, be denoted as ncd, if ncd<1, illustrate that this honeycomb does not have angle point to be positioned at cross section geometric model area, returns to the step 12 in claim 1, if ncd >=1, find out the line number that the 3rd column element value of matrix N Inf is greater than 0, and be again assigned to nin;

Step 12: judge the pattern of this honeycomb and carry out cellular node sequence;

First, remember that this honeycomb first node being arranged on honeycomb angle point in cross section geometric model area is positioned at pos at matrix N Inf ₁oK, and remember d ₁=pos ₁-1, remember that the node that last in cross section geometric model area of this honeycomb is arranged on honeycomb angle point is positioned at pos at matrix N Inf ₂oK, and remember d ₂=nin-pos ₂;

Then, if d ₁>d ₂, need the node of this honeycomb to become following order: the last point in former territory in territory first be the first row of homography NInf, and each node is arranged by clockwise order;

Step 13: find out the line number that the 3rd column element value of matrix N Inf is greater than 0, and be again assigned to nin, this represents that this honeycomb is through the operation of above-mentioned steps, and the node number nin that this honeycomb is positioned at cross section geometric model area upgrades;

When honeycomb is processed, need to ask the intersecting point coordinate of two line segments of 4 formations, concrete grammar is as follows:

First, suppose two line segments with end points coordinate be respectively P ₁(x ₁, y ₁), P ₂(x ₂, y ₂), P ₃(x ₃, y ₃), P ₄(x ₄, y ₄);

Then, structure following parameters:

$\left\{\begin{array}{c}{d}_1=-(x_3-x_4)(y_1-y_2)+(x_1-x_2)(y_3-y_4)\\ d_2=-(x_3-x_4)(y_1-y_2)+(x_1-x_2)(y_3-y_4)\end{array}\right.---\left(6\right)$

If d ₁d ₂=0, two line segments are parallel, without intersection point, otherwise, the intersection point P of the intersection point of two line segments or extended line ₀the coordinate figure of (x, y) is:

$\left\{\begin{array}{c}x=\frac{x_1[x_4(-y_2+y_3)+x_3(y_2-y_4)]+x_2[x_4(y_1-y_3)+x_3(-y_1+y_4)]}{d_1}\\ y=\frac{x_4(y_1-y_2)y_3+x_1{y}_2(y_3-y_4)+x_3(-y_1+y_2)y_4+x_2{y}_1(-y_3+y_4)}{d_2}\end{array}\right.---\left(7\right).$