Problem formulation
We consider a domain \(\Omega \in \mathbb{R}^3\) with a boudary \(\Gamma\).
This domain is subjected to body force \(g\), boundary traction \(t\) and prescribed
displacement \(g\). Its material bheviour is isotropic and linear elastic. The problem can be
expressed as:
Find the displacements \(u : \Omega \Rightarrow \mathbb{R}^3\) such that:
\[\begin{split}\begin{align}
{\rm div}~\sigma + b = 0 & {\rm \quad in~} \Omega \\
u = g & {\rm \quad on~} \Gamma_D \\
\sigma \cdot n = t & {\rm \quad on~} \Gamma_N
\end{align}\end{split}\]
With \(\Gamma_D \cup \Gamma_N = \Gamma\) and \(n\) is the outside normal on the bounday \(\Gamma\).
The stress tensor \(\sigma\) depends on the strain tensor \(\varepsilon\) through the relation:
\[\sigma = D \varepsilon\]
where \(D\) is the constitutive matrix of the linear elatic behaviour.
Galerkin weak form gives the following equation after integration:
\[\int_{\Omega} \delta \varepsilon^T \sigma \, {\rm d} \Omega - \int_{\Omega} \delta u^T \, b \, {\rm d} \Omega - \int_{\Gamma} \delta u \, t \, {\rm d} \Gamma = 0\]
We set up a discrete formulation from the Galerkin weak form, given by:
\[K \, u = F\]
where \(K\) is the the stiffness matrix, \(u\) the discrete displacements defined at control points and \(F\) the load vector.
Stiffness matrix and force vector can be assembled from elementary values \(K^e\) and \(F^e\).
Elementary stiffness matrix is given by:
\[K^e_{kl} = \int_{\overline{\Omega^e}} B^T_k \, D \, B_l \left| J \right| \, {\rm d}\overline{\Omega}\]
where \(k\) and \(l\) are control point indices of the element, \(B_i\) is the strain-displacement matrix, \(\left| J \right|\) is the Jacobian determinant and \(\overline{\Omega^e}\) is the parametric domain of the element.
Elementary force vector is given by:
\[F^e_k = \int_{\overline{\Omega^e}} N_k \, b \, \left| J \right| \, {\rm d}\overline{\Omega} + \int_{\overline{\Gamma^e_N}} N_k \, t \left| J \right| \, {\rm d}\overline{\Gamma}\]
where \(\overline{\Gamma}\) is the parametric domain of the boundary.
Sensitivities
In order to solve an optimization problem aimed at minimizing an objective function, we want to use a gradient-based algorithm.
To do this, we seek to express the gradients analytically.
We define an objective function \(f\), function of the design variables \(x\) and the state variable \(u\):
\[f := g \left( x, u\left( x \right) \right)\]
where u is the discrete solution of the linear system \(Ku=F\).
The derivative of the objectif function with respect to the design variables can be written as:
\[\frac{\mathrm{d} f}{\mathrm{d} x_i} = \frac{\partial g}{\partial x_i} + \frac{\partial g}{\partial u} \cdot \frac{\mathrm{d} u}{\mathrm{d} x_i}\]
We define \(u^*\) as the solution of the adjoint problem \(K u^* = \frac{\partial g}{\partial u}\). The derivative of the objective function becomes:
\[\frac{\mathrm{d} f}{\mathrm{d} x_i} = \frac{\partial g}{\partial x_i} + u^* \cdot \left( \frac{\partial F}{\partial x_i} - \frac{\partial K}{\partial x_i} u \right)\]
From optimization model to analysis model
TO DO
Embedded formulation
Table 1 Notations
|
Embedded solid |
Hull |
|---|
basis functions |
\(R\) |
\(N\) |
parameter |
\(\theta\) |
\(\xi\) |
physical points |
\(\xi\) |
\(X\) |
control points |
\(P\) |
\(Q\) |
In the embedded solid, the coordinate \(\xi\) is expressed as:
\[\xi = \sum_{a} R_a \left( \theta \right) P_a\]
And its derivative with respect to the coordinate \(\theta\):
\[\frac{\partial \xi}{\partial \theta} = \sum_{a} \frac{\partial R_a}{\partial \theta} P_a\]
and for specific directions \(i,j\):
\[\left( \frac{\partial \xi}{\partial \theta} \right)_{ij} = \frac{\partial \xi_i}{\partial \theta_j}\]
This quantity is stored in variable dxidtheta(i,j)
Then, we derive this quantity with respect to the coordinates of a particular control point of the embedded entity \(P_a\):
\[\left( \frac{\partial}{\partial P_a} \left( \frac{\partial \xi}{\partial \theta} \right) \right)_{ijk} = \frac{\partial}{\partial P_{a_k}} \left( \frac{\partial \xi_i}{\partial \theta_j} \right) = \frac{\partial R_a}{\partial \theta_j} \delta_{ik}\]
This quantity is stored in variable DdxidthetaDP(i,j,k)
The derivative of \(\frac{\partial \xi}{\partial \theta}\) with respect tyo the coordinates of the hull’s control points equals zero since this quantity does not depend on the control points \(Q\).
In the hull, the coordinate in the physical space is linked to the parametric coordinate \(\xi\) by the relation:
\[X = \sum_{a} N_a \left( \xi \right) Q_a\]
And its derivative with respect to \(\xi\):
\[\frac{\partial X}{\partial \xi} = \sum_a \frac{\partial N_a}{\partial \xi} Q_a\]
And for specific directions \(i,j\):
\[\left( \frac{\partial X}{\partial \xi} \right)_{ij} = \frac{\partial X_i}{\partial \xi_j}\]
This quantity is stored in variable dxdxi(i,j)
Its derivative with respect to a specific hull control point \(Q_a\) reads:
\[\left( \frac{\partial}{\partial Q_a} \left( \frac{\partial X}{\partial \xi} \right) \right)_{ijk} = \frac{\partial}{\partial Q_{a_k}} \left( \frac{\partial X_i}{\partial \xi_j} \right) = \frac{\partial N_a}{\partial \xi_j} \delta_{ik}\]
This quantity is stored in variable DdxdxiDQ(i,j,k)
To express the derivative with respect to embedded entity control point \(P_a\), we need to express the NURBS composition:
\[\frac{\partial X}{\partial \xi} = \sum_b \frac{\partial N_b \left( \sum_a R_a P_a \right)}{\partial \xi} Q_b\]
The derivative with respect to \(P_a\) reads:
\[\left( \frac{\partial}{\partial P_a} \left( \frac{\partial X}{\partial \xi} \right)\right)_{ijk} = \frac{\partial}{\partial P_{a_k}} \left( \frac{\partial X_i}{\partial \xi_j} \right) = R_a \frac{\partial^2 X_i}{\partial \xi_j \partial \xi_k}\]
This quantity is stored in variable DdxdxiDP(i,j,k)
As a summary
\[\begin{split}\begin{array}{c|c}
\left( \frac{\partial}{\partial Q_a} \left( \frac{\partial \xi}{\partial \theta} \right) \right)_{ijk} = 0 & \left( \frac{\partial}{\partial Q_a} \left( \frac{\partial X}{\partial \xi} \right) \right)_{ijk} = \frac{\partial N_a}{\partial \xi_j} \delta_{ik} \\
\left( \frac{\partial}{\partial P_a} \left( \frac{\partial \xi}{\partial \theta} \right) \right)_{ijk} = \frac{\partial R_a}{\partial \theta_j} \delta_{ik} & \left( \frac{\partial}{\partial P_a} \left( \frac{\partial X}{\partial \xi} \right) \right)_{ijk} = R_a \frac{\partial^2 X_i}{\partial \xi_j \partial \xi_k}
\end{array}\end{split}\]
Derivative of inverse mappings
In this part, we express the derivative of inverse mapping \(\frac{\partial \xi}{\partial X}\) and \(\frac{\partial \theta}{\partial \xi}\)
with respect to a quantity named \(\Lambda\) which can be eitehr the coordinates of control points of the hull or the embedded entity.
We start with:
\[\frac{\partial \xi}{\partial X} \cdot \frac{\partial X}{\partial \xi} = I\]
Which can be derived as:
\[\frac{\partial}{\partial \Lambda} \left( \frac{\partial \xi}{\partial X} \right) \cdot \frac{\partial X}{\partial \xi} + \frac{\partial \xi}{\partial X} \cdot \frac{\partial}{\partial \Lambda} \left( \frac{\partial X}{\partial \xi} \right) = 0\]
Multiplying by \(\frac{\partial \xi}{\partial X}\) gives:
\[\frac{\partial}{\partial \Lambda} \left( \frac{\partial \xi}{\partial X} \right) = - \frac{\partial \xi}{\partial X} \cdot \frac{\partial}{\partial \Lambda} \left( \frac{\partial X}{\partial \xi} \right) \cdot \frac{\partial \xi}{\partial X}\]
The same reads for the derivative of \(\frac{\partial \theta}{\partial \xi}\):
\[\frac{\partial}{\partial \Lambda} \left( \frac{\partial \theta}{\partial \xi} \right) = - \frac{\partial \theta}{\partial \xi} \cdot \frac{\partial}{\partial \Lambda} \left( \frac{\partial \xi}{\partial \theta} \right) \cdot \frac{\partial \theta}{\partial \xi}\]
Derivative of the Jacobian determinant
- There are several transformation to take into account:
Reference element space \(\overline{\xi}\) to embedded entity parametric space \(\theta\)
Embedded entity parametric space \(\theta\) to hull parametric space \(\xi\)
Hull parametric space \(\xi\) to physical space \(X\)
\[J = \frac{\partial X}{\partial \overline{\xi}} = \frac{\partial X}{\partial \xi} \cdot \frac{\partial \xi}{\partial \theta} \cdot \frac{\partial \theta}{\partial \overline{\xi}}\]
\(\frac{\partial \theta}{\partial \overline{\xi}}\) does not depend on control points corrdinates. Thus, derivative of Jacobian determinant with respect to control points coordinates reads:
\[\frac{\partial \left| J \right|}{\partial \Lambda} = \frac{\partial}{\partial \Lambda} \left( \left| \frac{\partial X}{\partial \xi} \right| \right) \cdot \left| \frac{\partial \xi}{\partial \theta} \right| \cdot \left| \frac{\partial \theta}{\partial \overline{\xi}} \right| + \left| \frac{\partial X}{\partial \xi} \right| \cdot \frac{\partial}{\partial \Lambda} \left( \frac{\partial \xi}{\partial \theta} \right) \cdot \left| \frac{\partial \theta}{\partial \overline{\xi}} \right|\]
Jacobi’s formula give the expression of the derivative of a matrix determinant:
\[\frac{\mathrm{d}}{\mathrm{d} t} \mathrm{det} \left( A \left( t \right) \right) = \mathrm{det} \left( A \left( t \right) \right) \cdot \mathrm{tr} \left( A \left(t \right)^{-1} \cdot \frac{\mathrm{d} A \left( t \right)}{\mathrm{d} t}\right)\]
Applying Jocobi’s formula top our case gives:
\[\begin{split}\begin{eqnarray}
\frac{\partial}{\partial \Lambda} \left( \left| \frac{\partial X}{\partial \xi} \right|\right) & = & \left| \frac{\partial X}{\partial \xi} \right| \cdot \mathrm{tr} \left( \frac{\partial \xi}{\partial X} \cdot \frac{\partial}{\partial \Lambda} \left( \frac{\partial X}{\partial \xi} \right) \right) \\
\frac{\partial}{\partial \Lambda} \left( \left| \frac{\partial \xi}{\partial \theta} \right|\right) & = & \left| \frac{\partial \xi}{\partial \theta} \right| \cdot \mathrm{tr} \left( \frac{\partial \theta}{\partial \xi} \cdot \frac{\partial}{\partial \Lambda} \left( \frac{\partial \xi}{\partial \theta} \right) \right)
\end{eqnarray}\end{split}\]
And the derivative of \(\left| J \right|\) reads:
\[\begin{split}\begin{eqnarray}
\frac{\partial \left| J \right|}{\partial \Lambda} & = & \left( \left| \frac{\partial X}{\partial \xi} \right| \cdot \mathrm{tr} \left( \frac{\partial \xi}{\partial X} \cdot \frac{\partial}{\partial \Lambda} \left( \frac{\partial X}{\partial \xi} \right) \right) \right) \cdot \left| \frac{\partial \xi}{\partial \theta} \right| \cdot \left| \frac{\partial \theta}{\partial \overline{\xi}} \right| \\
&& + \left| \frac{\partial X}{\partial \xi} \right| \cdot \left( \left| \frac{\partial \xi}{\partial \theta} \right| \cdot \mathrm{tr} \left( \frac{\partial \theta}{\partial \xi} \cdot \frac{\partial}{\partial \Lambda} \left( \frac{\partial \xi}{\partial \theta} \right) \right) \right) \cdot \left| \frac{\partial \theta}{\partial \overline{\xi}} \right| \\
& = & \left| J \right| \cdot \left[ \mathrm{tr} \left( \frac{\partial \xi}{\partial X} \cdot \frac{\partial}{\partial \Lambda} \left( \frac{\partial X}{\partial \xi} \right) \right) + \mathrm{tr} \left( \frac{\partial \theta}{\partial \xi} \cdot \frac{\partial}{\partial \Lambda} \left( \frac{\partial \xi}{\partial \theta} \right) \right) \right]
\end{eqnarray}\end{split}\]
Applying this to the cases of control points \(P\) and \(Q\) gives:
\[\begin{split}\begin{eqnarray}
\frac{\partial \left| J \right|}{\partial P} & = & \left| J \right| \cdot \left[ \mathrm{tr} \left( \frac{\partial \xi}{\partial X} \cdot \frac{\partial}{\partial P} \left( \frac{\partial X}{\partial \xi} \right) \right) + \mathrm{tr} \left( \frac{\partial \theta}{\partial \xi} \cdot \frac{\partial}{\partial P} \left( \frac{\partial \xi}{\partial \theta} \right) \right) \right] \\
\frac{\partial \left| J \right|}{\partial Q} & = & \left| J \right| \cdot \mathrm{tr} \left( \frac{\partial \xi}{\partial X} \cdot \frac{\partial}{\partial Q} \left( \frac{\partial X}{\partial \xi} \right) \right)
\end{eqnarray}\end{split}\]