Non-incremental strategies for simulating thermomechanical models with uncertainty

: Description of complex materials involves numerous computational challenges. Thus, the analysis of transient models needs for intensive computation. In some of our former works we proposed a technique based on the separated representation of the unknown field able to circumvent the curse of dimensionality in the treatment of highly multidimensional models. In this work, we are addressing the challenge related to the transient behaviour. For this purpose we propose a separated representation of transient models leading to a non-incremental strategy, allowing to impressive CPU time savings. The use of separated representations allows computing the solution of transient parametric models (the parameters are assumed as additional coordinates). The parameters include thermal coefficients but also thermal sources and/or initial conditions. The resulting curse of dimensionality is circumvented by using the separated representation that implies a complexity that scales linearly with the dimension of the space.


INTRODUCTION
The complexity of coupled thermomechanical models defined in domains involving very fine space and time meshes induce computational issues. On the other hand, some times the parametric nature of such models (because of different uncertainties) induces additional difficulties. The use of explicit incremental methods is restricted to very short time steps because of the stability requirements. On the other hand, implicit incremental methods need the solution of a linear system at each time step in the general non-linear case. In any case, the treatment of parametric models requires the solutions of many direct models. To circumvent these difficulties one could use reduced approximation bases within the incremental framework. Another possibility lies in the use of an alternative non-incremental technique. This work explores the use of a particular non-incremental technique based on the separated representation of the unknown field.

SOLVING PROBLEMS WITH UNCERTAINTY
Let's consider the general heat equation: with the following initial condition: where u is the temperature field and k the thermal conductivity. In general the thermal conductivity depends on the temperature, but in what follows we are assuming that the conductivity is constant but unknown.

ON A SEPARATED REPRESENTATION
To solve multidimensional problems with a low calculation cost, an efficient method has been recently developed. This technique was successfully applied for treating some multi-dimensional models encountered in the kinetic theory description of complex fluids [1] [2] [3] and for treating high resolution thermal homogenization [5]. For the sake of simplicity the method is described here in the 1D transient case defined by: For that purpose let us consider the weak formulation of the heat equation: The main idea is the assumption that the solution can be written in the following separated form: The introduction of u 0 in the previous expression enables to enforce an homogeneous initial condition, i.e.: Building-up such a solution needs an iterative procedure. At the iteration n we assume that the function T i and X i are known i n ∀ ≤ . Then we are looking for: where R and S have to be determined. We assume the weighting field: Thus, the resulting weak formulation writes: that defines a non-linear problem which needs an iterative process to be solved. An alternate directions strategy has given excellent results in our former studies [1], [2], [6]. This method is performed in two steps: Step 1: R being known, we are looking for S.
Step 2: S being known, we are looking for R.
Starting with an arbitrary tentative function R, we perform step 1 and then step 2, and again both steps until reaching convergence. Only step 1 will be described now, step 2 being very similar. Thus, R is assumed known and then * 0 R = . It results: The integrals on t Ω can be computed numerically. We denote by i α , i β and γ these integrals, leading to: This is a 1D steady-state weak formulation (3D in the general case) which can be solved for instance with the finite elements method, but if we extract the strong form, then it could be solved using any discretization technique.

Introducing new coordinates
In this section, we consider a model in which there is an uncertainty on one or many parameters. For illustration purposes, we want to solve the heat equation where the thermal conductivity k is a badly known, but constant, parameter. A way to overcome this difficulty is to take k as a new coordinate belonging to the uncertainty interval. The transient solution for a particular value of the conductivity can then be computed by restricting the general solution to each particular value of this extracoordinate. Obviously, the price to be paid is the increase of the model dimensionality. However, this is not a serious issue when one proceeds within the separated representation framework just described. The solution is now assumed as: The iterative method described in the previous section is then used to compute the solution.

