# Topology Optimization of Anisotropic Materials with Smooth Fiber Orientation
## 2 Model Derivation
### Density-Based Concurrent Optimization Framework  
The topology optimization framework proposed in this paper is an extension of the Solid Isotropic Material with Penalization (SIMP) method to orthotropic materials and combined with the CFAO method for concurrent optimization. Four-point element discretization is performed for a given design domain. For each element, element density (material distribution) and the fiber angle (element fiber orientation) are considered design variables. With compliance as the objective function, the volume fraction and fiber continuity rate are forced as constraints, and the mathematical optimization model is as follows:  
\[ \begin{split}&\min\,C(\mathbf{x},\theta)=\mathbf{U}^{\mathrm{T}} \cdot\mathbf{K}\cdot\mathbf{U}=\sum\limits_{e=1}^{N}x_{e}^{p}\mathbf{u}_{e}^{ \mathrm{T}}\cdot\mathbf{k}_{\theta}^{e}\cdot\mathbf{u}_{e}\\ &\text{s.t.}\ \ \mathbf{K}\cdot\mathbf{U}=\mathbf{F}\,\\ &\frac{V(\mathbf{x})}{V_{0}}\leq\text{volfrac}\,\\ & f(\theta)>0,\\ & x_{e}\in[0,1]\,\\ &\theta_{e}\in[-\frac{\pi}{2},\frac{\pi}{2}]\end{split}  \tag{1}\]where \(C\) represents the compliance, \(x\) is the density, \(\theta\) is the fiber angle, \(C\) is the function of \(x_{e}\) and \(\theta_{e}\), \(p\) is the penalty factor, \(\mathbf{U}\) and \(\mathbf{K}\) are the global displacement and global stiffness matrix, respectively, and \(x_{e}\), \(\mathbf{k}_{\theta}^{e}\), and \(\mathbf{u_{e}}\) represent the element density, element stiffness matrix, and node displacement matrix of element e. \(\mathbf{F}\) is the force vector, \(V(x)\) and \(V_{0}\) are the material volume and the design domain volume, respectively, volfrac is the predefined volume fraction, and \(f(\theta)\) is the fiber angle constraint function. In anisotropic materials, we calculate the compliance coefficient matrix in the natural coordinate system by introducing the transformation matrix \(\mathbf{T}\) as follows:  
\[\mathbf{D}=\mathbf{T}^{-1}\mathbf{C}\left(\mathbf{T}^{-1}\right)^{\mathsf{T}} \tag{2}\]  
where \(\mathbf{T}\) is the coordinate transformation matrix and \(\mathbf{C}\) is the compliance coefficient in the main coordinate system. \(\mathbf{C}\) and \(\mathbf{T}\) are calculated as follows:  
\[\mathbf{C}= \begin{bmatrix}\frac{E_{1}}{1-v_{12}v_{21}}&\frac{v_{21}E_{2}}{1-v _{12}v_{21}}&0\\ \frac{v_{21}E_{1}}{1-v_{12}v_{21}}&\frac{v_{21}}{1-v_{12}v_{21}}&0\\ 0&0&G_{12}\end{bmatrix}  \tag{3}\]  
\[\mathbf{T}= \begin{bmatrix}\cos^{2}\theta&\sin^{2}\theta&\sin 2\theta\\ \sin^{2}\theta&\cos^{2}\theta&-\sin 2\theta\\ -\cos\theta\sin\theta&\cos\theta\sin\theta&\cos^{2}\theta-\sin^{2}\theta\end{bmatrix}  \tag{4}\]  
Among them, fiber-reinforced composites are regarded as orthotropic materials, and there is \(\frac{v_{12}}{E_{2}}=\frac{v_{21}}{E_{1}}\). Therefore, in the plane stress problem, orthotropic materials have four independent elastic constants, \(E_{1},E_{2},v_{21}\), and \(G_{12}\), which represent the longitudinal elastic modulus, transverse elastic modulus, longitudinal Poisson's ratio, and shear modulus, respectively, and \(\theta\) is the fiber angle. Then, the element stiffness matrix of all elements is calculated as follows:  
\[\mathbf{k}_{\theta}=\iint\mathbf{B}^{\mathsf{T}}\mathbf{D}\mathbf{B}dxdy \tag{5}\]  
where \(\mathbf{k}_{\theta}\) a is the element stiffness matrix after counterclockwise rotation \(\mathsf{t}\), \(\mathbf{B}\) is the strain-stress matrix, and \(\mathbf{D}\) is the flexibility coefficient matrix in the natural coordinate system. Then, all element stiffness matrices are assembled into a global stiffness matrix \(\mathbf{K}\). According to the equilibrium equation \(\mathbf{F}=\mathbf{K}\mathbf{U}\) in Equation (1), the finite element analysis is carried out to solve the global displacement. Finally, a sensitivity analysis of design variables is performed.