Let us now apply the interface conditions to assemble two dynamic systems in various domains. For a broader overview, please refer to the paper “General Framework for Dynamic Substructuring” by our colleague Dennis de Klerk: link.
This section covers the implementation of the substructuring interface conditions for the physical domain. Let’s recall the governing equation (equilibrium between internal and external forces) in the physical domain and the coupling conditions:
$$ \mathbf{M}\ddot{\mathbf{u}} {(t)} + \mathbf{C}\dot{\mathbf{u}} {(t)} + \mathbf{K}{\mathbf{u}} {(t)} = \mathbf{f} {(t)} + \mathbf{g} {(t)} $$
$$ \mathbf{Bu}{(t)} = \mathbf{0} $$
$$ \mathbf{L}^\mathrm{T}\mathbf{g}{(t)} = \mathbf{0} $$
This includes the unknown interface forces $$ \mathbf{g}{(t) }$$ connecting the substructures. The time dependence $$ {(t)} $$ is omitted in the following equations for brevity. Depending on whether the motion $$ \mathbf{u}$$ or the forces $$ \mathbf{g}$$ are chosen as known at the interface, the primal or the dual-assembled system of equations is obtained, respectively:
Component Mode Synthesis (CMS) techniques are mainly used to reduce physical systems to describe a numerical model by its most dominant dynamic behavior. A benefit of CMS’s is that one can describe the physical domain’s internal DoFs by modes while the interface remains physical. Coupling is then performed on the physical DoFs like before.
Multiple reduction approaches $$ \mathbf{u} \approx \mathbf{R}\boldsymbol{\eta } $$ exist to reduce the physical DoFs $$\mathbf{u} $$ to a smaller set of generalized DoFs $$ \boldsymbol{\eta } $$. In the following paragraphs, some common model reduction methods are explained:
A basic approach is the Guyan reduction which performs static condensation. Let us consider a system without damping:
$$ \mathbf{M}\ddot{\mathrm{\mathbf{u}}} + \mathbf{Ku} = \mathbf{f} $$
The DoFs are divided into a set “b” on the boundary and an internal set “i”. If only external excitation on the boundary is assumed, this can be written as:
$$ \begin{bmatrix} \mathbf{M}_{\mathrm{bb}} & \mathbf{M}_{\mathrm{bi}}\\ \mathbf{M}_{\mathrm{ib}} & \mathbf{M}_{\mathrm{ii}} \end{bmatrix} \begin{bmatrix} \ddot{\mathbf{u}}_\mathrm{b}\\ \ddot{\mathbf{u}}_\mathrm{i} \end{bmatrix} \;\; + \;\; \begin{bmatrix} \mathbf{K}_{\mathrm{bb}} & \mathbf{K}_{\mathrm{bi}}\\ \mathbf{K}_{\mathrm{ib}} & \mathbf{K}_{\mathrm{ii}} \end{bmatrix} \begin{bmatrix} {\mathbf{u}}_\mathrm{b}\\ {\mathbf{u}}_\mathrm{i} \end{bmatrix} \;\; = \;\; \begin{bmatrix} {\mathbf{f}}\\ {\mathbf{0}} \end{bmatrix} $$
Furthermore, the internal set of response DoF $$ {\mathbf{u}}_\mathrm{i} $$ can be divided into a static and a dynamic part:
$$ {\mathbf{u}}_\mathrm{i} = {\mathbf{u}}_\mathrm{i}^{\mathrm{stat}} + {\mathbf{u}}_\mathrm{i}^{\mathrm{dyn}} $$
The static response can be exactly calculated by setting all accelerations $$ \ddot{\mathbf{u}} $$ to zero:
$$ \begin{bmatrix} \mathbf{K}_{\mathrm{bb}} & \mathbf{K}_{\mathrm{bi}}\\ \mathbf{K}_{\mathrm{ib}} & \mathbf{K}_{\mathrm{ii}} \end{bmatrix} \begin{bmatrix} {\mathbf{u}}_\mathrm{b}\\ {\mathbf{u}}_\mathrm{i} \end{bmatrix} \;\; = \;\; \begin{bmatrix} {\mathbf{f}}\\ {\mathbf{0}} \end{bmatrix} \rightarrow \mathbf{u}_\mathrm{i}^{\mathrm{stat}} = \underbrace {-\mathbf{K}_{\boldsymbol{\mathrm{ii}}}^{-1} \mathbf{K}_{\boldsymbol{\mathrm{ib}}}}_{\boldsymbol{\Phi}_\mathrm{c} }\mathbf{u}_\mathrm{b} $$
Finally, the Guyan reduction matrix $$ \mathbf{R}_\mathrm{G} $$ and the generalized coordinates $$ \boldsymbol{\eta}_\mathrm{G} $$ can be written as:
$$ \begin{bmatrix} \mathbf{u}_\mathrm{b}\\ \mathbf{u}_\mathrm{i} \end{bmatrix} = \begin{bmatrix} \mathbf{I}\\ \boldsymbol{\Phi}_\mathrm{c} \end{bmatrix} \mathbf{u}_\mathrm{b} = \mathbf{R}_\mathrm{G}\mathbf{u}_\mathrm{b} = \mathbf{R}_\mathrm{G} \boldsymbol{\eta}_\mathrm{G} $$
with the matrix $$ \boldsymbol{\Phi}_\mathrm{c} $$ representing the (static) interface constraint modes and the generalized DoF $$ \boldsymbol{\eta}_\mathrm{G} $$ containing only the boundary DoF $$ \mathbf{u}_\mathrm{b} $$.
The Hurty/Craig-Bampton method extends the Guyan reduction by adding several fixed-interface vibration modes representing the dynamic behavior of the internal dynamics. An advantage of this method is that this reduces interface locking when the reduction space $$ \mathbf{R} $$ is poor.
These fixed-interface vibration modes $$ \boldsymbol{\Phi}_\mathrm{i} $$ are computed as the undamped vibration modes of the subsystem constrained at all its boundary DoF $$ (\mathbf{u}_\mathrm{b} = \mathbf{0}) $$.
$$ (\mathbf{K}_{\mathrm{ii}} – \omega^\mathrm{2}_r \mathbf{M}_{\mathrm{ii}}) \boldsymbol{\varphi }_r = \mathbf{0} \;\;\;\; r=1,…,n $$
The eigenvectors $$ \boldsymbol{\varphi }_r $$ can be assembled in the fixed interface modal matrix $$ \boldsymbol{\Phi}_\mathrm{i} = [\boldsymbol{\varphi}_{\mathrm{i},1} \;\;\; \boldsymbol{\varphi}_{\mathrm{i},2} \;\;\; … \;\;\; \boldsymbol{\varphi}_{\mathrm{i},n}] $$. Using modal superposition, the internal dynamic responses can be approximated by:
$$ \mathbf{u}_\mathrm{i}^{\mathrm{dyn}} \approx \sum_{r=1}^{m} \boldsymbol{\varphi}_{\mathrm{i},r}\boldsymbol{\eta}_{\mathrm{i},r} = \boldsymbol{\Phi}_\mathrm{i}{\boldsymbol{\eta}}_{\mathrm{i}} $$
Where the internal dynamic response DoF $$ \mathbf{u}_\mathrm{i} $$ are reduced to internal generalized DoF $$ {\boldsymbol{\eta}}_{\mathrm{i}} $$. By limiting the displacements $$ \mathbf{u} $$ to the first $$ m\ll n $$ eigenmodes, one can drastically reduce the dimension of the original problem. The entire internal response $$ \mathbf{u}_\mathrm{i} $$ can be approximated by the addition of the static and dynamic internal responses:
$$ \mathbf{u}_\mathrm{i} \approx \boldsymbol{\Phi}_\mathrm{c}\mathbf{u}_\mathrm{b} + \boldsymbol{\Phi}_\mathrm{i}\boldsymbol{\eta} _\mathrm{i} $$
The reduction matrix $$ \mathbf{R}_{\mathrm{HCB}} $$ of the Hurty/Craig-Bampton method with its generalized coordinates $$ \boldsymbol{\eta} _{\mathrm{HCB}} $$ can then be expressed as:
$$ \begin{bmatrix} \mathbf{u}_\mathrm{b}\\ \mathbf{u}_\mathrm{i} \end{bmatrix} = \begin{bmatrix} \mathbf{I} & \mathbf{0}\\ \boldsymbol{\Phi}_\mathrm{c} & \boldsymbol{\Phi}_\mathrm{i} \end{bmatrix} \begin{bmatrix} \mathbf{u}_\mathrm{b}\\ \boldsymbol{\eta }_\mathrm{i} \end{bmatrix} = \mathbf{R}_{\mathrm{HCB}} \begin{bmatrix} \mathbf{u}_\mathrm{b}\\ \boldsymbol{\eta }_\mathrm{i} \end{bmatrix} = \mathbf{R}_{\mathrm{HCB}}\boldsymbol{\eta} _{\mathrm{HCB}} $$
With the matrix $$\boldsymbol{\Phi}_\mathrm{c }$$ representing the (static) interface constraint mode shapes and $$ \boldsymbol{\Phi}_\mathrm{i}$$ representing the (dynamic) fixed-interface vibration modes.
The Hintz-Herting reduction is similar to the Hurty/Craig-Bampton method, but it relies on the use of so-called inertia relief modes and residual flexibility modes as its reduced DoF (in addition to the boundary DoFs). The transformation matrix is given by:
$$ \mathbf{R}_{\mathrm{HH}} = \begin{bmatrix} \mathbf{I} & \mathbf{0} & \mathbf{0}\\ \boldsymbol{\Phi }_\mathrm{c} & \mathbf{H}_{\mathrm{ir}} & \boldsymbol{\Phi }_\mathrm{r} \end{bmatrix} $$
In this transformation matrix, $$ \boldsymbol{\Phi }_\mathrm{c} $$ again represents the (static) interface constraint mode shapes. $$ \mathbf{H}_{\mathrm{ir}} $$ represents the inertia relief modes and $$ \boldsymbol{\Phi }_\mathrm{r} $$ contains the residual flexibility modes.
Practical benefits:
This section focuses on the coupling of substructures whose dynamics $$\mathbf{u} $$ are approximated in generalized coordinates $$\boldsymbol{\eta}$$, with $$\mathbf{u} \approx \mathbf{R}\boldsymbol{\eta }$$. Let’s recall the governing equation in the reduced modal domain:
$$\mathbf{M}_\mathrm{m}\ddot{\boldsymbol{\eta} } + \mathbf{C}_\mathrm{m}\dot{\boldsymbol{\eta} }+ \mathbf{K}_\mathrm{m}{\boldsymbol{\eta} } = \mathbf{f}_\mathrm{m} + \mathbf{g}_\mathrm{m} \left\{\begin{matrix} \mathbf{M}_\mathrm{m} \overset{\Delta}{=} \mathbf{R}^\mathrm{T}\mathbf{MR}\\ \mathbf{C}_\mathrm{m} \overset{\Delta}{=} \mathbf{R}^\mathrm{T}\mathbf{CR}\\ \mathbf{K}_\mathrm{m} \overset{\Delta}{=} \mathbf{R}^\mathrm{T}\mathbf{KR}\\ \mathbf{f}_\mathrm{m} \overset{\Delta}{=} \mathbf{R}^\mathrm{T}\mathbf{f}\\ \mathbf{g}_\mathrm{m} \overset{\Delta}{=} \mathbf{R}^\mathrm{T}\mathbf{g} \end{matrix}\right. $$
Both the compatibility and force equilibrium equations need to be imposed on the subsystems in generalized coordinates $$\boldsymbol{\eta }$$. The compatibility condition can be transformed as:
$$\mathbf{Bu} = \mathbf{BR}\boldsymbol{\eta } = \mathbf{B}_\mathrm{m}\boldsymbol{\eta } = \mathbf{0}$$
The model reduction “weakens” the force equilibrium condition, i.e., the system response is computed only for forces that can be represented in the modal space:
$$\mathbf{L}_\mathrm{m}^\mathrm{T}\mathbf{g}_\mathrm{m} = \mathbf{0}$$
with the primal assembly operator being $$\mathbf{L}_\mathrm{m}^\mathrm{T} = \mathrm{null}(\mathbf{B}_\mathrm{m})$$. Again, this can be expressed in primal or dual form:
Frequency Based Substructuring (FBS) is characterized by assembling the FRFs of multiple substructures to obtain the FRFs of the entire system. Typically, the FRFs are obtained experimentally or from numerical models. Combining numerical and experimental models is referred to as “hybrid analysis”.
Let’s recall the dynamic equilibrium equation in the frequency domain:
$$[-\omega^\mathrm{2}\mathbf{M} + i\omega \mathbf{C} + \mathbf{K} ] \mathbf{u}(\omega ) = \mathbf{f}(\omega ) + \mathbf{g}(\omega )$$
$$ \mathbf{Z}(\omega)\mathbf{u}(\omega) = \mathbf{f}(\omega)+\mathbf{g}(\omega), \;\;\;\;\mathbf{u}(\omega)=\mathbf{Y}(\omega)(\mathbf{f}(\omega)+\mathbf{g}(\omega)) $$
This equation includes the unknown interface forces $$\mathbf{g}(\omega)$$ that couple the systems. To improve the readability, the frequency dependence $$(\omega)$$ is omitted in the following equations. For the two substructures A and B, the dynamic stiffness $$\mathbf{Z} $$ (or generally impedance) and receptance $$\mathbf{Y} $$ (or admittance in general) can be written in matrix notation:
$$\mathbf{Z}^\mathrm{A} \overset{\Delta}{=} \begin{bmatrix} \mathbf{Z}_{11}^\mathrm{A} & \mathbf{Z}_{12}^\mathrm{A}\\ \mathbf{Z}_{21}^\mathrm{A} & \mathbf{Z}_{22}^\mathrm{A} \end{bmatrix}, \;\;\;\;\; \mathbf{Y}^\mathrm{A} \overset{\Delta}{=} \begin{bmatrix} \mathbf{Y}_{11}^\mathrm{A} & \mathbf{Y}_{12}^\mathrm{A}\\ \mathbf{Y}_{21}^\mathrm{A} & \mathbf{Y}_{22}^\mathrm{A} \end{bmatrix}$$
$$\mathbf{Z}^\mathrm{B} \overset{\Delta}{=} \begin{bmatrix} \mathbf{Z}_{22}^\mathrm{B} & \mathbf{Z}_{23}^\mathrm{B}\\ \mathbf{Z}_{32}^\mathrm{B} & \mathbf{Z}_{33}^\mathrm{B} \end{bmatrix}, \;\;\;\;\; \mathbf{Y}^\mathrm{B} \overset{\Delta}{=} \begin{bmatrix} \mathbf{Y}_{22}^\mathrm{B} & \mathbf{Y}_{23}^\mathrm{B}\\ \mathbf{Y}_{32}^\mathrm{B} & \mathbf{Y}_{33}^\mathrm{B} \end{bmatrix}$$
To couple the subsystems, both matrices are put in block diagonal form:
$$\mathbf{Z} \overset{\Delta}{=} \begin{bmatrix} \mathbf{Z}^\mathrm{A} & \mathbf{0}\\ \mathbf{0} & \mathbf{Z}^\mathrm{B} \end{bmatrix}, \;\;\;\;\; \mathbf{Y} \overset{\Delta}{=} \begin{bmatrix} \mathbf{Y}^\mathrm{A} & \mathbf{0}\\ \mathbf{0} & \mathbf{Y}^\mathrm{B} \end{bmatrix}$$
Inserting this into the dynamic equation in impedance form, $$\mathbf{Z}$$ yields the starting point for the primal assembly, or, analogously in admittance notation, $$\mathbf{Y}$$ yielding the starting point for the dual assembly:
$$ \mathbf{Zu} = \mathbf{f} + \mathbf{g}, \;\;\; \mathbf{u} = \mathbf{Y}(\mathbf{f}+\mathbf{g}) $$
$$ \mathbf{Bu} = \mathbf{0} $$
$$\mathbf{L}^\mathrm{T}\mathbf{g}=\mathbf{0}$$
The coupling above is conducted “rigidly”, meaning that no gap is in between the substructures. Another coupling approach is compliant coupling. It implies that the force equilibrium condition between substructures still holds, but the compatibility condition is “weakened” by allowing a gap $${\Delta} \mathbf{u}$$ between the interface:
$$\mathbf{Bu} = \Delta \mathbf{u}$$
This can be represented by compliant damper elements (with stiffness $$\mathbf{K}_\mathrm{C}$$ and damping $$\mathbf{C}_\mathrm{C}$$) between substructures that carry no significant mass.
$$\mathbf{C}_\mathrm{C}\Delta \dot{\mathbf{u}} + \mathbf{K}_\mathrm{C}\Delta \mathbf{u} = \boldsymbol{\lambda} \rightarrow \mathbf{Z}_\mathrm{C}\Delta \mathbf{u}= \boldsymbol{\lambda}\rightarrow \Delta \mathbf{u} = \mathbf{Y}_\mathrm{C}\boldsymbol{\lambda}$$
In this case, the regular LM-FBS algorithm is augmented only for what concerns the interface term that gets inverted, namely in the following manner:
$$\boldsymbol{\lambda} = \underbrace{(\mathbf{BYB}^\mathrm{T} + \mathbf{Y}_\mathrm{C})^{-1}} _\text{interface stiffness} \;\;\;\;\;\; \underbrace{\mathbf{BYf}}_\text{incompatibility}$$
$$\mathbf{u}=\mathbf{Y}(\mathbf{f}-\mathbf{B}^\mathrm{T}(\mathbf{BYB}^\mathrm{T}+\mathbf{Y}_\mathrm{C})^{-1}\mathbf{BYf})$$
where $$ \mathbf{Y} $$ is a complete matrix containing the uncoupled admittance matrices on its diagonal and $$\mathbf{Y}_\mathrm{C}$$ is the (optional) interface flexibility. The procedure for compliant coupling of dynamic stiffness matrices is as follows:
$$\mathbf{f}=\mathbf{L}^\mathrm{T}(\mathbf{Z}+\mathbf{B}^\mathrm{T}\mathbf{Z}_\mathrm{C}\mathbf{B})\mathbf{Lu}$$
where, again, $$\mathbf{Z}$$ is a block-diagonal system matrix, $$\mathbf{Z}_\mathrm{C}$$ is the (optional) interface stiffness matrix.
A subclass of FBS methods is substructure decoupling. Decoupling is relevant for substructures that can only be measured in an assembly and when substructures were added for measurement purposes (e.g. Transmission Simulators).
For decoupling $$(\mathrm{AB – B=A})$$, a set of external forces $$ \mathbf{f} $$, motion $$ \mathbf{u} $$ and interface forces $$ \mathbf{g} $$ are defined:
$$\mathbf{u} \overset{\Delta}{=} \begin{bmatrix} \mathbf{u}_1^{\mathrm{AB}}\\ \mathbf{u}_2^{\mathrm{AB}}\\ \mathbf{u}_3^{\mathrm{AB}}\\ \mathbf{u}_2^{\mathrm{B}}\\ \mathbf{u}_3^{\mathrm{B}} \end{bmatrix}, \;\;\;\; \mathbf{f} \overset{\Delta}{=} \begin{bmatrix} \mathbf{f}_1^{\mathrm{AB}}\\ \mathbf{f}_2^{\mathrm{AB}}\\ \mathbf{f}_3^{\mathrm{AB}}\\ \mathbf{f}_2^{\mathrm{B}}\\ \mathbf{f}_3^{\mathrm{B}} \end{bmatrix}, \;\;\;\;\mathbf{g} \overset{\Delta}{=} \begin{bmatrix} \mathbf{g}_1^{\mathrm{AB}}\\ \mathbf{g}_2^{\mathrm{AB}}\\ \mathbf{g}_3^{\mathrm{AB}}\\ \mathbf{g}_2^{\mathrm{B}}\\ \mathbf{g}_3^{\mathrm{B}} \end{bmatrix}$$
Resulting from coupling, the impedance $$\mathbf{Z}^{\mathrm{AB}}$$ and admittance $$\mathbf{Y}^{\mathrm{AB}}$$ are given by:
$$\mathbf{Z}^{\mathrm{AB}} = \begin{bmatrix} \mathbf{Z}^\mathrm{A}_{11} &\mathbf{Z}^\mathrm{A}_{12} & \mathbf{0}\\ \mathbf{Z}^\mathrm{A}_{21}& \mathbf{Z}^\mathrm{A}_{22} + \mathbf{Z}^\mathrm{B}_{22} & \mathbf{Z}^\mathrm{B}_{23}\\ \mathbf{0} & \mathbf{Z}^\mathrm{B}_{32} & \mathbf{Z}^\mathrm{B}_{33} \end{bmatrix}, \;\;\;\;\; \mathbf{Y}^{\mathrm{AB}} = \begin{bmatrix} \mathbf{Y}^\mathrm{AB}_{11} &\mathbf{Y}^\mathrm{AB}_{12} & \mathbf{Y}^\mathrm{AB}_{13}\\ \mathbf{Y}^\mathrm{AB}_{21} &\mathbf{Y}^\mathrm{AB}_{22} & \mathbf{Y}^\mathrm{AB}_{23}\\ \mathbf{Y}^\mathrm{AB}_{31} &\mathbf{Y}^\mathrm{AB}_{32} & \mathbf{Y}^\mathrm{AB}_{33} \end{bmatrix} $$
The disassembly $$(\mathrm{AB – B})$$ equals the assembly of the system $$\mathrm{AB}$$ with the negative substructure $$\mathrm{B}$$ impedances/admittances:
$$\mathbf{Z} \overset{\Delta}{=} \begin{bmatrix} \mathbf{Z}^{\mathrm{AB}} & \mathbf{0}\\ \mathbf{0} & -\mathbf{Z}^\mathrm{B} \end{bmatrix}, \;\;\;\; \mathbf{Y} \overset{\Delta}{=} \begin{bmatrix} \mathbf{Y}^{\mathrm{AB}} & \mathbf{0}\\ \mathbf{0} & -\mathbf{Y}^\mathrm{B} \end{bmatrix}$$
Inserting this into the dynamic equation in impedance form $$\mathbf{Z}$$ yields the starting point for the primal disassembly, or analogously in admittance notation $$\mathbf{Y}$$ yielding the starting point for the dual disassembly:
$$\mathbf{Zu} = \mathbf{f}+\mathbf{g}, \;\;\; \mathbf{u}=\mathbf{Y}(\mathbf{f}+\mathbf{g})$$
Therefore, decoupling can be understood as coupling with negative system matrices. The only difference between coupling and decoupling is the definition of the global matrices:
$$\mathrm{Coupling:} \mathbf{Z} \overset{\Delta}{=} \begin{bmatrix} \mathbf{Z}^\mathrm{A} & \mathbf{0}\\ \mathbf{0} & \mathbf{Z}^\mathrm{B} \end{bmatrix}, \;\;\;\; \mathbf{Y} \overset{\Delta}{=} \begin{bmatrix} \mathbf{Y}^\mathrm{A} & \mathbf{0}\\ \mathbf{0} & \mathbf{Y}^\mathrm{B} \end{bmatrix}$$
$$\mathrm{Decoupling:} \mathbf{Z} \overset{\Delta}{=} \begin{bmatrix} \mathbf{Z}^\mathrm{AB} & \mathbf{0}\\ \mathbf{0} & -\mathbf{Z}^\mathrm{B} \end{bmatrix}, \;\;\;\; \mathbf{Y} \overset{\Delta}{=} \begin{bmatrix} \mathbf{Y}^\mathrm{AB} & \mathbf{0}\\ \mathbf{0} & -\mathbf{Y}^\mathrm{B} \end{bmatrix}$$
Substructuring can be applied to many domains if the interface dynamics satisfy the displacement compatibility and the force equilibrium. Depending on whether one chooses the motion $$ \mathbf{u} $$ or the forces $$ \mathbf{g} $$ as unknowns at the interface, a primal or dual assembly is obtained.
$$\mathbf{L}=\mathrm{null}(\mathbf{B}), \;\;\;\; \mathbf{B}^\mathrm{T}=\mathrm{null}(\mathbf{L}^\mathrm{T})$$
The primal approach relates to impedance addition, i.e., assembly in FE. It literally assembles the substructures by mapping the matching DoF to the same coordinate in a unique set of DoFs. From here on, matrix inversion and coordinate expansion may provide access to a force response solution in the full coordinate space.
Dual assembly finds a coupled solution in the same coordinate space as the unassembled system. To do so, additional Lagrange multipliers are introduced that represent the magnitude of the connecting forces. The dual approach employs subsystem admittances and is, therefore, a popular choice in the field of experimental substructuring. The frequency domain variant, commonly known as Lagrange Multiplier Frequency Based Substructuring (LM-FBS), is a highly versatile algorithm.
Analogous to coupling, the decoupling procedure can be performed in a primal and in a dual fashion. Decoupling only differs from coupling in the definition of the global system matrices.
Contact us for more VIBES
Contact our support team or call us on +31 85 822 50 49
Subscribe to our newsletter
Privacy statement. By submitting this form you agree to having your details registered in our database. At your request, we'll remove your details from our database immediately *
