Physical Solvers

Amongst the main components required in a topology optimization dedicated plateform, there is the so-called physical solver brick. Its purpose it to provide a numerical approximation of the solution to a partial differential equation corresponding to the problem at hand, given a set of boundary condition and on a given mesh, ergo a discretisation of a domain \(\Omega\).

The OpenPisco.PhysicalSolvers.SolverBase class aims at providing a base structure and a common interface to the physical solvers usable within the platform.

Auxiliary field/scalar

When solving a physical problem, at the very least, we expect the physical solver to compute the solution of the PDE at hand (be it a displacement, a temperature or something else). However, depending on the use case, we may want to rely on the physical solver to compute in the same time other physical quantities. This is the issue addressed by the auxiliary concept.

Note that the current solvers available in OpenPisco do not have the exact same capabilities, ergo some may be able to compute natively on their own (with or without a little help) a given physical quantity while the others may not.

The SolverBase class contains a dictionary named OpenPisco.PhysicalSolvers.SolverBase.SolverBase.auxiliaryFieldGeneration to store whether a solver is able to provide if required a given field evaluated on different support (namely on nodes, centroids or integration points for instance). In a derived class, such an attribute is redefined according to the solver capabilities.

Thus

• To check whether a solver is able to compute a named quantity on a given support, one can simply use OpenPisco.PhysicalSolvers.SolverBase.SolverBase.CanSupplyAuxiliaryField()

• Before solving the problem, one can ask the problem to generate a named quantity on a given support through OpenPisco.PhysicalSolvers.SolverBase.SolverBase.SetAuxiliaryFieldGeneration(). These quantities are to be computed when the problem is actually solved only.

• After Solving the physical problem, one can retrieve any named quantity on a given support after the problem has been solved but that were asked prior to the computation through OpenPisco.PhysicalSolvers.SolverBase.SolverBase.GetAuxiliaryField()

A similar logic is used for scalar quantities.

FieldsNames

All the fields and scalars mentionned are defined in OpenPisco.PhysicalSolvers.FieldsNames. The idea is to establish a mapping between each field/scalar variable and the actual field/scalar name (a string) used to denote a given field in physical solvers.

For instance, this is the name expected when one tries to retrieve a field/scalar from a file generated by a solver through the solver interface.

SolveByLevelSet

In the context of level set-based topology optimization, two settings are adressed:

• The conform case, or body-fitted approach, where the domain \(\Omega\) is simply bounded by the 0-th isusurface of the level set function. It can be extracted with OpenPisco.PhysicalSolvers.SolverBase.SolverBase.ExtractInteriorMesh()

• The fix mesh setting, or non conform, where we do not have an explicit mesh of the domain \(\Omega\) and therefore have to rely on the Ersatz material approximation (the default value for this pseudo-material stiffness is given in OpenPisco.PhysicalSolvers.SolverBase.SolverBase.eVoid)

Depending of the type of level set-based topology optimization performed, the mesh given to the physical solver (so-called OpenPisco.PhysicalSolvers.SolverBase.SolverBase.cleanMesh) must be defined accordingly in the derived classes. Note that it does not mean all physical analysis supported by the physical solvers available are supposed to handle both cases. The last step is to call the actual physical solver to solve the physical problem.

As of now, 4 interfaces with finite element solvers are available (three of them relies on external tools while the last is a native solver). Even though they do share a common interface, the core principles behind their respective implementation is different. The Macro & template files related to three of them are in the ExternalTools module.

Physical solvers available

External tools

Albeit slightly different regarding their core inner working, these solvers still share a common interface.

Code_Aster

Code_Aster 1 is an open-source finite element solver developed by EDF, for the analysis of structures and thermomechanics. Its core implementation is in Fortran and is also endowed as a domain specific language as a commande language. It can handle many behavior laws (elasticity, thermoelasticity, plasticity, fracture, hyperelasticity), finite elements (continuous 2D/3D, structural elements, kinematic couplings) and types of analysis (linear/nonlinear static, quasi-static, modal, dynamic on a modal basis or direct). It is also possible to natively compute quantities at Gauss points or at nodes.

