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 through OpenPisco.PhysicalSolvers.SolverBase.SolverBase.SetAuxiliaryScalar() and OpenPisco.PhysicalSolvers.SolverBase.SolverBase.GetAuxiliaryScalar().

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:

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:

Zset#

Zset 3 is…To be described

OpenFOAM#

OpenFOAM 4 is an open-source finite volume solver developed by The OpenFOAM Foundation Ltd, for computational fluid dynamics (CFD). Its core implementation is in C++. It can handle many CFD cases, as noted in Ref. 5, OpenFOAM is capable of “the prediction of fluid motion and forces by computation using numerical analysis, generally extended to include heat, thermodynamics, chemistry and solids”.

Internal tools#

Muscat solver 6 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.UnstructuredFEAMecaSolver.UnstructuredFEAMecaProblem

Thermal (conduction) 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

  • OpenPisco.PhysicalSolvers.UnstructuredFEAThermalSolver.UnstructuredFEAThermalProblem

Thermal (conduction) analysis based on 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 fixed temperature

  • \(\rho\) is the material density

  • \(C_p\) is the thermal capacity

  • \(\lambda_c\) is the thermal conductivity

Associated implementation:

  • OpenPisco.PhysicalSolvers.GeneralAsterPhysicalSolver.AsterThermalEigenValues

Thermoelasticity (conduction) 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) &= {\cal A} e(u) + {\cal A} e^{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:

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

  • \(e^{ther}\) is the thermal strain equal to

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

  • \({\cal A}\) is the fourth order elasticity tensor equal to

\[\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

  • OpenPisco.PhysicalSolvers.UnstructuredFEAMecaSolver.UnstructuredFEAMecaThermalStrainProblem

Static buckling analysis#

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

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

Incompressible steady-state lamninar fluid flow analysis#

The Navier-Stokes equations represent the governing equations for incompressible steady-state laminar fluid flow, they are expressed as follows:

\[\begin{split}\begin{equation} \label{eq:steady_state_NS_incompressible} \left\{ \begin{array}{rllll} \text{div}(v) &=&0 &\text{in} &\Omega_f,\\ -\text{div}(\sigma_f(v,p)) + \rho\nabla vv &=& f_f &\text{in} &\Omega_f,\\ v &=& v_0 &\text{on} &{\partial}\Omega^D_f,\\ \sigma_f(v,p) n &=& 0 &\text{on} &{\partial}\Omega^N_f,\\ v &=& 0 &\text{on} &\Gamma,\\ \end{array} \right. \end{equation}\end{split}\]

where:

  • \(\rho\) is the material density

  • \(v\) is the veloctiy

  • \(p\) is the pressure

  • \(\sigma\) is the fluid stress tensor for a Newtonian fluid

\[\begin{equation} \sigma (v,p) = 2\mu e(v) - pI & with & e(v) = (\nable v + \nable v^T)/2 \end{equation}\]

  • \(\mu\) is the dynamic viscosity

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

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

Associated implementation:

  • OpenPisco.PhysicalSolvers.OpenFoamSteadyState

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

OpenFOAM

Linear elasticity analysis

True

True

True

True

False

Thermo linear elasticity analysis

True

False

False

False

False

Static buckling analysis

True

False

False

False

False

Modal analysis

True

False

False

False

False

Harmonic analysis

True

False

False

False

False

Thermal analysis

True

False

False

False

False

Thermal analysis eigenvalues

True

False

False

False

False

Incompressible steady-state laminar Newtonian flow

False

False

False

False

True

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.

strain rate#

The strain rate tensor is given by \(e(v)=(\nabla v + (\nabla v)^T)/2\) where \(v\) is a velocity.

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}.\]

potential_energy#

The potential 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_potential_energy#

It the potential 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

N,C,IP

potential_energy

N,C,IP

N

N

N,C,IP

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_potential_energy

True

True

True

True
1

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

2

https://freefem.org/

3

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

4

https://openfoam.org/

5

https://cfd.direct/openfoam/computational-fluid-dynamics/

6

https://gitlab.com/drti/muscat

7

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