Results
We consider the heat problem in the domain ] [ 1,1 x Ω = − with, for the sake of simplicity, homogeneous essential boundary conditions. The initial condition is: 2 0 ( ) 1 u x x = − and we want to compute the temperature until t max =0.5s. and for 1 10 k ≤ ≤ .
The temperature field versus the thermal conductivity is depict in Figure 1 at time t=0.1s. and in Figure 2 at time t=0.5s. A reduced number of sums allow computing the solution with a reasonable accuracy. Obviously, the accuracy can be improved by using more sums in the separated approximation. In the limit case, the computed solution approaches the one computed by using the standard finite element method.

REDUCED SEPARATED REPRESENTATION
Though the separated representation allows considerable time savings, one could envisage additional CPU time savings. For instance the heat equation defined in large 2D or 3D domains needs the solution of large linear systems at step 1 of the iterative strategy presented above. We analyze the solution of the space step: where B i are the basis functions whose number is n x and s i are the weights associated with these basis functions. With this approximation, the virtual field yields: Introducing these notations in Eq. (11) it results: that leads to the linear system: If n x is smaller than the number of nodes distributed in the space, then significant CPU time savings can be accomplished.

Remark:
Such an approximation may also be performed on the time space for the function R(t).

BUILDING-UP THE APPROXIMATION
The main remaining question is how to build-up an optimal approximation basis. A solution is to build the basis functions by invoking the proper orthogonal decomposition as was performed in [4]. The main issue is the necessity to perform some previous simulations. We have opted for an "a priori" strategy. The idea is to check (at each iteration of the computation of S (respect. R)) if the basis needs to be enriched, and if that is the case, to propose an optimal enrichment.

Checking the accuracy of the basis
We assume that the basis functions are known for To check the accuracy of the basis, we compute the residual function that can be easily derived from Eq. : If the residual norm is greater than a given small enough value: 2 Res Res the basis must be enriched. On the contrary, if the norm of the residual is small enough, we continue to the next iteration keeping the same reduced basis.

Enrichment of the basis
The main aim is to define the enriched function B′ and the corresponding weight s′ in order to cancel the residual. From Eq. (17) we can write: that taking into account Eq. (17) leads to: which is an ordinary differential equation, that could be solved under the normality constraint 1 B′ = .
Obviously, the solution of that equation implies the solution of a linear system whose size corresponds to the number of nodes used in the discretization, however, this solution is only performed when the reduced basis needs for an enrichment, and not at each iteration as in the standard scheme. As soon as the new function B' is computed, it is added to the approximation basis, s' is added to the weights list and n x is increased of one unit. In the worst case, in which an enrichment must be performed at each iteration, the CPU time would be the same than the standard strategy. However, the analyzed cases prove that the CPU time can be reduced in one order of magnitude.

EXAMPLES
Let's consider the same example as in section 2.2.2 but with t max =1 and k varying linearly in time such as k=2t. Both, the time and the physical spaces are approximated employing 1000 nodes. The basis enrichment criterion in Eq. (18) is set to The computed solution is depicted in Figure 3.

Figure 3: Temperature versus time and position computed using a reduced separated representation
A reference solution u ref was computed using an implicit finite differences scheme. The finite difference scheme needed 120 seconds to achieve the results, whereas the standard separated representation required 13 seconds. The reduced separated representation only needed 7 seconds. The error (in the usual Lebesgue norm) is depicted Figure 4. In general the error is lower than 1% except in the boundaries neighbourhood where the solution being close to 0, the error in percent is artificially higher.

Figure 4: Error of the reduced separated representation method versus time and position
In this case the CPU time saving is not significant because the space domain is 1D, but the higher is the dimension of the space, the more significant will be the expected CPU time savings when one proceed within the reduced separated representation framework.

CONCLUSIONS
The separated representation has been successfully applied as a general solver for parametric linear and nonlinear transient models. A further reduction was accomplished by using reduced bases for accounting each model coordinate. This technique opens new possibilities to perform efficient simulations of complex thermomechanical models.