The interfacing with OpenPisco is done with the following entities:

• The class OpenPisco.ExternalTools.Aster.AsterInterface.AsterInterface, whose main responsabilities is to

  • Select and retrieve the proper study among the template available based on the type of physical analysis considered, using OpenPisco.ExternalTools.Aster.AsterInterface.AsterInterface.GetStudyFilesFromPath(). The studies are all located in OpenPisco.ExternalTools.Aster where

    • The *.comm files are scripts using Code_Aster domain specific high-level language to define and solve a physical problem

    • The *.export are the actual file given to the executable. It contains input/output filenames required for/produce by the physical analysis, including the *.comm files ordered adequately for the analysis to run properly

  • Call the executable to solve the physical problem at hand with OpenPisco.ExternalTools.Aster.AsterInterface.AsterInterface.ExecuteAster()

  • Read the files generated by Code_Aster in order, for instance, to retrieve fields through OpenPisco.ExternalTools.Aster.AsterInterface.AsterInterface.GetField()

• Classes derived from OpenPisco.ExternalTools.Aster.AsterSolverBase.AsterSolverBase, for instance in OpenPisco.PhysicalSolvers.GeneralAsterPhysicalSolver, where each of them are associated to a specific physical analysis.

Their role is to:

• Generate a proper python script to be used by the studies (*.comm) with OpenPisco.ExternalTools.Aster.AsterSolverBase.AsterSolverBase.WriteParametersInput(). Amongst others, and depending on the type of analysis, it can be employed to provide information to the solver regarding

  • material properties

  • boundary regions

  • boundary conditions

• Specify the auxiliary fields/scalars to be generated with OpenPisco.ExternalTools.Aster.AsterSolverBase.AsterSolverBase.WriteAuxiliaryInputs()

• Solve the physical problem with OpenPisco.ExternalTools.Aster.AsterSolverBase.AsterSolverBase.Solve()

Freefem++

Freefem++ 2 is a high level finite element software using a C++ idiom domain specific language. It aims to ease the setting of a physical problem by relying on its variational formulation through this idiom.

The interfacing with OpenPisco consists in:

• The class OpenPisco.ExternalTools.FreeFem.FreeFemInterface.FreeFemInterface, whose main responsabilities is to

  • Call the executable to solve the physical problem at hand with OpenPisco.ExternalTools.FreeFem.FreeFemInterface.FreeFemInterface.ExecuteFreefem()

  • Read the files generated by FreeFem in order, for instance, to retrieve fields through OpenPisco.ExternalTools.FreeFem.FreeFemInterface.FreeFemInterface.GetNodalField()

• The class OpenPisco.ExternalTools.FreeFem.FreeFemProblemWriter.FreeFemProblemWriter, whose main responsabilities is to generate the Freefem scripts by combining the physical problem information in OpenPisco.ExternalTools.FreeFem.FreeFemProblemWriter.FreeFemProblemWriter.WriteFreefemProblem(), such as

  • material properties

  • boundary regions

  • boundary conditions

Zset

Zset 3 is…To be described

Internal tools

Muscat solver 4 is the internal finite element solver available within OpenPisco. It is implemented in Python, and Cython/C++ for performance critical routines and relies on the problem variational formulation. The interfacing with OpenPisco consist in a single module OpenPisco.PhysicalSolvers.UnstructuredFEAGenericProblem.

Note that, unlike the external tools:

• There is no need to performed any output operation like writing the mesh for instance, since OpenPisco already use Muscat routine to handle the mesh

• There is no executable involved

Physical analysis available

We provide now a list of the physical analysis supported by the platform, it is mean to be exhaustive. Let \(\Omega \subset \mathbb{R}^3\) such that \(\partial \Omega = \Gamma_D \cup \Gamma_N \cup \Gamma\).

Linear elasticity analysis

We denote by \(u\) the solution in displacement of the following problem

