Integration of Direction Fields with Standard Options in Finite Element Programs

The two-dimensional differential equation y’ = f(x,y) can be interpreted as a direction field. Commercial Finite Element (FE) programs can be used for this integration task without additional programming, provided that these programs have options for the calculation of orthotropic heat conduction problems. The differential equation to be integrated with arbitrary boundaries is idealized as an FE model with thermal 2D elements. Its orthotropic thermal conductivities are specified as k1 = 1 and k2 = 0. In doing so, k1 is parallel to y ́, and k2 is oriented perpendicular to this. For this extreme case, it is shown that the isotherms are identical to the solution of y’ = f(x,y). The direction fields, for example, can be velocity vectors in fluid mechanics or principal stress directions in structural mechanics. In the case of the latter, possibilities for application in the construction of fiber-reinforced plastics (FRP) arise, since fiber courses, which follow the local principal stress directions, make use of the superior stiffness and strength of the fibers.


Introduction
Orthotropic materials can have extremely different thermal conductivities, for example, in printed circuit boards.The tracks conducting current are made of copper (Cu) and are electrically insulated against the plastic plate.Due to its low thermal conductivity (k 2 ), the Cu conducting track (k 1 ) is heated to a constant, relatively high temperature.It can be considered as an isothermal line.The thermal heat flow in the plastic due to conduction is very low.For example, if there are two Cu conducting tracks with constant but different temperatures, then a steady temperature gradient independent of k 2 prevails in the plastic between them, even if k 2 →0.The Cu-tracks can be described by their position y = y(x), and they control the temperature distribution in the insulator.Alternatively, the position can be described by y'(x) together with their starting positions.If these tracks are infinitely densely distributed, the direction field of the Cu-conductance k 1 is described by y´= f(x,y).Insulation works perpendicular to the Cu-tracks when k 2 << k 1 .Should the (infinitely densely) guided conduction paths be directed parallel to an arbitrary direction field y' = f(x,y), then the following hypothesis shall be mathematically verified: Hypothesis: The isotherms of an orthotropic steady-state 2D thermal conduction problem with the thermal conductivities k 1 and k 2 are tangential to an arbitrarily prescribed direction field y´= f(x,y) provided that the local orientation of k 1 follows the direction field y´, and perfect insulation exists perpendicular to this (k 2 = 0).
The integration method, which can be derived from the hypothesis, is the topic of the dissertation by the author [1]; however, this was published without mathematical verification.Two-dimensional direction fields of the type y´= f(x,y) can describe physically different processes.In fluid mechanics, the vectorial direction field y´= v(x,y)/u(x,y) of the velocity components u and v are tangential to the streamlines.The general solution of the direction field to be integrated provides the flow field in the form of streamlines.This task is solved by default in Computational Fluid Dynamics (CFD) codes and is therefore not pursued further here.In stress analysis, the visualization of principal stress (PS) lines is not present in commercial FE programs.The visualization of stream and tensor lines is the subject of numerous investigations.A representative selection is given by [2][3][4][5][6].
Beyer [2] investigated in his dissertation principal stress trajectories arising in civil engineering and in constructions of fiber-reinforced plastics.He established an algorithm for their visualization in two and three dimensions.The procedure is similar to the method shown in Figure 1; however, improved accuracy is gained with his variable and iterative step control.
tangential to the streamlines.The general solution of the direction field to be integrated provides the flow field in the form of streamlines.This task is solved by default in Computational Fluid Dynamics (CFD) codes and is therefore not pursued further here.In stress analysis, the visualization of principal stress (PS) lines is not present in commercial FE programs.The visualization of stream and tensor lines is the subject of numerous investigations.A representative selection is given by [2][3][4][5][6].
Beyer [2] investigated in his dissertation principal stress trajectories arising in civil engineering and in constructions of fiber-reinforced plastics.He established an algorithm for their visualization in two and three dimensions.The procedure is similar to the method shown in Figure 1; however, improved accuracy is gained with his variable and iterative step control.Delmarcelle [3,4] investigated symmetric tensor fields by studying their topological structure.The basic constituents of tensor topology are the degenerate points where eigenvalues (or principal stress values) are equal to each other.Separatrices and degenerate points are the basic constituents of tensor field topology.A separatice is the tensor line (PS line) between two degenerate points.Knowledge of the topological structure of tensor fields is important, as standard integration methods need this information.
With standard integration methods the uniform distribution of tensor lines depends critically on the placement of starting points (elements), see Figure 1.Jobard's algorithm [5] for the placement of evenly-spaced PS lines (in his paper: streamlines) is based on a separating distance between adjacent PS lines to control their density.If a PS line is too close to another one, it is abruptly cut.With respect to fiber placement in FRP construction this may have adverse consequences.Trichoche [6] extends the work of Delmarcelle [3,4] to broaden the scope of topology-based vector and tensor field visualization.One focus in his work is on turbulent flows.The streamline distribution is extremely complex for turbulent flows with its local, small-scale details.Another focus is instabilities in structural problems.Visualization of stress trajectories should permit detection and identification of bifurcation to provide insight into the time-dependent structural evolution.These issues are beyond the scope of this paper.Delmarcelle [3,4] investigated symmetric tensor fields by studying their topological structure.The basic constituents of tensor topology are the degenerate points where eigenvalues (or principal stress values) are equal to each other.Separatrices and degenerate points are the basic constituents of tensor field topology.A separatice is the tensor line (PS line) between two degenerate points.Knowledge of the topological structure of tensor fields is important, as standard integration methods need this information.
With standard integration methods the uniform distribution of tensor lines depends critically on the placement of starting points (elements), see Figure 1.Jobard's algorithm [5] for the placement of evenly-spaced PS lines (in his paper: streamlines) is based on a separating distance between adjacent PS lines to control their density.If a PS line is too close to another one, it is abruptly cut.With respect to fiber placement in FRP construction this may have adverse consequences.Trichoche [6] extends the work of Delmarcelle [3,4] to broaden the scope of topology-based vector and tensor field visualization.One focus in his work is on turbulent flows.The streamline distribution is extremely complex for turbulent flows with its local, small-scale details.Another focus is instabilities in structural problems.
Visualization of stress trajectories should permit detection and identification of bifurcation to provide insight into the time-dependent structural evolution.These issues are beyond the scope of this paper.
All these papers lead to specialized algorithms, which are not public.Therefore, it is desirable to investigate the possibilities of integrating vectorial and tensorial direction fields within FE programs without programming efforts or additional software.PS lines in FRP constructions are of particular interest since fiber courses, which follow the local PS directions, make use of the superior stiffness and strength properties of the fibers [1,[7][8][9][10].
To demonstrate some disadvantages with conventional standard integration methods, Figure 1 shows an integration procedure developed by the author [1] (p. 16) applied to the bottom of a shell structure.
Figure 1b shows the PS directions of the blue principal stress trajectories on the bottom of the shell.In the middle plane and on the top of the shell there are other PS directions since membrane and bending stresses occur simultaneously.These directions stem from an eigenvalue problem of the stress tensor, whose eigenvectors (principal stress directions) remain indefinite in the sign.The PS directions can therefore switch the sign from one element to its neighboring element.The integration pattern in Figure 1c takes this problematic situation into account by replacing the direction vector with a bidirectional element.The bidirectional element intersects two opposing element sides and thereby allows the integration to be continued in the two neighboring elements with their directional elements.Each start element (marked in red) in Figure 1a provides two individual fiber courses with the procedure just described, corresponding to the two PS directions.A regular placement of the start elements in no way guarantees a uniform distribution of the PS lines.A homogenization of the fiber courses is achieved only with considerable effort.A further disadvantage of this procedure is immediately apparent; each individual trajectory (fiber) must be calculated separately.Should the designer change the fiber volume ratio for reasons of strength, then the procedure must be repeated from the beginning.This problem is substantially simplified through the hypothesis set out above, which can be modeled by means of the orthotropic heat conduction in FE programs.Since the two-dimensional direction field of the principal stresses describes a planar stress condition, shells with an arbitrary boundary can also be processed.An additional advantage is that all of the PS trajectories (fiber courses) can be calculated by means of a single thermal analysis.If the fiber volume ratio is subsequently increased, correspondingly more isotherms (fiber courses) are extracted from the existing continuous temperature field.
The recommended calculation alternative for the integration of y´= f(x,y) is based on the drastically changing properties of Fourier's law for k 1 /k 2 →∞.The following section provides a mathematical verification of the hypothesis set out above.The practicability of the integration method is based on independence from thermal boundary conditions; the isotherms always follow the direction field in nonsingular areas.However, the (weighted) distance of the isotherms from each other is influenced by the boundary conditions and can be controlled by them.In Section 3, the procedure is based on some rules and these are demonstrated practically through examples.The influence of the orthotropy ratio k 1 /k 2 , as well as the influence of singularities on the course of the isotherms, is investigated.

The Integration of Direction Fields by Means of the Orthotropic Heat Equation
The direction field is given through an ordinary differential equation of the first order in the Cartesian coordinate system (x,y): In accordance with the hypothesis in the previous section, it shall be demonstrated that orthotropic heat conduction can be used for the integration of this equation.

The Anisotropic Fourier Heat Conduction Law and the Heat Conduction Equation
Verification of the hypothesis: In the two-dimensional case, Fourier's anisotropic heat conduction law with respect to a Cartesian coordinate system is given by Equation (2), and the steady heat balance without internal heat sources by Equation (3), [11]: These equations are combined, resulting in the anisotropic heat conduction equation: Orthotropic (principal) thermal conductivities k 1 and k 2 being oriented according to the direction field of Equation ( 1) should be transformed into the Cartesian system.The transformation equations for the 2D heat conduction tensor are analogous to a 2D stress tensor, and can be illustrated using trigonometric relationships with Mohr's (Stress) Circle, as shown in Figure 2.

The Anisotropic Fourier Heat Conduction Law and the Heat Conduction Equation
Verification of the hypothesis: In the two-dimensional case, Fourier's anisotropic heat conduction law with respect to a Cartesian coordinate system is given by Equation (2), and the steady heat balance without internal heat sources by Equation (3), [11]: These equations are combined, resulting in the anisotropic heat conduction equation: Orthotropic (principal) thermal conductivities k1 and k2 being oriented according to the direction field of Equation ( 1) should be transformed into the Cartesian system.The transformation equations for the 2D heat conduction tensor are analogous to a 2D stress tensor, and can be illustrated using trigonometric relationships with Mohr's (Stress) Circle, as shown in Figure 2.
The hypothesis in Section 1 is based on extreme orthotropic conductivity ratios, k1/k2→∞ and therefore can be validated only for k2 = 0. Setting k2 = 0 in Equation ( 5) and substituting the simplified expression in Equation (2) we have: (5c) The hypothesis in Section 1 is based on extreme orthotropic conductivity ratios, k 1 /k 2 →∞ and therefore can be validated only for k 2 = 0. Setting k 2 = 0 in Equation ( 5) and substituting the simplified expression in Equation (2) we have: The heat conduction Equation ( 4) is solved in the FE programs.Substitution of Equation ( 5) into (4), together with k 2 = 0, results in: The solution T(x,y) of Equation ( 7) can be visualized by isotherms, which represent a direction field y T ´.An isotherm is characterized by a constant temperature C. Its equation in implicit form is given by F(x,y) = T(x,y) − C = 0.A directional element of this curve has the slope: The question arises as to which properties Equation ( 6) shows if the direction fields from Equations ( 1) and ( 8) are equated: Use of "f " from Equation ( 9) in the numerators of Equation ( 6) results in: The first and second term for q x and q y cancel out, and the right hand sides in Equations (10a) and (10b) become zero.Back-substitution of (∂T/∂x)/(∂T/∂y) = −f into Equation (10) results in: For this extreme case (k 2 = 0), both heat flux components [q x , q y ] = q are identical to zero.Therefore, the heat flux components q 1 and q 2 parallel to the orthotropic (principal) thermal conductivities (k 1 , k 2 ) in Figure 2a are also zero.The resulting isotherms do not contradict their definition; the heat flux along an isotherm is zero and, perpendicular to this, the heat flux is also zero due to k 2 = 0, and the equating of Equations ( 1) and ( 8) was justified: If the heat fluxes are zero for any direction field y´= f(x,y), then their divergence div(q) in the form of Equation ( 3) or ( 4) is also zero.The latter is the heat conduction equation, which is solved in FE programs for steady-state problems without internal heat sources.From Fourier's law with the unknowns T and q, a directly solvable equation has been formed for T alone provided that k 2 = 0. Whether the orthotropic heat conduction problem in Equation ( 11) is solved, which is not generally available as a calculation option in FE programs, or Equation ( 7) is solved by means of an FE program, the same isotherms result and are tangential to the given direction field y´= f(x,y).Whenever k 2 > 0, the proof does not work.

The Relationship between the Isotherms of Extremely Orthotropic Heat Conduction and the Characteristics of a Partial Linear Differential Equation of the First Order
A linear partial differential equation of the first order is: and is discussed in detail in the literature [12][13][14].The required function T(x,y) depicts a surface in three-dimensional space.The theory of the "characteristic method" [12][13][14] for solving Equation (12) implies that its characteristics are the contour lines of T(x,y) and can be calculated with the ordinary differential equation: with A = 0.The similarity of Equations ( 11) and ( 12) is obvious for A = 1: Along the characteristics, the solution variable T is constant.Therefore, they can be interpreted as isotherms.Equation ( 14) physically describes the convective transport of the T variable.Equation ( 12) is valid for a partial differential equation with two independent variables.If x and y are used, then there is a time-independent situation in 2D space.If x is used and t instead of y, then there is a time-dependent situation in 1D space, the latter associates the transport character of Equation (12).A simple transport equation reads [12]: with time t and the convection velocity u = dx/dt.The variable T represents temperature here; it could also stand for mass, momentum, or energy.In order to demonstrate the particular properties of Fourier heat conduction (k 2 = 0) in a simple example, f = y' = 1 is substituted in Equation ( 11) and u = 1 in Equation (15).Therefore, both equations are formally identical.The location variable y in Equation ( 11) corresponds to the time variable t in Equation ( 15), y and t are interchangeable.By simple integration of Equation ( 13), dy/dx = f = 1, the isotherms are represented by the straight lines y = x + C. Figure 3 shows the corresponding solution for Equations ( 11) and ( 15), which was obtained with the FE program ABAQUS (Dassault Systèmes, Vélizy-Villacoublay, France).The initial condition for T(t = 0, x 0 ≤ x ≤ x 1 ) and T(y = 0, x 0 ≤ x ≤ x 1 ) is given by a parabolic temperature distribution.The orientation of k 1 is given by the vector (1, 1) in the entire solution area.Extreme orthotropy is achieved by zero-setting of k 2 whereby k 2 is perpendicular to k 1 .
The isotherms in Figure 3 depict both the general solution for the two-dimensional Equation (11) in the x-y plane as well as the time-dependent solution for the one-dimensional Equation (15) in the x-t plane.The isotherms follow the direction field y´= 1, not only inside the rectangular area but also on the boundary.
The hypothesis in Section 1 restricts usage to steady state problems.If the FE user investigates the behavior of the (extreme) orthotropic heat conduction beyond this restriction also in the time domain (with three independent variables: x, y, t), then the isotherms are no longer parallel to the specified direction field for finite times.The parabolic start temperature distribution propagates, corresponding to the thermal diffusivity a = k 1 /(cρ) with finite velocity (c = specific heat, ρ = density).As can be seen from Figure 4, the illustratively plotted 0.01•T max -isotherms at time t 1 < t 2 < t 3 are not parallel to the direction field y´= 1.This only applies for the steady-state case when t→∞.

Temperature Boundary Conditions for Extremely Orthotropic Heat Conduction
The Equations ( 4) and ( 7) are solved in FE programs.These are second order and require Dirichlet or Neumann boundary conditions.Non-specified conditions at the boundary include ∂T/∂n = q = 0 (n = normal vector perpendicular to the boundary).In case of extreme orthotropy when k 1 > 0 and k 2 = 0, the latter condition is always fulfilled according to (11) both internally as well as on the boundary; therefore, in principle, no further conditions need to be specified.However, without specification of a reference temperature, the FE solution is singular.If a reference temperature T ref is specified at an arbitrary node, we only obtain the trivial solution T(x,y) = T ref .Therefore, at least two different temperatures at two different nodes need to be specified in order to obtain a non-trivial solution.All temperature boundary conditions, which are specified at two or more nodes, provide isotherms that are tangential to the given direction field y´= f(x,y).The difference between a two-node and a multiple-node boundary condition manifests itself through different weightings with respect to the location of the isotherms.In order to obtain appropriate solutions, the direction field must be viewed in connection with the physical task, see Figure 5.This figure shows a disc (half-model) with an elliptical hole under pure bending.

Temperature Boundary Conditions for Extremely Orthotropic Heat Conduction
The Equations ( 4) and ( 7) are solved in FE programs.These are second order and require Dirichlet or Neumann boundary conditions.Non-specified conditions at the boundary include ∂T/∂n = q = 0 (n = normal vector perpendicular to the boundary).In case of extreme orthotropy when k1 > 0 and k2 = 0, the latter condition is always fulfilled according to Equation (11) both internally as well as on the boundary; therefore, in principle, no further conditions need to be specified.However, without specification of a reference temperature, the FE solution is singular.If a reference temperature Tref is specified at an arbitrary node, we only obtain the trivial solution T(x,y) = Tref.Therefore, at least two different temperatures at two different nodes need to be specified in order to obtain a non-trivial solution.All temperature boundary conditions, which are specified at two or more nodes, provide isotherms that are tangential to the given direction field y´ = f(x,y).The difference between a two-node and a multiple-node boundary condition manifests itself through different weightings with respect to the location of the isotherms.In order to obtain appropriate solutions, the direction field must be viewed in connection with the physical task, see Figure 5.This figure shows a disc (half-model) with an elliptical hole under pure bending.From a practical point of view, FRP constructions should have a uniform fiber placement.This is achieved with an orthotropic heat conduction analysis, specifying a linear temperature distribution for y > 0 along the right boundary (Figure 5a).The isotherms, as representatives of the σ1 trajectories, are then distributed equidistantly.For y < 0, T = 0 is defined.The σ2 trajectories are calculated in a further orthotropic heat conduction analysis.These require a linear temperature distribution for y < 0 and T = 0 for y > 0. Both isotherm images are superimposed and result in the two fiber layers in Figure 5a.
Should the objective be concentrated on load-related fiber placement, then a linear weighted distribution of isotherms is achieved by defining a parabolic temperature distribution for y > 0 (σ1 trajectories) and y < 0 (σ2 trajectories) along the right boundary; see Figure 5b.
The two-node temperature boundary condition is likewise possible.For example, to determine the σ1 trajectories, specification of T(L,h/2) = Tmax and T(L,0) = 0 provides similar results, though the isotherms are not distributed equidistantly on the right boundary; see also Section 3.1.From a practical point of view, FRP constructions should have a uniform fiber placement.This is achieved with an orthotropic heat conduction analysis, specifying a linear temperature distribution for y > 0 along the right boundary (Figure 5a).The isotherms, as representatives of the σ 1 trajectories, are then distributed equidistantly.For y < 0, T = 0 is defined.The σ 2 trajectories are calculated in a further orthotropic heat conduction analysis.These require a linear temperature distribution for y < 0 and T = 0 for y > 0. Both isotherm images are superimposed and result in the two fiber layers in Figure 5a.
Should the objective be concentrated on load-related fiber placement, then a linear weighted distribution of isotherms is achieved by defining a parabolic temperature distribution for y > 0 (σ 1 trajectories) and y < 0 (σ 2 trajectories) along the right boundary; see Figure 5b.
The two-node temperature boundary condition is likewise possible.For example, to determine the σ 1 trajectories, specification of T(L,h/2) = T max and T(L,0) = 0 provides similar results, though the isotherms are not distributed equidistantly on the right boundary; see also Section 3.1.
In Section 2.1 it was shown that in cases of extreme orthotropy (k 1 > 0, k 2 = 0), the differential Equation ( 7) of the second order and the Fourier heat conduction Equation (11) of the first order physically describe the same situation.In Section 2.2, it was shown that the Fourier heat conduction Equation ( 11) can be traced back to the general linear differential Equation (12).The integration of the corresponding characteristic Equation ( 13) provides characteristics identical to the isotherms in the Fourier heat conduction Equation (11).The theory of the characteristics also includes the treatment of boundary conditions [14], the rules of which can be applied completely to the Fourier heat conduction Equation ( 14). Figure 6 shows the most important rule.In Section 2.1 it was shown that in cases of extreme orthotropy (k1 > 0, k2 = 0), the differential Equation ( 7) of the second order and the Fourier heat conduction Equation (11) of the first order physically describe the same situation.In Section 2.2, it was shown that the Fourier heat conduction Equation ( 11) can be traced back to the general linear differential Equation (12).The integration of the corresponding characteristic Equation ( 13) provides characteristics identical to the isotherms in the Fourier heat conduction Equation (11).The theory of the characteristics also includes the treatment of boundary conditions [14], the rules of which can be applied completely to the Fourier heat conduction Equation ( 14). Figure 6 shows the most important rule.(11).The specification of arbitrary boundary conditions along P2P3 is not permitted since this curve is a characteristic or isoline (y = const.).
According to this, Equation ( 12) is integrated through the specification of the solution variables T along a curve Γ.However, only the specification on the P1P2 curve is permitted.Since a solution variable is constant on a characteristic, arbitrary specifications along P2P3 are not permitted.Furthermore, arbitrary specifications along P3P4 lead to contradictory definitions with the specifications along P1P2.The curve Γ must not necessarily lie on the exterior boundary, but can also cross the integration area.

Procedure and Application Examples
The direction field y´ = f(x,y) to be integrated very often originates from FE analysis in the form of PS directions.The procedure for their integration shall be determined with the help of some rules: 1. Calculation of the PS directions in all elements (disc, plate, or shell).If this is not provided by the respective FE program used, they can be calculated via the stress components [σxx σyy σxy].
For example, y´ = dy/dx = tanα = σxy/(σxx − σ2); see also Figure 2c.Plates and shells have variable stresses across the thicknesses.The direction field therefore must be evaluated for every "thickness integration point".2. Replacement of the structural elements with thermal elements.(11).The specification of arbitrary boundary conditions along P2P3 is not permitted since this curve is a characteristic or isoline (y = const.).
According to this, Equation ( 12) is integrated through the specification of the solution variables T along a curve Γ.However, only the specification on the P1P2 curve is permitted.Since a solution variable is constant on a characteristic, arbitrary specifications along P2P3 are not permitted.Furthermore, arbitrary specifications along P3P4 lead to contradictory definitions with the specifications along P1P2.The curve Γ must not necessarily lie on the exterior boundary, but can also cross the integration area.

Procedure and Application Examples
The direction field y´= f(x,y) to be integrated very often originates from FE analysis in the form of PS directions.The procedure for their integration shall be determined with the help of some rules: 1.
Calculation of the PS directions in all elements (disc, plate, or shell).If this is not provided by the respective FE program used, they can be calculated via the stress components [σ xx σ yy σ xy ].For example, y´= dy/dx = tanα = σ xy /(σ xx − σ 2 ); see also Figure 2c.Plates and shells have variable stresses across the thicknesses.The direction field therefore must be evaluated for every "thickness integration point".

2.
Replacement of the structural elements with thermal elements.

3.
Transmission of the PS directions (Step 1) to local systems, depending on the FE program used, (ABAQUS: *Orientation).
Definition of thermal boundary conditions according to Section 2.3.

6.
The ratios k 1 /k 2 > 10 4 and k 2 /k 1 > 10 4 provide the PS trajectories for the first and for the second principal stress, respectively.
It should be emphasized that the integration of the PS directions for the calculation of the real fiber course is a separate task, independent of the structural analysis.The PS direction field of the structural calculation is adopted as material orientation for the thermal conductivities in the thermal analysis.The boundary conditions are independent in both analyses and have nothing in common with each other.

Influence of the Orthotropy Ratio k 1 /k 2 on the Integration Accuracy
The influence of the orthotropy ratio k 1 /k 2 on the integration accuracy is demonstrated by the following example.Figure 7 shows a perforated disc under axial tension.Due to symmetry, a quarter model is sufficient (U x = 0, U y = 0).The supplementary material accompanying this paper includes the ABAQUS input for this example.
It should be emphasized that the integration of the PS directions for the calculation of the real fiber course is a separate task, independent of the structural analysis.The PS direction field of the structural calculation is adopted as material orientation for the thermal conductivities in the thermal analysis.The boundary conditions are independent in both analyses and have nothing in common with each other.

Influence of the Orthotropy Ratio k1/k2 on the Integration Accuracy
The influence of the orthotropy ratio k1/k2 on the integration accuracy is demonstrated by the following example.Figure 7 shows a perforated disc under axial tension.Due to symmetry, a quarter model is sufficient (Ux = 0, Uy = 0).The supplementary material accompanying this paper includes the ABAQUS input for this example.The task is to place the fibers in the direction of the largest principal stress for an FRP construction.The PS directions portray the direction field y´ = f(x,y), which shall be integrated.The elements of the mechanical model are replaced with heat conduction elements whose heat conduction k1 is oriented parallel to the previously calculated PS directions of σ1.The k2 direction is orthogonal to k1.Now the ratio k1/k2 is successively increased, thus the isotherms get more and more tangential to the PS direction of σ1.The influence of the two-node temperature boundary condition at points P2 and P3 (Figure 7d) on the course of isotherms is no longer perceptible for k1/k2 > 10 4 .Under such simple boundary conditions, it must be ensured that their positions are meaningful.If both temperatures are determined at points P1 and P2, then, practically, two different temperatures are defined on one isotherm.
The relatively uniform distribution of the PS trajectories is remarkable, although only a very simple temperature boundary condition was selected.Should the desired outcome be a distribution of isotherms exactly equidistant between P2 and P3, the temperature along this line must be specified as linearly increasing.Since the isotherms represent the fiber course, the constant line load along P3P2 The task is to place the fibers in the direction of the largest principal stress for an FRP construction.The PS directions portray the direction field y´= f(x,y), which shall be integrated.The elements of the mechanical model are replaced with heat conduction elements whose heat conduction k 1 is oriented parallel to the previously calculated PS directions of σ 1 .The k 2 direction is orthogonal to k 1 .Now the ratio k 1 /k 2 is successively increased, thus the isotherms get more and more tangential to the PS direction of σ 1 .The influence of the two-node temperature boundary condition at points P2 and P3 (Figure 7d) on the course of isotherms is no longer perceptible for k 1 /k 2 > 10 4 .Under such simple boundary conditions, it must be ensured that their positions are meaningful.If both temperatures are determined at points P1 and P2, then, practically, two different temperatures are defined on one isotherm.
The relatively uniform distribution of the PS trajectories is remarkable, although only a very simple temperature boundary condition was selected.Should the desired outcome be a distribution of isotherms exactly equidistant between P2 and P3, the temperature along this line must be specified as linearly increasing.Since the isotherms represent the fiber course, the constant line load along P3P2 is uniformly carried by the fibers.The maximum load in the perforated disc at points A and B in Figure 7e is approximately three times the applied tension, which is well reflected by the fiber course.
At this point, the transport character of the Fourier heat conduction equation in Figure 7d and Equations ( 14) and ( 15) shall be emphasized again.The specification of a linear temperature distribution on the boundary P3P2 is transported along its characteristics without diffusion effects, and these characteristics follow the user defined direction field y´= f(x,y).

The Influence of Singularities on the Course of the Isotherms
In cases of stress singularities, the following distinctions are made: An isotropic (or neutral) point with an undefined PS direction occurs if the two principal stresses σ 1 and σ 2 are identical.If both principal stresses are zero, then a classic singularity is present.Points with concentrated loads are neither one nor the other.Their isoclines are particularly concentrated at the point of force transmission.What does the course of stress trajectories (represented by isotherms) look like in the vicinity of these points?This can be studied using an example of a circular ring under diametrically opposed single forces.Figure 8a shows the course of the isoclines and stress trajectories that were determined photoelastically by Frocht [15]. is uniformly carried by the fibers.The maximum load in the perforated disc at points A and B in Figure 7e is approximately three times the applied tension, which is well reflected by the fiber course.At this point, the transport character of the Fourier heat conduction equation in Figure 7d and Equations ( 14) and ( 15) shall be emphasized again.The specification of a linear temperature distribution on the boundary P3P2 is transported along its characteristics without diffusion effects, and these characteristics follow the user defined direction field y´ = f(x,y).

The Influence of Singularities on the Course of the Isotherms
In cases of stress singularities, the following distinctions are made: An isotropic (or neutral) point with an undefined PS direction occurs if the two principal stresses σ1 and σ2 are identical.If both principal stresses are zero, then a classic singularity is present.Points with concentrated loads are neither one nor the other.Their isoclines are particularly concentrated at the point of force transmission.What does the course of stress trajectories (represented by isotherms) look like in the vicinity of these points?This can be studied using an example of a circular ring under diametrically opposed single forces.Figure 8a shows the course of the isoclines and stress trajectories that were determined photoelastically by Frocht [15].The experimental results were verified through an FE stress analysis (ABAQUS) followed by an orthotropic heat conduction analysis in order to determine the trajectories, as shown in Figure 8b.Modeling a quarter of the disc is sufficient due to the symmetry.The stress calculation provides the PS directions.In the subsequent heat conduction calculation, these PS directions serve as local systems for the thermal conductivities k1 and k2.Two of these systems are illustratively plotted in the fourth quadrant of Figure 8b.The system k1 = 1, k2 = 0 provides the isotherms as σ1 trajectories.With k1 = 0, k2 = 1, the σ2 trajectories are determined.The calculation of the isoclines is helpful in order to localize the singular points A, D, E, H, and K (Figure 8b, left).The isoclines result from Equation ( 16), i.e., the contour lines of the stress expression on the right side must be visualized; see also Figure 2c.
In the trajectories image (right), the singular points appear with the designations B, C, F, G, and L, and are either of the "non-interlocking" type (B, C, F, G) or "interlocking" type (L).This evaluation The experimental results were verified through an FE stress analysis (ABAQUS) followed by an orthotropic heat conduction analysis in order to determine the trajectories, as shown in Figure 8b.Modeling a quarter of the disc is sufficient due to the symmetry.The stress calculation provides the PS directions.In the subsequent heat conduction calculation, these PS directions serve as local systems for the thermal conductivities k 1 and k 2 .Two of these systems are illustratively plotted in the fourth quadrant of Figure 8b.The system k 1 = 1, k 2 = 0 provides the isotherms as σ 1 trajectories.With k 1 = 0, k 2 = 1, the σ 2 trajectories are determined.The calculation of the isoclines is helpful in order to localize the singular points A, D, E, H, and K (Figure 8b, left).The isoclines result from Equation ( 16), i.e., the contour lines of the stress expression on the right side must be visualized; see also Figure 2c.
In the trajectories image (right), the singular points appear with the designations B, C, F, G, and L, and are either of the "non-interlocking" type (B, C, F, G) or "interlocking" type (L).This evaluation is useful since a uniform distribution of the isotherms (fibers) is desirable, which is only achievable with knowledge of the singularities.This is shown in Figure 9a for the σ 2 trajectories (with σ 2 < σ 1 ).A uniform distribution is achieved through a linearly increasing temperature along RF, with continuation along JL. with knowledge of the singularities.This is shown in Figure 9a for the σ2 trajectories (with σ2 < σ1).A uniform distribution is achieved through a linearly increasing temperature along RF, with continuation along JL.Specifying a linear increase in temperature from R to J is incorrect since the σ2 trajectory along RF turns into the σ1 trajectory at point F along FJ.The FJ boundary segment itself is an isotherm and no variable temperature specification is allowed there.
The existence of a singular point on a load-free boundary of a 2D structure always results in this change of trajectories.This property is responsible for the fact that efforts need to be made regarding the specification of suitable temperature boundary conditions in order to obtain uniformly distributed isotherms.Should the temperature be defined only along RF, then only a part of the structure will show isotherms.A further boundary segment with temperature boundary conditions must then be applied in order to capture the rest of the structure.In Figure 9a, this is segment JL.The procedure is the same when looking for σ1 trajectories.The temperature is defined as linearly increasing along PB with continuation along LN.
The aim of calculating strictly uniformly distributed isotherms (trajectories) can only be achieved locally.Figure 9a shows uniformly distributed σ1 trajectories along PB; however, these end up unevenly distributed on the PR boundary.
On the perforated disc (Figure 7), the trajectory change mentioned above does not occur along the right boundary of P2P3.This boundary is free from singularities.The specification of two temperatures at points P2 and P3 is sufficient.We can generally use this simple two-node boundary condition (Section 2.3) for all 2D discs and 3D shells and get-depending on the number of singularities-correspondingly weighted isotherm distributions.Figure 9b shows this simple procedure for the circular ring with its singularities.The isotherms shown there are less uniformly distributed than in Figure 9a, but this is adequate for a first overview.
The circular ring problem (Figure 8) contains the most important forms of singularities: Isotropic points of the "interlocking" and "non-interlocking" type, as well as concentrated point loads.Errors in representation of isotherms are inevitable since they do not fulfil y´ = f(x,y) in the near vicinity of the singularity.In fact, the isotherms in close proximity there evade in a sideways direction Specifying a linear increase in temperature from R to J is incorrect since the σ 2 trajectory along RF turns into the σ 1 trajectory at point F along FJ.The FJ boundary segment itself is an isotherm and no variable temperature specification is allowed there.
The existence of a singular point on a load-free boundary of a 2D structure always results in this change of trajectories.This property is responsible for the fact that efforts need to be made regarding the specification of suitable temperature boundary conditions in order to obtain uniformly distributed isotherms.Should the temperature be defined only along RF, then only a part of the structure will show isotherms.A further boundary segment with temperature boundary conditions must then be applied in order to capture the rest of the structure.In Figure 9a, this is segment JL.The procedure is the same when looking for σ 1 trajectories.The temperature is defined as linearly increasing along PB with continuation along LN.
The aim of calculating strictly uniformly distributed isotherms (trajectories) can only be achieved locally.Figure 9a shows uniformly distributed σ 1 trajectories along PB; however, these end up unevenly distributed on the PR boundary.
On the perforated disc (Figure 7), the trajectory change mentioned above does not occur along the right boundary of P2P3.This boundary is free from singularities.The specification of two temperatures at points P2 and P3 is sufficient.We can generally use this simple two-node boundary condition (Section 2.3) for all 2D discs and 3D shells and get-depending on the number of singularities-correspondingly weighted isotherm Figure 9b shows this simple procedure for the circular ring with its singularities.The isotherms shown there are less uniformly distributed than in Figure 9a, but this is adequate for a first overview.
The circular ring problem (Figure 8) contains the most important forms of singularities: Isotropic points of the "interlocking" and "non-interlocking" type, as well as concentrated point loads.Errors in representation of isotherms are inevitable since they do not fulfil y´= f(x,y) in the near vicinity of the singularity.In fact, the isotherms in close proximity there evade in a sideways direction depending on the model mesh size.The model in Figure 8b, in comparison to Figure 8a, is clearly fine enough to minimize this error.The element size is virtually the same all over; 60 elements were used (with four nodes respectively) along JN.Through further mesh refinement in the area of singularities, the error can be kept arbitrary small.

The Principal Stress Trajectories in Shell Structures
A conventional method for the calculation of principal stress trajectories in a shell was shown in Figure 1.Disadvantages resulted from the uneven fiber placement, which originates from the separate calculation of each individual fiber.In integrating the direction field, it is not practically possible to select optimal start positions (start elements) beforehand.
In contrast to this is the method of integration by means of an orthotropic heat conduction calculation.Even with a two-node temperature boundary condition with Q 0 : T = 0 • C, Q 1 : T = 1 • C in Figure 10a, a better fiber distribution can be attained than is depicted in Figure 1a.This temperature distribution, shown in Figure 10a, provides the approximate position of three local temperature extremes at the positions P 0 , P 1 and P 2 .

The Principal Stress Trajectories in Shell Structures
A conventional method for the calculation of principal stress trajectories in a shell was shown in Figure 1.Disadvantages resulted from the uneven fiber placement, which originates from the separate calculation of each individual fiber.In integrating the direction field, it is not practically possible to select optimal start positions (start elements) beforehand.
In contrast to this is the method of integration by means of an orthotropic heat conduction calculation.Even with a two-node temperature boundary condition with Q0: T = 0 °C, Q1: T = 1 °C in Figure 10a, a better fiber distribution can be attained than is depicted in Figure 1a.This temperature distribution, shown in Figure 10a, provides the approximate position of three local temperature extremes at the positions P0, P1 and P2.By repeating the heat conduction analysis with a three-node temperature boundary condition (P0: T = 0 °C, P1: T = 1 °C, P2: T = 2 °C, Figure 10b), the uniformity of the isotherms is improved.A specification of the linear varying temperatures along boundary segments analogous to the previous section would yield further improvement.However, in light of the result in Figure 10b, it does not appear to be crucial.
The right half of Figure 8b highlights another aspect: The continuous temperature field is infinitely dense with respect to the isotherms, so an arbitrary number of contours can be visualized, avoiding calculation of trajectories one after another (Figure 1).

Optimization of Fiber Placement in FRP Constructions
In order to demonstrate a useful application beyond the visualization of tensor lines, the fiber By repeating the heat conduction analysis with a three-node temperature boundary condition (P 0 : T = 0 • C, P 1 : T = 1 • C, P 2 : T = 2 • C, Figure 10b), the uniformity of the isotherms is improved.A specification of the linear varying temperatures along boundary segments analogous to the previous section would yield further improvement.However, in light of the result in Figure 10b, it does not appear to be crucial.right half of Figure 8b highlights another aspect: The continuous temperature field is infinitely dense with respect to the isotherms, so an arbitrary number of contours can be visualized, avoiding calculation of trajectories one after another (Figure 1).

Optimization of Fiber Placement in FRP Constructions
In order to demonstrate a useful application beyond the visualization of tensor lines, the fiber placement in FRP constructions is highlighted.Fiber courses, which follow local PS directions, use the superior stiffness and strength properties of the fibers.In plane stress structures the fiber courses due to σ 1 and σ 2 are perpendicular to each other, forming a curvilinear cross-ply (CCP) laminate.With regard to optimality, attention must be paid to the algebraic signs of principal stresses.
Fiber placement in the direction of principal stresses is optimal according to net theory [16], and this is also valid for CCP-laminates if ignoring the support of the matrix.To be on the safe side, the carrying capacity of the construction can be calculated by using a fiber failure criterion.A design based on an inter-fiber failure must additionally take into account the contribution of the matrix.In this case CCP-laminates are optimal only in areas with principal stresses showing different algebraic signs.However, if the algebraic signs are the same, then the curvilinear balanced angle-ply (CBAP) laminate is the better choice according to classical lamination theory [17].This design differs from the (right-angled) CCP-laminate through a correction angle ±β regarding the σ 1 direction.The direction field y´for the largest principal stress σ 1 is then represented by two direction fields: y´+ β and y´− β, as shown in Figure 11.signs.However, if the algebraic signs are the same, then the curvilinear balanced angle-ply (CBAP) laminate is the better choice according to classical lamination theory [17].This design differs from the (right-angled) CCP-laminate through a correction angle ±β regarding the σ1 direction.The direction field y´ for the largest principal stress σ1 is then represented by two direction fields: y´ + β and y´ − β, as shown in Figure 11.For example, the perforated disc (Figure 7e) contains the areas σ1·σ2 < 0 (CCP-laminate is optimal) and σ1·σ2 > 0 (CBAP-laminate is optimal).In the latter case the modified direction fields can also be integrated without restriction by means of the orthotropic heat conduction analysis.The correction angle β is dependent on the material and the principal stress ratio σ1/σ2 [17].It should be noted that the orthotropic heat conduction is suitable for the integration of arbitrary direction fields even if they are modified for reasons determined by the engineer.Reference is made to Moldenhauer [1] (pp. 54-60] and [18] for further details and consideration of particular aspects of load-related optimization (layer thickness distribution, loading condition range, and nonlinearity of stress distribution over the shell thickness).

Discussion and Conclusions
The method presented in this paper is suitable for integrating arbitrary vectorial and tensorial 2D direction fields by means of commercial FE programs, provided that the FE programs have For example, the perforated disc (Figure 7e) contains the areas σ 1 •σ 2 < 0 (CCP-laminate is optimal) and σ 1 •σ 2 > 0 (CBAP-laminate is optimal).In the latter case the modified direction fields can also be integrated without restriction by means of the orthotropic heat conduction analysis.The correction angle β is dependent on the material and the principal stress ratio σ 1 /σ 2 [17].It should be noted that the orthotropic heat conduction is suitable for the integration of arbitrary direction fields even if they are modified for reasons determined by the engineer.Reference is made to Moldenhauer [1] (pp. 54-60] and [18] for further details and consideration of particular aspects of load-related optimization (layer thickness distribution, loading condition and nonlinearity of stress distribution over the shell thickness).

Discussion and Conclusions
The method presented in this paper is suitable for integrating arbitrary vectorial and tensorial 2D direction fields by means of commercial FE programs, provided that the FE programs have options for the analysis of orthotropic heat conduction.The calculated isotherms represent the integral curves of these direction fields inside and on the edge of the 2D structure (disc or shell).The procedure is based on some rules (Section 3), which concern only the standard input of the FE program.Additional programming is not required.Two temperatures at two different positions shall be specified as a minimum for thermal boundary conditions.The isotherms arising from this can still run in a non-uniform way.Uniformity is achieved while linearly varying temperatures along a boundary segment are specified (Figure 5a).Should the direction field or the differential equation y´= f(x,y) exhibit singularities on the boundary, then the definition of the boundary conditions must be carefully adjusted (Figure 9).In this case, a preliminary visualization of the isoclines, Equation ( 16), and/or an analysis with the two-node temperature boundary condition is helpful.The precision of the integration depends on the element size of the FE model.
In comparison with the traditional integration methods, significant improvements arise with the new procedure:

•
The integration of the direction fields can take place with standard FE programs with options for an orthotropic heat conduction analysis.No additional programming is necessary.

•
Any number of isotherms can be extracted from the integrated continuous temperature field.

•
Quick and accurate results are obtained with a simple two-node temperature boundary condition; however, the distribution of isotherms may be non-uniform if singular points are present in the direction field.

•
As the computed temperature field is continuous, valuable additional information can be extracted from this field.For example, the density of the isotherms (tensor lines) can be computed by displaying the temperature gradients, see [1] (p. 96).

•
Applying the proposed integration method to stress tensors it should be noted that principal stress direction fields can originate from linear or nonlinear analyses, see [1] (p. 74).If principal shear stress (or strain) directions are evaluated in the plastic range then the trajectories can be regarded as slip-lines, see [1] (pp. 41-44).

•
Civil engineering: Aligning reinforcement in concrete structures parallel to PS directions is a meaningful tool to effectively increase the low tension strength of concrete, see also [1] (p. 78).
Mathematical aspects were the focus of this paper.An important application seems to be on FRP constructions.Aligning fiber courses along PS lines make the most of the superior stiffness and strength properties of fibers.Calculation of these fiber courses by integrating the PS directions followed by the evaluation of this optimized FRP design can be done with one and the same FE program.To evaluate and optimize an FRP construction based on PS lines, some practical aspects should be highlighted.
As with all FRP problems, the optimization success is greatest when a special load case is investigated.In reality, loads and boundary conditions change frequently, and the principal stress directions change accordingly.However, the associated degradation of the optimization method affects alternative optimization strategies as well.If one load case dominates other load cases, the optimization effort can be worthwhile nevertheless.A CCP-or CBAP-laminate can be combined with textile preforms and can be stitched onto standard laminates.This technique is known as tailored fiber placement (TFP).The layer with the fiber course aligned to the trajectories covers the critical load case, and the standard laminate the other load cases.In 2D, the TFP stitching technique has been used for a long time [10].In 3D, textile preforming or tape laying are possible options.

Figure 1 .
Figure 1.Shell under lateral pressure p = 0.15 MPa, bedded elastically (Kf = 1000 N/mm); typical dimensions Δx, Δy, and Δz, modulus E = 50 GPa, thickness t = 2 mm.(a) stress trajectories, standard integration begins in the red start elements; (b) PS directions from a linear static analysis; (c) detail, standard integration begins in element A. Improved accuracy is gained with mesh refinement.

Figure 1 .
Figure 1.Shell under lateral pressure p = 0.15 MPa, bedded elastically (K f = 1000 N/mm); typical dimensions ∆x, ∆y, and ∆z, modulus E = 50 GPa, thickness t = 2 mm.(a) stress trajectories, standard integration begins in the red start elements; (b) PS directions from a linear static analysis; (c) detail, standard integration begins in element A. Improved accuracy is gained with mesh refinement.

Figure 3 .
Figure 3.The Fourier heat conduction equation with k1/k2→∞ is hyperbolic and has transport character.The diagram is valid for the x-y plane as well as for the x-t plane.As an example, the direction field y´ = f(x,y) = 1 is integrated with the orthotropic heat conduction (k2 = 0).

Figure 4 .
Figure 4.The steady-state problem from Figure 3 analyzed in the time domain.The isotherms are parallel to the given direction field y´ = 1 only for the final steady state when t∞.

Figure 3 .
Figure 3.The Fourier heat conduction equation with k 1 /k 2 →∞ is hyperbolic and has transport character.The diagram is valid for the x-y plane as well as for the x-t plane.As an example, the direction field y´= f(x,y) = 1 is integrated with the orthotropic heat conduction (k 2 = 0).

Figure 3 .
Figure 3.The Fourier heat conduction equation with k1/k2→∞ is hyperbolic and has transport character.The diagram is valid for the x-y plane as well as for the x-t plane.As an example, the direction field y´ = f(x,y) = 1 is integrated with the orthotropic heat conduction (k2 = 0).

Figure 4 .
Figure 4.The steady-state problem from Figure 3 analyzed in the time domain.The isotherms are parallel to the given direction field y´ = 1 only for the final steady state when t∞.

Figure 4 .
Figure 4.The steady-state problem from Figure 3 analyzed in the time domain.The isotherms are parallel to the given direction field y´= 1 only for the final steady state when t→∞.

Figure 5 .
Figure 5. FRP disc with hole under bending: (a) Linearly increasing temperature boundary conditions for σ1 and σ2 trajectories yield a uniform distribution of fibers (isotherms); (b) parabolic temperature boundary conditions yield a weighted distribution of fibers, which carries a bending moment more effectively.

Figure 5 .
Figure 5. FRP disc with hole under bending: (a) Linearly increasing temperature boundary conditions for σ 1 and σ 2 trajectories yield a uniform distribution of fibers (isotherms); (b) parabolic temperature boundary conditions yield a weighted distribution of fibers, which carries a bending moment more effectively.

Figure 6 .
Figure 6.The solution curves T(x,y) of the direction field y´ = f(x,y) are simultaneously the characteristics of the Fourier heat conduction Equation(11).The specification of arbitrary boundary conditions along P2P3 is not permitted since this curve is a characteristic or isoline (y = const.).

Figure 6 .
Figure 6.The solution curves T(x,y) of the direction field y´= f(x,y) are simultaneously the characteristics of the Fourier heat conduction Equation(11).The specification of arbitrary boundary conditions along P2P3 is not permitted since this curve is a characteristic or isoline (y = const.).

Figure 9 .
Figure 9. Circular ring analogous to Figure 8: (a) Determination of the σ2 trajectory through specification of a linearly increasing temperature along RF and JL; (b) Simplified analysis with a twonode temperature boundary condition; T0 and TL for the σ1 trajectories, T0 and TL for the σ2 trajectories.

Figure 9 .
Figure 9. Circular ring analogous to Figure 8: (a) Determination of the σ 2 trajectory through specification of a linearly increasing temperature along RF and JL; (b) Simplified analysis with a two-node temperature boundary condition; T 0 and T L for the σ 1 trajectories, T 0 and T L for the σ 2 trajectories.

Figure 10 .
Figure 10.Shell structure analogous to Figure 1: (a) Isotherms with a two-node temperature boundary condition at Q0 and Q1; (b) Isotherms with a three-node temperature boundary condition at P0, P1, P2, right: Temperature field for extraction of any number of isotherms (fibers).

Figure 10 .
Figure 10.Shell structure analogous to Figure 1: (a) Isotherms with a two-node temperature boundary condition at Q 0 and Q 1 ; (b) Isotherms with a three-node temperature boundary condition at P 0, P 1 , P 2 , right: Temperature field for extraction of any number of isotherms (fibers).

Figure 11 .
Figure 11.Perforated disc (detail) similar to Figure 7e.In areas with the same principal stress algebraic sign, the CBAP-laminate is optimal: (a) Correction angle ±ß regarding the σ1 direction; (b) Integration of the modified direction field.

Figure 11 .
Figure 11.Perforated disc (detail) similar to Figure 7e.In areas with the same principal stress algebraic sign, the CBAP-laminate is optimal: (a) Correction angle ±ß regarding the σ 1 direction; (b) Integration of the modified direction field.