\[\begin{split}\begin{equation} \left\{ \begin{aligned} \text {div} (\sigma(u)) + f &= 0 \qquad &\text{ in } \Omega,\\ \sigma(u) \cdot n &= g \qquad &\text{ on } \Gamma_N, \\ \sigma(u) \cdot n &= 0 \qquad &\text{ on } \Gamma, \\ u &= 0\qquad &\text{ on } \Gamma_D.\\ \end{aligned} \right. \end{equation}\end{split}\]

where:

• \(\sigma\) is the stress tensor satisfying the Hooke’s law

\[\begin{equation} \sigma (u) = 2\mu e(u) + \lambda \text{ tr}(e(u))I. \end{equation}\]

• \(e(u)\) is the strain tensor in the small deformation framework

• \(\lambda\) and \(\mu\) are the Lamé coefficient and can also be expressed with respect to \(E\), the Young modulus, and \(\kappa\), the Poisson ratio,

\[\begin{equation*} \lambda=\frac{\kappa E}{(1+\kappa)(1-2\kappa)},\qquad \mu=\frac{E}{2(1+\kappa)}, \end{equation*}\]

• \(g\) is the density of surfaces forces acting on \(\Gamma_N\)

• \(f\) is the density of volumic forces acting on \(\Omega\)

• the body is clamped on \(\Gamma_D\)

Associated implementation:

• OpenPisco.PhysicalSolvers.GeneralAsterPhysicalSolver.AsterStatic

• OpenPisco.PhysicalSolvers.GeneralZSetSolverMeca.GeneralZSetPhysicalProblemMeca

• OpenPisco.PhysicalSolvers.StaticFreeFemSolverMeca.StaticFreeFemSolverMeca

• OpenPisco.PhysicalSolvers.UnstructuredFEAGenericSolver.UnstructuredFEAGenericProblem

Thermal analysis

We are looking for the temperature \(T\), solution of the following problem

\[\begin{split}\begin{equation} \left\{ \begin{aligned} - \rho C_p \frac{\partial T}{\partial t} - \text {div} (\lambda \nabla T) &= S \qquad &\text{ in } \Omega \times \mathbb{R}^{+},\\ \lambda \nabla T \cdot n &= -h (T-T_{ext}) \qquad &\text{ on } \Gamma_N \times \mathbb{R}^{+}, \\ \lambda \nabla T \cdot n &= 0 \qquad &\text{ on } \Gamma \times \mathbb{R}^{+}, \\ T &= T_d\qquad &\text{ on } \Gamma_D \times \mathbb{R}^{+},\\ T(t=0) &= T_{init}\qquad &\text{ in } \Omega.\\ \end{aligned} \right. \end{equation}\end{split}\]

where:

• \(S\) is the thermal source

• \(T_{ext}\) is the emperature of the fluid in contact with the structure

• \(T_d\) is the prescribed temperature

• \(h\) is the convection coefficient

• \(\rho\) is the material density

• \(C_p\) is the thermal capacity

• \(\lambda_c\) is the thermal conductivity

• \(T_{init}\) is the initial temperature

• the temperature is prescribed at \(T_d\) on \(\Gamma_D\)

Associated implementation:

• OpenPisco.PhysicalSolvers.GeneralAsterPhysicalSolver.AsterThermal

Thermal analysis eigenvalues

The unknown temperature \(T\) of following thermal problem

\[\begin{split}\begin{equation} \left\{ \begin{aligned} - \rho C_p \frac{\partial T}{\partial t} - \text {div} (\lambda \nabla T) &= S \qquad &\text{ in } \Omega \times \mathbb{R}^{+},\\ T &= T_d\qquad &\text{ on } \partial \Omega \times \mathbb{R}^{+},\\ T(t=0) &= T_{init}\qquad &\text{ in } \Omega.\\ \end{aligned} \right. \end{equation}\end{split}\]

can be computed using a spectral decomposition involving the following eigenvalue problem

\[\begin{split}\begin{equation} \left\{ \begin{aligned} \text {div} (\lambda_c \nabla T_i) &= \mu_i T_i \qquad &\text{ in } \Omega \times \mathbb{R}^{+},\\ T &= T_d\qquad &\text{ on } \partial \Omega \times \mathbb{R}^{+},\\ \end{aligned} \right. \end{equation}\end{split}\]

where:

• \(\mu_i\) is the i-th eigenvalue

• \(T_i\) is the i-th eigenmode

• \(T_d\) is the prescribed temperature

• \(\rho\) is the material density

• \(C_p\) is the thermal capacity

• \(\lambda_c\) is the thermal conductivity

• the temperature is prescribed at \(T_d\) on \(\partial \Omega\)

Associated implementation:

• OpenPisco.PhysicalSolvers.GeneralAsterPhysicalSolver.AsterThermalEigenValues

Thermo linear elasticity analysis

We investigate a weakly coupled thermoelastic problem. As such, the first step is to perform a classic thermal analysis to find the temperature \(T\) (see thermal analysis for more details). The second step is, for a given temperature field, to find the displacement governed by the following problem

\[\begin{split}\begin{equation} \left\{ \begin{aligned} \text {div} (\sigma(u,\Delta T)) + f &= 0 \qquad &\text{ in } \Omega,\\ \sigma(u,\Delta T) &= \sigma^{elas}(u) + \sigma^{ther}(\Delta T) \qquad &\text{ in } \Omega,\\ \sigma(u) \cdot n &= g \qquad &\text{ on } \Gamma_N, \\ \sigma(u) \cdot n &= 0 \qquad &\text{ on } \Gamma, \\ u &= 0\qquad &\text{ on } \Gamma_D.\\ \end{aligned} \right. \end{equation}\end{split}\]

where:

• \(\sigma^{elas}\) is the elastic contribution to the stress tensor, satisfying the Hooke’s law

\[\begin{equation} \sigma (u) = 2\mu e(u) + \lambda \text{ tr}(e(u))I. \end{equation}\]

• \(e(u)\) is the strain tensor

• \(\sigma^{ther}\) is the thermal contribution to the stress tensor, satisfying the following relation

\[\begin{equation} \sigma^{ther}(\Delta T) = -\alpha {\cal A} (\Delta T I). \end{equation}\]

• \({\cal A}\) is the fourth order elasticity tensor given by

\[\begin{equation} {\cal A}_{ijkh}=\lambda\,\delta_{ij}\delta_{kh}+\mu(\delta_{ik}\delta_{jh}+\delta_{ih}\delta_{jk}) \end{equation}\]

with \(\delta_{ij}\) the kronecker symbol.

• \(\lambda\) and \(\mu\) are the Lamé coefficients

• \(g\) is the density of surfaces forces

• \(f\) is the density of volumic forces

• the body is clamped on \(\Gamma_D\)

Associated implementation:

• OpenPisco.PhysicalSolvers.GeneralAsterPhysicalSolver.AsterThermoElastic

Static buckling analysis

For this analysis, we refer to the following theoretical documentation 5.

Modal analysis

We are looking for the eigenpulsations \(\omega_j > 0\) and its associated eigenmodes \(r_j\) such that

\[\begin{split}\begin{equation}\label{eq:systemdirectmodal} \left\{ \begin{array}{ccccc} -\text{div}(\sigma (r_j))&=&\omega_j^2\rho r_j&\text{in}&\Omega,\\ r_j&=&0&\text{on}&\Gamma_D,\\ \sigma (r_j) \cdot n&=&0&\text{on}&\Gamma,\\ \end{array} \right. \end{equation}\end{split}\]

where:

• \(\sigma\) is the stress tensor satisfying the Hooke’s law

• \(\rho\) is material the density

• \(j\) is the mode id

• the body is clamped on \(\Gamma_D\)

Associated implementation:

• OpenPisco.PhysicalSolvers.GeneralAsterPhysicalSolver.AsterModal

Harmonic analysis

We consider the harmonic response problem of a visco-elastic body depicted by the following problem

\[\begin{split}\begin{equation} \label{u} \left\{ \begin{array}{rlll} \rho \ddot{u} + c\dot{u} - \text{div}(\sigma (u)) &= & 0 & \text{in }\Omega \times \mathbb{R}^{+}, \\ u &= & 0 & \text{on }\Gamma_D\times \mathbb{R}^{+},\\ \sigma(u) \cdot n &= &g(t) & \text{on }\Gamma_N\times \mathbb{R}^{+},\\ \sigma(u) \cdot n &= &0 & \text{on }\Gamma_R \cup \Gamma_M \times\mathbb{R}^{+}. \end{array} \right. \end{equation}\end{split}\]

where:

• \(\sigma\) is the stress tensor satisfying the Hooke’s law

• \(\rho\) is the material density

• \(c\) is the damping operator

• \(g\) is the density of surfaces forces acting on \(\Gamma_N\)

• the body is clamped on \(\Gamma_D\)

Associated implementation:

• OpenPisco.PhysicalSolvers.GeneralAsterPhysicalSolver.AsterHarmonicForced

Physical analysis summary

Note that the following table only specify whether a specific analysis is available with OpenPisco when using a given physical solver. “False” does not mean it is not possible to actually run an analysis with a solver; at most it means the interface to do so is not available.

Table 2 Physical analysis supported w.r.t. tools

Analysis	Code_Aster	Freefem++	Zset	Muscat solver
Linear elasticity analysis	True	True	True	True
Thermo linear elasticity analysis	True	False	False	False
Static buckling analysis	True	False	False	False
Modal analysis	True	False	False	False
Harmonic analysis	True	False	False	False
Thermal analysis	True	False	False	False
Thermal analysis eigenvalues	True	False	False	False

Fields available

stress

We refer to the cauchy stress tensor, denoted by \(\boldsymbol{\sigma}\).

strain

In the small deformation framework, the strain tensor is given by \(e(u)=(\nabla u + (\nabla u)^T)/2\) where \(u\) is a displacement.

von_mises

\[\begin{align} \sigma_\text{v} &= \sqrt{\frac{3}{2} s_{ij}s_{ij}} \end{align}\]

where \(s_{ij}\) are the components of the Stress deviator tensor \(\boldsymbol{\sigma}^\text{dev}\) defined by

\(\boldsymbol{\sigma}^\text{dev} = \boldsymbol{\sigma} - \frac{\operatorname{tr}\left(\boldsymbol{\sigma}\right)}{3} \mathbf{I}\). where

\[\operatorname{tr}\left(\boldsymbol{\sigma}\right)=\sum_{n=1}^{3} \sigma_{ii}.\]

elastic_energy

The elastic energy is given by \(\frac{1}{2}\boldsymbol{\sigma}(u):e(u)\)

Scalars available

volume

The volume of the domain is given by \(\int_{\Omega}{ dx}\)

mass

The volume of the domain is given by \(\int_{\Omega}{\rho dx}\) where \(\rho\) is the material density.

int_elastic_energy

It is simply the elastic energy integrated on the whole domain, that is \(\int_{\Omega}{\frac{1}{2}\boldsymbol{\sigma}(u):e(u) dx}\)

Fields/Scalars summary

Similarly to the physical analysis summary, an empty cell is to be understood as “the interface with this solver can not generate this field/scalar, it is not available as of now”.

We now denote by N,C,IP (ergo nodes, centroids, integration points) the field support type available for a field with respect to each solver. If any of this notation exist in a cell in the table below, it means the physical solver can generate this particular field on this specific support.

Table 3 Fields available w.r.t. tools

Analysis	Code_Aster	Freefem++	Zset	Muscat
stress	N,C,IP	N
strain	N,C,IP
von_mises	N,C,IP	N
elastic_energy	N,C,IP	N	N	N

Next, we introduce the scalars in a similar way

Table 4 Scalars available w.r.t. tools

Analysis	Code_Aster	Freefem++	Zset	Muscat
volume	False	False	False	False
mass	False	False	False	False
int_elastic_energy	True	True	True	True

1

https://www.code-aster.org/

2

https://freefem.org/

3

http://www.zset-software.com/

4

https://gitlab.com/drti/muscat

5

https://www.code-aster.org/V2/doc/v14/en/man_r/r7/r7.05.01.pdf